Merge branch 'master' into feat-spherical-isometries

Author: MikaelSlevinsky

.gitignore

@@ -1,5 +1,5 @@
 docs/build/
-docs/site/
+docs/src/generated
 deps/build.log
 deps/libfasttransforms.*
 .DS_Store

Project.toml

@@ -1,6 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.10.1"
+version = "0.10.2"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

README.md

@@ -1,6 +1,6 @@
 # FastTransforms.jl
 
-[![Build Status](https://github.com/JuliaApproximation/FastTransforms.jl/workflows/CI/badge.svg)](https://github.com/JuliaApproximation/FastTransforms.jl/actions?query=workflow%3ACI) [![Travis](https://travis-ci.org/JuliaApproximation/FastTransforms.jl.svg?branch=master)](https://travis-ci.org/JuliaApproximation/FastTransforms.jl) [![codecov](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/stable) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/dev)
+[![Build Status](https://github.com/JuliaApproximation/FastTransforms.jl/workflows/CI/badge.svg)](https://github.com/JuliaApproximation/FastTransforms.jl/actions?query=workflow%3ACI) [![Travis](https://travis-ci.com/JuliaApproximation/FastTransforms.jl.svg?branch=master)](https://travis-ci.com/JuliaApproximation/FastTransforms.jl) [![codecov](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/stable) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/dev)
 
 `FastTransforms.jl` allows the user to conveniently work with orthogonal polynomials with degrees well into the millions.
 
@@ -157,18 +157,8 @@ julia> @time norm(ipaduatransform(paduatransform(v)) - v)/norm(v)
 
 # References:
 
-   [1]  B. Alpert and V. Rokhlin. <a href="http://dx.doi.org/10.1137/0912009">A fast algorithm for the evaluation of Legendre expansions</a>, *SIAM J. Sci. Stat. Comput.*, **12**:158—179, 1991.
+   [1]  D. Ruiz—Antolín and A. Townsend. <a href="https://doi.org/10.1137/17M1134822">A nonuniform fast Fourier transform based on low rank approximation</a>, *SIAM J. Sci. Comput.*, **40**:A529–A547, 2018.
 
-   [2]  N. Hale and A. Townsend. <a href="http://dx.doi.org/10.1137/130932223">A fast, simple, and stable Chebyshev—Legendre transform using an asymptotic formula</a>, *SIAM J. Sci. Comput.*, **36**:A148—A167, 2014.
+   [2]  R. M. Slevinsky. <a href="https://doi.org/10.1016/j.acha.2017.11.001">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, *Appl. Comput. Harmon. Anal.*, **47**:585—606, 2019.
 
-   [3]  J. Keiner. <a href="http://dx.doi.org/10.1137/070703065">Computing with expansions in Gegenbauer polynomials</a>, *SIAM J. Sci. Comput.*, **31**:2151—2171, 2009.
-
-   [4]  D. Ruiz—Antolín and A. Townsend. <a href="https://arxiv.org/abs/1701.04492">A nonuniform fast Fourier transform based on low rank approximation</a>, arXiv:1701.04492, 2017.
-
-   [5]  R. M. Slevinsky. <a href="https://doi.org/10.1093/imanum/drw070">On the use of Hahn's asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform</a>, *IMA J. Numer. Anal.*, **38**:102—124, 2018.
-
-   [6]  R. M. Slevinsky. <a href="https://doi.org/10.1016/j.acha.2017.11.001">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, *Appl. Comput. Harmon. Anal.*, **47**:585—606, 2019.
-
-   [7]  R. M. Slevinsky, <a href="https://arxiv.org/abs/1711.07866">Conquering the pre-computation in two-dimensional harmonic polynomial transforms</a>, arXiv:1711.07866, 2017.
-
-   [8]  A. Townsend, M. Webb, and S. Olver. <a href="https://doi.org/10.1090/mcom/3277">Fast polynomial transforms based on Toeplitz and Hankel matrices</a>, in press at *Math. Comp.*, 2017.
+   [3]  R. M. Slevinsky, <a href="https://arxiv.org/abs/1711.07866">Conquering the pre-computation in two-dimensional harmonic polynomial transforms</a>, arXiv:1711.07866, 2017.

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 Documenter = "~0.24"

docs/make.jl

@@ -1,14 +1,41 @@
-using Documenter, FastTransforms
+using Documenter, FastTransforms, Literate
+
+const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
+const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")
+
+examples = [
+    "chebyshev.jl",
+    "disk.jl",
+    "nonlocaldiffusion.jl",
+    "padua.jl",
+    "sphere.jl",
+    "spinweighted.jl",
+    "triangle.jl",
+]
+
+for example in examples
+    example_filepath = joinpath(EXAMPLES_DIR, example)
+    Literate.markdown(example_filepath, OUTPUT_DIR; execute=true)
+end
 
 makedocs(
-			doctest = false,
-			format = Documenter.HTML(),
-			sitename = "FastTransforms.jl",
-			authors = "Richard Mikael Slevinsky",
-			pages = Any[
-					"Home" => "index.md"
-					]
-			)
+            doctest = false,
+            format = Documenter.HTML(),
+            sitename = "FastTransforms.jl",
+            authors = "Richard Mikael Slevinsky",
+            pages = Any[
+                    "Home" => "index.md",
+                    "Examples" => [
+                        "generated/chebyshev.md",
+                        "generated/disk.md",
+                        "generated/nonlocaldiffusion.md",
+                        "generated/padua.md",
+                        "generated/sphere.md",
+                        "generated/spinweighted.md",
+                        "generated/triangle.md",
+                        ],
+                    ]
+        )
 
 
 deploydocs(

examples/chebyshev.jl

@@ -1,52 +1,46 @@
-#############
+# # Chebyshev transform
 # This demonstrates the Chebyshev transform and inverse transform,
 # explaining precisely the normalization and points
-#############
 
 using FastTransforms
-
-# first kind points -> first kind polynomials
 n = 20
+
+# First kind points $\to$ first kind polynomials
 p_1 = chebyshevpoints(Float64, n, Val(1))
 f = exp.(p_1)
 f̌ = chebyshevtransform(f, Val(1))
-
 f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
 f̃(0.1) ≈ exp(0.1)
 
-# first kind polynomials -> first kind points
+# First kind polynomials $\to$ first kind points
 ichebyshevtransform(f̌, Val(1)) ≈ exp.(p_1)
 
-# second kind points -> first kind polynomials
+# Second kind points $\to$ first kind polynomials
 p_2 = chebyshevpoints(Float64, n, Val(2))
 f = exp.(p_2)
 f̌ = chebyshevtransform(f, Val(2))
-
 f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
 f̃(0.1) ≈ exp(0.1)
 
-# first kind polynomials -> second kind points
+# First kind polynomials $\to$ second kind points
 ichebyshevtransform(f̌, Val(2)) ≈ exp.(p_2)
 
-
-# first kind points -> second kind polynomials
-n = 20
+# First kind points $\to$ second kind polynomials
 p_1 = chebyshevpoints(Float64, n, Val(1))
 f = exp.(p_1)
 f̌ = chebyshevutransform(f, Val(1))
 f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * f̌
 f̃(0.1) ≈ exp(0.1)
 
-# second kind polynomials -> first kind points
+# Second kind polynomials $\to$ first kind points
 ichebyshevutransform(f̌, Val(1)) ≈ exp.(p_1)
 
-
-# second kind points -> second kind polynomials
+# Second kind points $\to$ second kind polynomials
 p_2 = chebyshevpoints(Float64, n, Val(2))[2:n-1]
 f = exp.(p_2)
 f̌ = chebyshevutransform(f, Val(2))
 f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * f̌
 f̃(0.1) ≈ exp(0.1)
 
-# second kind polynomials -> second kind points
+# Second kind polynomials $\to$ second kind points
 ichebyshevutransform(f̌, Val(2)) ≈ exp.(p_2)

examples/disk.jl

@@ -1,47 +1,56 @@
-#############
+# # Holomorphic integration on the unit disk
 # In this example, we explore integration of a harmonic function:
-#
-#   f(x,y) = (x^2-y^2+1)/[(x^2-y^2+1)^2+(2xy+1)^2],
-#
+# ```math
+#   f(x,y) = \frac{x^2-y^2+1}{(x^2-y^2+1)^2+(2xy+1)^2},
+# ```
 # over the unit disk. In this case, we know from complex analysis that the
-# integral of a holomorphic function is equal to π × f(0,0).
-# We analyze the function on an N×M tensor product grid defined by:
-#
-#   rₙ = cos[(n+1/2)π/2N], for 0 ≤ n < N, and
-#
-#   θₘ = 2π m/M, for 0 ≤ m < M;
-#
+# integral of a holomorphic function is equal to $\pi \times f(0,0)$.
+# We analyze the function on an $N\times M$ tensor product grid defined by:
+# ```math
+# \begin{aligned}
+# r_n & = \cos\left[(n+\tfrac{1}{2})\pi/2N\right],\quad{\rm for}\quad 0\le n < N,\quad{\rm and}\\
+# \theta_m & = 2\pi m/M,\quad{\rm for}\quad 0\le m < M;
+# \end{aligned}
+# ```
 # we convert the function samples to Chebyshev×Fourier coefficients using
 # `plan_disk_analysis`; and finally, we transform the Chebyshev×Fourier
-# coefficients to disk harmonic coefficients using `plan_disk2cxf`.
+# coefficients to Zernike polynomial coefficients using `plan_disk2cxf`.
 #
-# For the storage pattern of the arrays, please consult the documentation.
-#############
+# For the storage pattern of the arrays, please consult the
+# [documentation](https://MikaelSlevinsky.github.io/FastTransforms).
 
 using FastTransforms, LinearAlgebra
 
+# Our function $f$ on the disk:
 f = (x,y) -> (x^2-y^2+1)/((x^2-y^2+1)^2+(2x*y+1)^2)
 
+# The Zernike polynomial degree:
 N = 5
 M = 4N-3
 
+# The radial grid:
 r = [sinpi((N-n-0.5)/(2N)) for n in 0:N-1]
-θ = (0:M-1)*2/M # mod π.
+
+# The angular grid (mod $\pi$):
+θ = (0:M-1)*2/M
 
 # On the mapped tensor product grid, our function samples are:
 F = [f(r*cospi(θ), r*sinpi(θ)) for r in r, θ in θ]
 
+# We precompute a Zernike--Chebyshev×Fourier plan:
 P = plan_disk2cxf(F)
+
+# And an FFTW Chebyshev×Fourier analysis plan on the disk:
 PA = plan_disk_analysis(F)
 
 # Its Zernike coefficients are:
 U = P\(PA*F)
 
-# The Zernike coefficients are useful for integration. The integral of f(x,y)
-# over the disk should be π/2 by harmonicity. The coefficient of Z_0^0
-# multiplied by √π is:
+# The Zernike coefficients are useful for integration. The integral of $f(x,y)$
+# over the disk should be $\pi/2$ by harmonicity. The coefficient of $Z_0^0$
+# multiplied by `√π` is:
 U[1, 1]*sqrt(π)
 
-# Using an orthonormal basis, the integral of [f(x,y)]^2 over the disk is
+# Using an orthonormal basis, the integral of $[f(x,y)]^2$ over the disk is
 # approximately the square of the 2-norm of the coefficients:
 norm(U)^2

examples/nonlocaldiffusion.jl

@@ -1,3 +1,39 @@
+# # Nonlocal diffusion on $\mathbb{S}^2$
+# This example calculates the spectrum of the nonlocal diffusion operator:
+# ```math
+# \mathcal{L}_\delta u = \int_{\mathbb{S}^2} \rho_\delta(|\mathbf{x}-\mathbf{y}|)\left[u(\mathbf{x}) - u(\mathbf{y})\right] \,\mathrm{d}\Omega(\mathbf{y}),
+# ```
+# defined in Eq. (2) of
+#
+# R. M. Slevinsky, H. Montanelli, and Q. Du, [A spectral method for nonlocal diffusion operators on the sphere](https://doi.org/10.1016/j.jcp.2018.06.024), *J. Comp. Phys.*, **372**:893--911, 2018.
+#
+# In the above, $0<\delta<2$, $-1<\alpha<1$, and the kernel:
+# ```math
+# \rho_\delta(|\mathbf{x}-\mathbf{y}|) = \frac{4(1+\alpha)}{\pi \delta^{2+2\alpha}} \frac{\chi_{[0,\delta]}(|\mathbf{x}-\mathbf{y}|)}{|\mathbf{x}-\mathbf{y}|^{2-2\alpha}},
+# ```
+# where $\chi_I(\cdot)$ is the indicator function on the set $I$.
+#
+# This nonlocal operator is diagonalized by spherical harmonics:
+# ```math
+# \mathcal{L}_\delta Y_\ell^m(\mathbf{x}) = \lambda_\ell(\alpha, \delta) Y_\ell^m(\mathbf{x}),
+# ```
+# and its eigenfunctions are given by the generalized Funk--Hecke formula:
+# ```math
+# \lambda_\ell(\alpha, \delta) = \frac{(1+\alpha) 2^{2+\alpha}}{\delta^{2+2\alpha}}\int_{1-\delta^2/2}^1 \left[P_\ell(t)-1\right] (1-t)^{\alpha-1} \,\mathrm{d} t.
+# ```
+# In the paper, the authors use Clenshaw--Curtis quadrature and asymptotic evaluation of Legendre polynomials to achieve $\mathcal{O}(n^2\log n)$ complexity for the evaluation of the first $n$ eigenvalues. With a change of basis, this complexity can be reduced to $\mathcal{O}(n\log n)$.
+#
+# First, we represent:
+# ```math
+# P_n(t) - 1 = \sum_{j=0}^{n-1} \left[P_{j+1}(t) - P_j(t)\right] = -\sum_{j=0}^{n-1} (1-t) P_j^{(1,0)}(t).
+# ```
+# Then, we represent $P_j^{(1,0)}(t)$ with Jacobi polynomials $P_i^{(\alpha,0)}(t)$ and we integrate using [DLMF 18.9.16](https://dlmf.nist.gov/18.9.16):
+# ```math
+# \int_x^1 P_i^{(\alpha,0)}(t)(1-t)^\alpha\,\mathrm{d}t = \left\{ \begin{array}{cc} \frac{(1-x)^{\alpha+1}}{\alpha+1} & \mathrm{for~}i=0,\\ \frac{1}{2i}(1-x)^{\alpha+1}(1+x)P_{i-1}^{(\alpha+1,1)}(x), & \mathrm{for~}i>0.\end{array}\right.
+# ```
+# The code below implements this algorithm, making use of the Jacobi--Jacobi transform `plan_jac2jac`.
+# For numerical stability, the conversion from Jacobi polynomials $P_j^{(1,0)}(t)$ to $P_i^{(\alpha,0)}(t)$ is divided into conversion from $P_j^{(1,0)}(t)$ to $P_k^{(0,0)}(t)$, before conversion from $P_k^{(0,0)}(t)$ to $P_i^{(\alpha,0)}(t)$.
+
 using FastTransforms, LinearAlgebra
 
 function oprec!(n::Integer, v::AbstractVector, alpha::Real, delta2::Real)
@@ -13,19 +49,6 @@ function oprec!(n::Integer, v::AbstractVector, alpha::Real, delta2::Real)
     return v
 end
 
-"""
-This example calculates the spectrum of the nonlocal diffusion operator:
-
-```math
-ℒ_δ u = ∫_𝕊² ρ_δ(|𝐱-𝐲|)[u(𝐱) - u(𝐲)] dΩ(𝐲),
-```
-
-defined in Eq. (2) of
-
-    R. M. Slevinsky, H. Montanelli, and Q. Du, A spectral method for nonlocal diffusion operators on the sphere, J. Comp. Phys., 372:893--911, 2018.
-
-available at https://doi.org/10.1016/j.jcp.2018.06.024
-"""
 function evaluate_lambda(n::Integer, alpha::T, delta::T) where T
     delta2 = delta*delta
     scl = (1+alpha)*(2-delta2/2)
@@ -60,7 +83,11 @@ function evaluate_lambda(n::Integer, alpha::T, delta::T) where T
     return lambda
 end
 
-lambda = evaluate_lambda(1024, -0.5, 1.0)
-lambdabf = evaluate_lambda(1024, parse(BigFloat, "-0.5"), parse(BigFloat, "1.0"))
+# The spectrum in `Float64`:
+lambda = evaluate_lambda(10, -0.5, 1.0)
+
+# The spectrum in `BigFloat`:
+lambdabf = evaluate_lambda(10, parse(BigFloat, "-0.5"), parse(BigFloat, "1.0"))
 
+# The $\infty$-norm relative error:
 norm(lambda-lambdabf, Inf)/norm(lambda, Inf)

examples/padua.jl

@@ -1,16 +1,20 @@
-#############
+# # Padua transform
 # This demonstrates the Padua transform and inverse transform,
 # explaining precisely the normalization and points
-#############
 
 using FastTransforms
 
+# We define the Padua points and extract Cartesian components:
 N = 15
 pts = paduapoints(N)
-x = pts[:,1]; y = pts[:,2]
+x = pts[:,1]
+y = pts[:,2];
 
+# We take the Padua transform of the function:
 f = (x,y) -> exp(x + cos(y))
-f̌ = paduatransform(f.(x , y))
+f̌ = paduatransform(f.(x , y));
+
+# and use the coefficients to create an approximation to the function $f$:
 f̃ = (x,y) -> begin
     j = 1
     ret = 0.0
@@ -21,6 +25,8 @@ f̃ = (x,y) -> begin
     ret
 end
 
+# At a particular point, is the function well-approximated?
 f̃(0.1,0.2) ≈ f(0.1,0.2)
 
+# Does the inverse transform bring us back to the grid?
 ipaduatransform(f̌) ≈ f̃.(x,y)