start triangular harmonics

Author: MikaelSlevinsky

LICENSE.md

@@ -1,6 +1,8 @@
 The FastTransforms.jl package is licensed under the MIT "Expat" License:
 
-> Copyright (c) 2017: Richard Mikael Slevinsky.
+> Copyright (c) 2016-2017: Richard Mikael Slevinsky and other contributors:
+>
+> https://github.com/MikaelSlevinsky/FastTransforms.jl/graphs/contributors
 >
 > Permission is hereby granted, free of charge, to any person obtaining
 > a copy of this software and associated documentation files (the

README.md

@@ -166,7 +166,7 @@ As with other fast transforms, `plan_sph2fourier` saves effort by caching the pr
 
    [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.
+   [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.
 
    [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.
 

src/FastTransforms.jl

@@ -29,6 +29,10 @@ export SlowSphericalHarmonicPlan, FastSphericalHarmonicPlan, ThinSphericalHarmon
 export sph2fourier, fourier2sph, plan_sph2fourier
 export sphones, sphzeros, sphrand, sphrandn, sphevaluate
 
+export SlowTriangularHarmonicPlan
+export tri2cheb, cheb2tri, plan_tri2cheb
+export triones, trizeros, trirand, trirandn, trievaluate
+
 # Other module methods and constants:
 #export ChebyshevJacobiPlan, jac2cheb, cheb2jac
 #export sqrtpi, pochhammer, stirlingseries, stirlingremainder, Aratio, Cratio, Anαβ
@@ -69,6 +73,7 @@ jac2jac(x...)=th_jac2jac(x...)
 
 include("hierarchical.jl")
 include("SphericalHarmonics/SphericalHarmonics.jl")
+include("TriangularHarmonics/TriangularHarmonics.jl")
 
 include("gaunt.jl")
 

src/SphericalHarmonics/fastplan.jl

@@ -77,7 +77,7 @@ function Base.At_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
         At_mul_B_col_J!(Y, BF[J-1], B, 2J)
         At_mul_B_col_J!(Y, BF[J-1], B, 2J+1)
     end
-    zero_spurious_modes!(Y)
+    sph_zero_spurious_modes!(Y)
 end
 
 Base.Ac_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, FP, X)

src/SphericalHarmonics/slowplan.jl

@@ -177,7 +177,7 @@ function Base.At_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
         A_mul_B_col_J!!(Y, p1inv, B, J)
         A_mul_B_col_J!!(Y, p1inv, B, J+1)
     end
-    zero_spurious_modes!(At_mul_B!(RP, Y))
+    sph_zero_spurious_modes!(At_mul_B!(RP, Y))
 end
 
 Base.Ac_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, SP, X)

src/SphericalHarmonics/sphfunctions.jl

@@ -1,4 +1,4 @@
-function zero_spurious_modes!(A::AbstractMatrix)
+function sph_zero_spurious_modes!(A::AbstractMatrix)
     M, N = size(A)
     n = N÷2
     @inbounds for j = 1:n

src/SphericalHarmonics/thinplan.jl

@@ -148,7 +148,7 @@ function Base.At_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
         end
     end
 
-    zero_spurious_modes!(Y)
+    sph_zero_spurious_modes!(Y)
 end
 
 Base.Ac_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, TP, X)

src/TriangularHarmonics/TriangularHarmonics.jl

@@ -0,0 +1,24 @@
+@compat abstract type TriangularHarmonicPlan{T} end
+
+function *(P::TriangularHarmonicPlan, X::AbstractMatrix)
+    A_mul_B!(zero(X), P, X)
+end
+
+function \(P::TriangularHarmonicPlan, X::AbstractMatrix)
+    At_mul_B!(zero(X), P, X)
+end
+
+include("trifunctions.jl")
+include("slowplan.jl")
+
+function plan_tri2cheb(A::AbstractMatrix, α, β, γ; opts...)
+    M, N = size(A)
+    if M ≤ 1023
+        SlowTriangularHarmonicPlan(A, α, β, γ)
+    else
+        ThinTriangularHarmonicPlan(A, α, β, γ; opts...)
+    end
+end
+
+tri2cheb(A::AbstractMatrix, α, β, γ; opts...) = plan_tri2cheb(A, α, β, γ; opts...)*A
+cheb2tri(A::AbstractMatrix, α, β, γ; opts...) = plan_tri2cheb(A, α, β, γ; opts...)\A

src/TriangularHarmonics/slowplan.jl

@@ -0,0 +1,131 @@
+immutable Pnmp1toPlm{T} <: AbstractRotation{T}
+    rotations::Vector{Givens{T}}
+end
+
+function Pnmp1toPlm{T}(::Type{T}, n::Int, m::Int, α::T, β::T, γ::T)
+    G = Vector{Givens{T}}(n)
+    @inbounds for ℓ = 1:n
+        c = sqrt((2m+β+γ+2)/(ℓ+2m+β+γ+2)*(2ℓ+2m+α+β+γ+2)/(ℓ+2m+α+β+γ+2))
+        s = sqrt(ℓ/(ℓ+2m+β+γ+2)*(ℓ+α)/(ℓ+2m+α+β+γ+2))
+        G[n+1-ℓ] = Givens(ℓ, ℓ+1, c, s)
+    end
+    Pnmp1toPlm(G)
+end
+
+@inline length(P::Pnmp1toPlm) = length(P.rotations)
+@inline getindex(P::Pnmp1toPlm, i::Int) = P.rotations[i]
+
+function Base.A_mul_B!(C::Pnmp1toPlm, A::AbstractVecOrMat)
+    @inbounds for i = 1:length(C)
+        A_mul_B!(C.rotations[i], A)
+    end
+    A
+end
+
+function Base.A_mul_B!(A::AbstractMatrix, C::Pnmp1toPlm)
+    @inbounds for i = length(C):-1:1
+        A_mul_B!(A, C.rotations[i])
+    end
+    A
+end
+
+immutable TriRotationPlan{T} <: AbstractRotation{T}
+    layers::Vector{Pnmp1toPlm{T}}
+end
+
+function TriRotationPlan{T}(::Type{T}, n::Int, α::T, β::T, γ::T)
+    layers = Vector{Pnmp1toPlm{T}}(n)
+    @inbounds for m = 0:n-1
+        layers[m+1] = Pnmp1toPlm(T, n-m, m, α, β, γ)
+    end
+    TriRotationPlan(layers)
+end
+
+function Base.A_mul_B!(P::TriRotationPlan, A::AbstractMatrix)
+    M, N = size(A)
+    @inbounds for m = N-1:-1:1
+        layer = P.layers[m]
+        for ℓ = (m+1):N
+            @simd for i = 1:length(layer)
+                G = layer[i]
+                a1, a2 = A[G.i1,ℓ], A[G.i2,ℓ]
+                A[G.i1,ℓ] = G.c*a1 + G.s*a2
+                A[G.i2,ℓ] = G.c*a2 - G.s*a1
+            end
+        end
+    end
+    A
+end
+
+function Base.At_mul_B!(P::TriRotationPlan, A::AbstractMatrix)
+    M, N = size(A)
+    @inbounds for m = 1:N-1
+        layer = P.layers[m]
+        for ℓ = (m+1):N
+            @simd for i = length(layer):-1:1
+                G = layer[i]
+                a1, a2 = A[G.i1,ℓ], A[G.i2,ℓ]
+                A[G.i1,ℓ] = G.c*a1 - G.s*a2
+                A[G.i2,ℓ] = G.s*a1 + G.c*a2
+            end
+        end
+    end
+    A
+end
+
+Base.Ac_mul_B!(P::TriRotationPlan, A::AbstractMatrix) = At_mul_B!(P, A)
+
+
+immutable SlowTriangularHarmonicPlan{T} <: TriangularHarmonicPlan{T}
+    RP::TriRotationPlan{T}
+    p::NormalizedLegendreToChebyshevPlan{T}
+    pinv::ChebyshevToNormalizedLegendrePlan{T}
+    B::Matrix{T}
+end
+
+function SlowTriangularHarmonicPlan{T}(A::Matrix{T}, α, β, γ)
+    @assert β == γ == -half(T)
+    @assert α == zero(T)
+    M, N = size(A)
+    n = N
+    RP = TriRotationPlan(T, n-1, α, β, γ)
+    p = plan_normleg2cheb(A)
+    pinv = plan_cheb2normleg(A)
+    B = zeros(A)
+    SlowTriangularHarmonicPlan(RP, p, pinv, B)
+end
+
+function Base.A_mul_B!(Y::Matrix, SP::SlowTriangularHarmonicPlan, X::Matrix)
+    RP, p, B = SP.RP, SP.p, SP.B
+    copy!(B, X)
+    A_mul_B!(RP, B)
+    M, N = size(X)
+    for J = 1:N
+        A_mul_B_col_J!!(Y, p, B, J)
+    end
+    @inbounds for J = 1:N
+        nrm = sqrt(π/(2-δ(J-1,0)))
+        @simd for I = 1:M
+            Y[I,J] *= nrm
+        end
+    end
+    Y
+end
+
+function Base.At_mul_B!(Y::Matrix, SP::SlowTriangularHarmonicPlan, X::Matrix)
+    RP, pinv, B = SP.RP, SP.pinv, SP.B
+    copy!(B, X)
+    M, N = size(X)
+    @inbounds for J = 1:N
+        nrm = sqrt((2-δ(J-1,0))/π)
+        @simd for I = 1:M
+            B[I,J] *= nrm
+        end
+    end
+    for J = 1:N
+        A_mul_B_col_J!!(Y, pinv, B, J)
+    end
+    tri_zero_spurious_modes!(At_mul_B!(RP, Y))
+end
+
+Base.Ac_mul_B!(Y::Matrix, SP::SlowTriangularHarmonicPlan, X::Matrix) = At_mul_B!(Y, SP, X)

src/TriangularHarmonics/trifunctions.jl

@@ -0,0 +1,58 @@
+function tri_zero_spurious_modes!(A::AbstractMatrix)
+    M, N = size(A)
+    @inbounds for j = 2:N
+        @simd for i = M-j+2:M
+            A[i,j] = 0
+        end
+    end
+    A
+end
+
+function trirand{T}(::Type{T}, m::Int, n::Int)
+    A = zeros(T, m, n)
+    for j = 1:n
+        for i = 1:m+1-j
+            A[i,j] = rand(T)
+        end
+    end
+    A
+end
+
+function trirandn{T}(::Type{T}, m::Int, n::Int)
+    A = zeros(T, m, n)
+    for j = 1:n
+        for i = 1:m+1-j
+            A[i,j] = randn(T)
+        end
+    end
+    A
+end
+
+function triones{T}(::Type{T}, m::Int, n::Int)
+    A = zeros(T, m, n)
+    for j = 1:n
+        for i = 1:m+1-j
+            A[i,j] = one(T)
+        end
+    end
+    A
+end
+
+trizeros{T}(::Type{T}, m::Int, n::Int) = zeros(T, m, n)
+
+doc"""
+Pointwise evaluation of triangular harmonic:
+
+```math
+\tilde{P}_{\ell,m}^{(\alpha,\beta,\gamma)}(x,y).
+```
+"""
+trievaluate(x, y, L, M, α, β, γ) = trievaluate(x, L, M, α, β, γ)*trievaluate(x, y, M, β, γ)
+
+function trievaluate(x::Number, L::Integer, M::Integer, α::Number, β::Number, γ::Number)
+
+end
+
+function trievaluate(x::Number, y::Number, M::Integer, β::Number, γ::Number)
+
+end