clarify the different A_mul_B!'s, work on SHT API and docs

Author: MikaelSlevinsky

README.md

@@ -96,16 +96,16 @@ julia> @time norm(ipaduatransform(paduatransform(v))-v)
 
 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 = rand(Float64, 251, 501); FastTransforms.zero_spurious_modes!(F);
+julia> F = sphrandn(Float64, 256, 256);
 
 julia> G = sph2fourier(F);
 
 julia> H = fourier2sph(G);
 
 julia> norm(F-H)
-7.422366861016818e-14
+4.950645831278297e-14
 
-julia> F = rand(Float64, 1024, 2047); FastTransforms.zero_spurious_modes!(F);
+julia> F = sphrandn(Float64, 1024, 1024);
 
 julia> G = sph2fourier(F; sketch = :none);
 Pre-computing thin plan...100%|██████████████████████████████████████████████████| Time: 0:00:04
@@ -114,7 +114,7 @@ julia> H = fourier2sph(G; sketch = :none);
 Pre-computing thin plan...100%|██████████████████████████████████████████████████| Time: 0:00:04
 
 julia> norm(F-H)
-1.5062262753260893e-12
+1.1510623098225283e-12
 
 ```
 

docs/src/index.md

@@ -4,24 +4,15 @@
 
 In numerical analysis, it is customary to expand a function in a basis:
 ```math
-f(x) = \sum_{\ell = 0}^{n} f_{\ell} \phi_{\ell}(x).
+f(x) \sim \sum_{\ell=0}^{\infty} f_{\ell} \phi_{\ell}(x).
 ```
-Orthogonal polynomials are examples of convenient bases for sufficiently smooth functions. To perform some operation, it may be necessary to convert our representation to one in a new basis, say, ``\{\psi_m(x)\}_{m\ge0}``.
-
-If each function ``\phi_{\ell}`` is expanded in the basis ``\psi_m``, such as:
-```math
-\phi_{\ell}(x) \sim \sum_{m=0}^{\infty} c_{m,\ell}\psi_{m}(x),
-```
-then our original function has the alternative representation:
-```math
-f(x) \sim \sum_{m = 0}^{\infty} g_m \psi_m(x),
-```
-where ``g_m`` are defined by the sum:
+It may be more convenient to transform our representation to one in a new basis, say, ``\{\psi_m(x)\}_{m\ge0}``:
 ```math
-g_m = \sum_{\ell = 0}^{n} c_{m,\ell} f_{\ell}.
+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``.
 
-This is the classical connection problem. In many cases of interest, the both representations are finite and we seek a fast method (faster than ``\mathcal{O}(n^2)``) to transform the coefficients ``f_{\ell}`` to ``g_m``. These are the fast transforms.
+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
 
@@ -99,6 +90,10 @@ gaunt
 paduapoints
 ```
 
+```@docs
+sphevaluate
+```
+
 ## Internal Methods
 
 ```@docs

src/FastTransforms.jl

@@ -3,10 +3,11 @@ module FastTransforms
 
 using Base, ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter, Compat
 
-import Base: *, \, size, view, A_mul_B!, At_mul_B!, Ac_mul_B!
+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
@@ -21,6 +22,7 @@ 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

src/SphericalHarmonics/Butterfly.jl

