Support Spherical Harmonics

README.md

@@ -2,24 +2,56 @@
 
 [![Build Status](https://travis-ci.org/MikaelSlevinsky/FastTransforms.jl.svg?branch=master)](https://travis-ci.org/MikaelSlevinsky/FastTransforms.jl) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://MikaelSlevinsky.github.io/FastTransforms.jl/stable) [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://MikaelSlevinsky.github.io/FastTransforms.jl/latest)
 
-The aim of this package is to provide new classes of fast transforms with low
-pre-computation. One approach is based on the use of asymptotic formulae to
-relate the transforms to a small number of fast Fourier transforms. Another
-approach is based on a Toeplitz-dot-Hankel decomposition of the matrix of
-connection coefficients. Both new classes of fast transforms do not
-require large pre-computation for fast execution and they are designed
-to work on expansions of functions with any degree of regularity.
-
-The Chebyshev—Jacobi transform and its inverse are implemented. This
-allows the fast conversion of Chebyshev expansion coefficients to Jacobi expansion
-coefficients and back.
+The aim of this package is to provide fast orthogonal polynomial transforms that are designed for expansions of functions with any degree of regularity. There are multiple approaches to the classical connection problem, though the user does not need to know the specifics.
+
+One approach is based on the use of asymptotic formulae to relate the transforms to a small number of fast Fourier transforms. Another approach is based on a Toeplitz-dot-Hankel decomposition of the matrix of connection coefficients. Alternatively, the matrix of connection coefficients may be decomposed hierarchically à la Fast Multipole Method.
+
+## The Chebyshev—Legendre Transform
+
+The Chebyshev—Legendre transform allows the fast conversion of Chebyshev expansion coefficients to Legendre expansion coefficients and back.
+
 ```julia
 julia> Pkg.add("FastTransforms")
 
 julia> using FastTransforms
 
 julia> c = rand(10001);
 
+julia> leg2cheb(c);
+
+julia> cheb2leg(c);
+
+julia> norm(cheb2leg(leg2cheb(c))-c)
+5.564168202018823e-13
+```
+
+The implementation separates pre-computation into a type of plan. This type is constructed with either `plan_leg2cheb` or `plan_cheb2leg`. Let's see how much faster it is if we pre-compute.
+
+```julia
+julia> p1 = plan_leg2cheb(c);
+
+julia> p2 = plan_cheb2leg(c);
+
+julia> @time leg2cheb(c);
+  0.082615 seconds (11.94 k allocations: 31.214 MiB, 6.75% gc time)
+
+julia> @time p1*c;
+  0.004297 seconds (6 allocations: 78.422 KiB)
+
+julia> @time cheb2leg(c);
+  0.110388 seconds (11.94 k allocations: 31.214 MiB, 8.16% gc time)
+
+julia> @time p2*c;
+  0.004500 seconds (6 allocations: 78.422 KiB)
+```
+
+## The Chebyshev—Jacobi Transform
+
+The Chebyshev—Jacobi transform allows the fast conversion of Chebyshev expansion coefficients to Jacobi expansion coefficients and back.
+
+```julia
+julia> c = rand(10001);
+
 julia> @time norm(icjt(cjt(c,0.1,-0.2),0.1,-0.2)-c,Inf)
   0.258390 seconds (431 allocations: 6.278 MB)
 1.4830359162942841e-12
@@ -43,15 +75,61 @@ is valid for the half-open square `(α,β) ∈ (-1/2,1/2]^2`. Therefore, the fas
 when the parameters are inside. If the parameters `(α,β)` are not exceptionally beyond the square,
 then increment/decrement operators are used with linear complexity (and linear conditioning) in the degree.
 
-The Padua transform and its inverse are also implemented thanks to
-[Michael Clarke](https://github.com/MikeAClarke). These are optimized methods
-designed for computing the bivariate Chebyshev coefficients by interpolating a
-bivariate function at the Padua points on `[-1,1]^2`.
+## The Padua Transform
+
+The Padua transform and its inverse are implemented thanks to [Michael Clarke](https://github.com/MikeAClarke). These are optimized methods designed for computing the bivariate Chebyshev coefficients by interpolating a bivariate function at the Padua points on `[-1,1]^2`.
+
+```julia
+julia> n = 200;
+
+julia> N = div((n+1)*(n+2),2);
+
+julia> v = rand(N); # The length of v is the number of Padua points
+
+julia> @time norm(ipaduatransform(paduatransform(v))-v)
+0.006571 seconds (846 allocations: 1.746 MiB)
+3.123637691861415e-14
+
+```
+
+## The Spherical Harmonic Transform
+
+Let `F` be a matrix of spherical harmonic expansion coefficients arranged by increasing order in absolute value, alternating between negative and positive orders. Then `sph2fourier` converts the representation into a bivariate Fourier series, and `fourier2sph` converts it back.
+```julia
+julia> F = sphrandn(Float64, 256, 256);
+
+julia> G = sph2fourier(F);
+
+julia> H = fourier2sph(G);
+
+julia> norm(F-H)
+4.950645831278297e-14
+
+julia> F = sphrandn(Float64, 1024, 1024);
+
+julia> G = sph2fourier(F; sketch = :none);
+Pre-computing...100%|███████████████████████████████████████████| Time: 0:00:04
+
+julia> H = fourier2sph(G; sketch = :none);
+Pre-computing...100%|███████████████████████████████████████████| Time: 0:00:04
+
+julia> norm(F-H)
+1.1510623098225283e-12
+
+```
+
+As with other fast transforms, `plan_sph2fourier` saves effort by caching the pre-computation. Be warned that for dimensions larger than `1,000`, this is no small feat!
 
 # References:
 
-   1.	N. Hale and A. Townsend. <a href="http://dx.doi.org/10.1137/130932223">A fast, simple, and stable Chebyshev—Legendre transform using and asymptotic formula</a>, SIAM J. Sci. Comput., 36:A148—A167, 2014.
+   [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.
+
+   [2]  N. Hale and A. Townsend. <a href="http://dx.doi.org/10.1137/130932223">A fast, simple, and stable Chebyshev—Legendre transform using and asymptotic formula</a>, *SIAM J. Sci. Comput.*, **36**:A148—A167, 2014.
+
+   [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]  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>, in press at *IMA J. Numer. Anal.*, 2017.
 
-   2.	R. M. Slevinsky. <a href="http://arxiv.org/abs/1602.02618">On the use of Hahn's asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform</a>, arXiv:1602.02618, 2016.
+   [5]  R. M. Slevinsky. <a href="https://arxiv.org/abs/1705.05448">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, arXiv:1705.05448, 2017.
 
-   3.	A. Townsend, M. Webb, and S. Olver. <a href="http://arxiv.org/abs/1604.07486">Fast polynomial transforms based on Toeplitz and Hankel matrices</a>, arXiv:1604.07486, 2016.
+   [6]  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.

REQUIRE

@@ -1,3 +1,6 @@
 julia 0.5
 ToeplitzMatrices 0.2
-Compat 0.17
+HierarchicalMatrices 0.0.1
+LowRankApprox 0.0.2
+ProgressMeter 0.3.4
+Compat 0.17

docs/make.jl

@@ -18,6 +18,6 @@ deploydocs(
     julia  = "0.5",
     osname = "linux",
     target = "build",
-    deps   = nothing,
+    deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
     make   = nothing
     )

docs/src/index.md

@@ -1,8 +1,37 @@
 # FastTransforms.jl Documentation
 
+## Introduction
+
+In numerical analysis, it is customary to expand a function in a basis:
+```math
+f(x) \sim \sum_{\ell=0}^{\infty} f_{\ell} \phi_{\ell}(x).
+```
+It may be necessary to transform our representation to one in a new basis, say, ``\{\psi_m(x)\}_{m\ge0}``:
+```math
+f(x) \sim \sum_{m=0}^{\infty} g_m \psi_m(x).
+```
+In many cases of interest, both representations are of finite length ``n`` and we seek a fast method (faster than ``\mathcal{O}(n^2)``) to transform the original coefficients ``f_{\ell}`` to the new coefficients ``g_m``.
+
+A similar problem arises when we wish to evaluate ``f`` at a set of points ``\{x_m\}_{m=0}^n``. We wish to transform coefficients of ``f`` to values at the set of points in fewer than ``\mathcal{O}(n^2)`` operations.
 
 ## Fast Transforms
 
+```@docs
+leg2cheb
+```
+
+```@docs
+cheb2leg
+```
+
+```@docs
+plan_leg2cheb
+```
+
+```@docs
+plan_cheb2leg
+```
+
 ```@docs
 cjt
 ```
@@ -39,6 +68,18 @@ plan_paduatransform!
 plan_ipaduatransform!
 ```
 
+```@docs
+sph2fourier
+```
+
+```@docs
+fourier2sph
+```
+
+```@docs
+plan_sph2fourier
+```
+
 ## Other Exported Methods
 
 ```@docs
@@ -49,6 +90,10 @@ gaunt
 paduapoints
 ```
 
+```@docs
+sphevaluate
+```
+
 ## Internal Methods
 
 ```@docs

src/FastTransforms.jl

@@ -1,17 +1,29 @@
 __precompile__()
 module FastTransforms
 
-using Base, ToeplitzMatrices, Compat
+using Base, ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter, Compat
 
-import Base: *
-import Base: view
+import Base: *, \, size, view
+import Base: getindex, setindex!, Factorization, length
+import Base.LinAlg: BlasFloat, BlasInt
+import HierarchicalMatrices: HierarchicalMatrix, unsafe_broadcasttimes!
+import HierarchicalMatrices: A_mul_B!, At_mul_B!, Ac_mul_B!
+import LowRankApprox: ColPerm
 
 export cjt, icjt, jjt, plan_cjt, plan_icjt
 export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac
+export normleg2cheb, cheb2normleg, normleg12cheb2, cheb22normleg1
+export plan_leg2cheb, plan_cheb2leg
+export plan_normleg2cheb, plan_cheb2normleg
+export plan_normleg12cheb2, plan_cheb22normleg1
 export gaunt
 export paduatransform, ipaduatransform, paduatransform!, ipaduatransform!, paduapoints
 export plan_paduatransform!, plan_ipaduatransform!
 
+export SlowSphericalHarmonicPlan, FastSphericalHarmonicPlan, ThinSphericalHarmonicPlan
+export sph2fourier, fourier2sph, plan_sph2fourier
+export sphones, sphzeros, sphrand, sphrandn, sphevaluate
+
 # Other module methods and constants:
 #export ChebyshevJacobiPlan, jac2cheb, cheb2jac
 #export sqrtpi, pochhammer, stirlingseries, stirlingremainder, Aratio, Cratio, Anαβ
@@ -42,12 +54,15 @@ include("cjt.jl")
 
 include("toeplitzhankel.jl")
 
-leg2cheb(x...)=th_leg2cheb(x...)
-cheb2leg(x...)=th_cheb2leg(x...)
+#leg2cheb(x...)=th_leg2cheb(x...)
+#cheb2leg(x...)=th_cheb2leg(x...)
 leg2chebu(x...)=th_leg2chebu(x...)
 ultra2ultra(x...)=th_ultra2ultra(x...)
 jac2jac(x...)=th_jac2jac(x...)
 
+include("hierarchical.jl")
+include("SphericalHarmonics/SphericalHarmonics.jl")
+
 include("gaunt.jl")