@@ -166,7 +166,7 @@ function A_mul_B!{T}(y::AbstractVecOrMat{T}, A::IDPackedV{T}, P::ColumnPermutati
     k, n = size(A)
     At_mul_B!(P, x, jstart)
     copy!(y, istart, x, jstart, k)
-    HierarchicalMatrices.A_mul_B!(y, A.T, x, istart, jstart+k)
+    A_mul_B!(y, A.T, x, istart, jstart+k)
     A_mul_B!(P, x, jstart)
     y
 end
@@ -176,7 +176,7 @@ for f! in (:At_mul_B!, :Ac_mul_B!)
         function $f!{T}(y::AbstractVecOrMat{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int)
             k, n = size(A)
             copy!(y, istart, x, jstart, k)
-            HierarchicalMatrices.$f!(y, A.T, x, istart+k, jstart)
+            $f!(y, A.T, x, istart+k, jstart)
             A_mul_B!(P, y, istart)
             y
         end
@@ -185,9 +185,9 @@ end
 
 ### A_mul_B!, At_mul_B!, and  Ac_mul_B! for a Butterfly factorization.
 
-A_mul_B!{T}(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) = A_mul_B_col_J!(u, B, b, 1)
-At_mul_B!{T}(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) = At_mul_B_col_J!(u, B, b, 1)
-Ac_mul_B!{T}(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) = Ac_mul_B_col_J!(u, B, b, 1)
+Base.A_mul_B!{T}(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) = A_mul_B_col_J!(u, B, b, 1)
+Base.At_mul_B!{T}(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) = At_mul_B_col_J!(u, B, b, 1)
+Base.Ac_mul_B!{T}(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) = Ac_mul_B_col_J!(u, B, b, 1)
 
 function A_mul_B_col_J!{T}(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int)
     L = length(B.factors) - 1
@@ -237,7 +237,7 @@ function A_mul_B_col_J!{T}(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::I
     for i = 1:tL
         ml = mu+1
         mu += size(columns[i], 1)
-        HierarchicalMatrices.A_mul_B!(u, columns[i], temp1, ml+COLSHIFT, inds[i])
+        A_mul_B!(u, columns[i], temp1, ml+COLSHIFT, inds[i])
     end
 
     u
@@ -267,7 +267,7 @@ for f! in (:At_mul_B!,:Ac_mul_B!)
             for i = 1:tL
                 ml = mu+1
                 mu += size(columns[i], 1)
-                HierarchicalMatrices.$f!(temp1, columns[i], b, inds[i], ml+COLSHIFT)
+                $f!(temp1, columns[i], b, inds[i], ml+COLSHIFT)
             end
 
             ii, jj = tL, 1

src/SphericalHarmonics/SphericalHarmonics.jl

@@ -1,15 +1,3 @@
-function zero_spurious_modes!(A::AbstractMatrix)
-    M, N = size(A)
-    n = N÷2
-    @inbounds for j = 1:n
-        @simd for i = M-j+1:M
-            A[i,2j] = 0
-            A[i,2j+1] = 0
-        end
-    end
-    A
-end
-
 @compat abstract type SphericalHarmonicPlan{T} end
 
 function *(P::SphericalHarmonicPlan, X::AbstractMatrix)
@@ -20,6 +8,7 @@ function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
     At_mul_B!(zero(X), P, X)
 end
 
+include("sphfunctions.jl")
 include("slowplan.jl")
 include("Butterfly.jl")
 include("fastplan.jl")

src/SphericalHarmonics/fastplan.jl

@@ -38,7 +38,7 @@ end
 
 FastSphericalHarmonicPlan(A::Matrix; opts...) = FastSphericalHarmonicPlan(A, round(Int, log2(size(A, 1)+1)-6); opts...)
 
-function A_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
+function Base.A_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
     RP, BF, p1, p2, B = FP.RP, FP.BF, FP.p1, FP.p2, FP.B
     fill!(B, zero(eltype(B)))
     M, N = size(X)
@@ -60,7 +60,7 @@ function A_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
     Y
 end
 
-function At_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
+function Base.At_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
     RP, BF, p1inv, p2inv, B = FP.RP, FP.BF, FP.p1inv, FP.p2inv, FP.B
     copy!(B, X)
     M, N = size(X)
@@ -82,6 +82,6 @@ function At_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
     zero_spurious_modes!(Y)
 end
 
-Ac_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, FP, X)
+Base.Ac_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, FP, X)
 
 allranks(FP::FastSphericalHarmonicPlan) = mapreduce(allranks,vcat,FP.BF)

src/SphericalHarmonics/slowplan.jl

@@ -1,5 +1,7 @@
 import Base.LinAlg: Givens, AbstractRotation
 
+### These three A_mul_B! should be in Base, but for the time being they do not add methods to Base.A_mul_B!; they add methods to the internal A_mul_B!.
+
 function A_mul_B!{T<:Real}(G::Givens{T}, A::AbstractVecOrMat)
     m, n = size(A, 1), size(A, 2)
     if G.i2 > m
@@ -56,14 +58,14 @@ end
 @inline length(P::Pnmp2toPlm) = length(P.rotations)
 @inline getindex(P::Pnmp2toPlm, i::Int) = P.rotations[i]
 
-function A_mul_B!(C::Pnmp2toPlm,A::AbstractVecOrMat)
+function Base.A_mul_B!(C::Pnmp2toPlm,A::AbstractVecOrMat)
     @inbounds for i = 1:length(C)
         A_mul_B!(C.rotations[i], A)
     end
     A
 end
 
-function A_mul_B!(A::AbstractMatrix, C::Pnmp2toPlm)
+function Base.A_mul_B!(A::AbstractMatrix, C::Pnmp2toPlm)
     @inbounds for i = length(C):-1:1
         A_mul_B!(A, C.rotations[i])
     end
@@ -83,7 +85,7 @@ function RotationPlan{T}(::Type{T}, n::Int)
     RotationPlan(layers)
 end
 
-function A_mul_B!(P::RotationPlan, A::AbstractMatrix)
+function Base.A_mul_B!(P::RotationPlan, A::AbstractMatrix)
     M, N = size(A)
     @inbounds for m = N÷2-2:-1:0
         layer = P.layers[m+1]
@@ -102,7 +104,7 @@ function A_mul_B!(P::RotationPlan, A::AbstractMatrix)
     A
 end
 
-function At_mul_B!(P::RotationPlan, A::AbstractMatrix)
+function Base.At_mul_B!(P::RotationPlan, A::AbstractMatrix)
     M, N = size(A)
     @inbounds for m = 0:N÷2-2
         layer = P.layers[m+1]
@@ -121,7 +123,7 @@ function At_mul_B!(P::RotationPlan, A::AbstractMatrix)
     A
 end
 
-Ac_mul_B!(P::RotationPlan, A::AbstractMatrix) = At_mul_B!(P, A)
+Base.Ac_mul_B!(P::RotationPlan, A::AbstractMatrix) = At_mul_B!(P, A)
 
 
 immutable SlowSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
@@ -145,7 +147,7 @@ function SlowSphericalHarmonicPlan{T}(A::Matrix{T})
     SlowSphericalHarmonicPlan(RP, p1, p2, p1inv, p2inv, B)
 end
 
-function A_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
+function Base.A_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
     RP, p1, p2, B = SP.RP, SP.p1, SP.p2, SP.B
     copy!(B, X)
     A_mul_B!(RP, B)
@@ -162,7 +164,7 @@ function A_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
     Y
 end
 
-function At_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
+function Base.At_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
     RP, p1inv, p2inv, B = SP.RP, SP.p1inv, SP.p2inv, SP.B
     copy!(B, X)
     M, N = size(X)
@@ -178,4 +180,4 @@ function At_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
     zero_spurious_modes!(At_mul_B!(RP, Y))
 end
 
-Ac_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, SP, X)
+Base.Ac_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, SP, X)

test/sphericalharmonictestfunctions.jl renamed to src/SphericalHarmonics/sphfunctions.jl

@@ -1,4 +1,16 @@
-function sphrand{T}(::Type{T}, m, n)
+function zero_spurious_modes!(A::AbstractMatrix)
+    M, N = size(A)
+    n = N÷2
+    @inbounds for j = 1:n
+        @simd for i = M-j+1:M
+            A[i,2j] = 0
+            A[i,2j+1] = 0
+        end
+    end
+    A
+end
+
+function sphrand{T}(::Type{T}, m::Int, n::Int)
     A = zeros(T, m, 2n-1)
     for i = 1:m
         A[i,1] = rand(T)
@@ -12,7 +24,7 @@ function sphrand{T}(::Type{T}, m, n)
     A
 end
 
-function sphrandn{T}(::Type{T}, m, n)
+function sphrandn{T}(::Type{T}, m::Int, n::Int)
     A = zeros(T, m, 2n-1)
     for i = 1:m
         A[i,1] = randn(T)
@@ -26,6 +38,22 @@ function sphrandn{T}(::Type{T}, m, n)
     A
 end
 
+function sphones{T}(::Type{T}, m::Int, n::Int)
+    A = zeros(T, m, 2n-1)
+    for i = 1:m
+        A[i,1] = one(T)
+    end
+    for j = 1:n
+        for i = 1:m-j
+            A[i,2j] = one(T)
+            A[i,2j+1] = one(T)
+        end
+    end
+    A
+end
+
+sphzeros{T}(::Type{T}, m::Int, n::Int) = zeros(T, m, 2n-1)
+
 function normalizecolumns!(A::AbstractMatrix)
     m, n = size(A)
     @inbounds for j = 1:n
@@ -54,7 +82,14 @@ function maxcolnorm(A::AbstractMatrix)
     norm(nrm, Inf)
 end
 
-function sphevaluatepi(θ::Number,L::Integer,M::Integer)
+doc"""
+Pointwise evaluation of spherical harmonic ``Y_{\ell}^m(\theta,\varphi)``.
+"""
+sphevaluate(θ, φ, L, M) = sphevaluatepi(θ/π, φ/π, L, M)
+
+sphevaluatepi(θ::Number, φ::Number, L::Integer, M::Integer) = sphevaluatepi(θ,L,M)*sphevaluatepi(φ,M)
+
+function sphevaluatepi(θ::Number, L::Integer, M::Integer)
     ret = one(θ)/sqrt(two(θ))
     if M < 0 M = -M end
     c, s = cospi(θ), sinpi(θ)
@@ -76,3 +111,5 @@ function sphevaluatepi(θ::Number,L::Integer,M::Integer)
         return ret
     end
 end
+
+sphevaluatepi(φ::Number, M::Integer) = complex(cospi(M*φ),sinpi(M*φ))/sqrt(two(φ)*π)

src/SphericalHarmonics/thinplan.jl

@@ -42,7 +42,7 @@ end
 
 ThinSphericalHarmonicPlan(A::Matrix; opts...) = ThinSphericalHarmonicPlan(A, round(Int, log2(size(A, 1)+1)-6); opts...)
 
-function A_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
+function Base.A_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
     RP, BF, p1, p2, B = TP.RP, TP.BF, TP.p1, TP.p2, TP.B
     copy!(B, X)
     M, N = size(X)
@@ -98,7 +98,7 @@ function A_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
 end
 
 
-function At_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
+function Base.At_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
     RP, BF, p1inv, p2inv, B = TP.RP, TP.BF, TP.p1inv, TP.p2inv, TP.B
     copy!(B, X)
     M, N = size(X)
@@ -153,7 +153,7 @@ function At_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
     zero_spurious_modes!(Y)
 end
 
-Ac_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, TP, X)
+Base.Ac_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, TP, X)
 
 allranks(TP::ThinSphericalHarmonicPlan) = mapreduce(i->allranks(TP.BF[i]),vcat,sort!([LAYERSKELETON-1:LAYERSKELETON:length(TP.BF);LAYERSKELETON:LAYERSKELETON:length(TP.BF)]))