Support ApproxFun's (unfortunate...) reverse lexigraphical order

Author: dlfivefifty

src/PaduaTransform.jl

@@ -1,22 +1,29 @@
 """
 Pre-plan an Inverse Padua Transform.
 """
-immutable IPaduaTransformPlan{IDCTPLAN,T}
+# lex indicates if its lexigraphical (i.e., x, y) or reverse (y, x)
+immutable IPaduaTransformPlan{lex,IDCTPLAN,T}
     cfsmat::Matrix{T}
     padvals::Vector{T}
     idctplan::IDCTPLAN
 end
 
-function plan_ipaduatransform{T}(::Type{T},N::Integer)
+IPaduaTransformPlan{T,lex}(cfsmat::Matrix{T},padvals::Vector{T},idctplan,::Type{Val{lex}}) =
+    IPaduaTransformPlan{lex,typeof(idctplan),T}(cfsmat,padvals,idctplan)
+
+function plan_ipaduatransform{T}(::Type{T},N::Integer,lex)
     n=Int(cld(-3+sqrt(1+8N),2))
     if N ≠ div((n+1)*(n+2),2)
         error("Padua transforms can only be applied to vectors of length (n+1)*(n+2)/2.")
     end
     IPaduaTransformPlan(Array{T}(n+2,n+1),Array{T}(N),
-        FFTW.plan_r2r!(Array{T}(n+2,n+1),FFTW.REDFT00))
+        FFTW.plan_r2r!(Array{T}(n+2,n+1),FFTW.REDFT00),lex)
 end
 
-plan_ipaduatransform{T}(v::AbstractVector{T}) = plan_ipaduatransform(eltype(v),length(v))
+
+plan_ipaduatransform{T}(::Type{T},N::Integer) = plan_ipaduatransform(T,N,Val{true})
+
+plan_ipaduatransform{T}(v::AbstractVector{T},lex...) = plan_ipaduatransform(eltype(v),length(v),lex...)
 
 """
 Inverse Padua Transform maps the 2D Chebyshev coefficients to the values of the interpolation polynomial at the Padua points.
@@ -31,12 +38,12 @@ function *{T}(P::IPaduaTransformPlan,v::AbstractVector{T})
     return paduavals
 end
 
-ipaduatransform(v::AbstractVector) = plan_ipaduatransform(v)*v
+ipaduatransform(v::AbstractVector,lex...) = plan_ipaduatransform(v,lex...)*v
 
 """
 Creates (n+2)x(n+1) Chebyshev coefficient matrix from triangle coefficients.
 """
-function trianglecfsmat(P::IPaduaTransformPlan,cfs::AbstractVector)
+function trianglecfsmat(P::IPaduaTransformPlan{true},cfs::AbstractVector)
     N=length(cfs)
     n=Int(cld(-3+sqrt(1+8N),2))
     cfsmat=fill!(P.cfsmat,0)
@@ -55,6 +62,25 @@ function trianglecfsmat(P::IPaduaTransformPlan,cfs::AbstractVector)
     return cfsmat
 end
 
+function trianglecfsmat(P::IPaduaTransformPlan{false},cfs::AbstractVector)
+    N=length(cfs)
+    n=Int(cld(-3+sqrt(1+8N),2))
+    cfsmat=fill!(P.cfsmat,0)
+    m=1
+    for d=1:n+1
+        @inbounds for k=d:-1:1
+            j=d-k+1
+            cfsmat[k,j]=cfs[m]
+            if m==N
+                return cfsmat
+            else
+                m+=1
+            end
+        end
+    end
+    return cfsmat
+end
+
 """
 Vectorizes the function values at the Padua points.
 """
@@ -77,22 +103,26 @@ end
 """
 Pre-plan a Padua Transform.
 """
-immutable PaduaTransformPlan{DCTPLAN,T}
+immutable PaduaTransformPlan{lex,DCTPLAN,T}
     vals::Matrix{T}
     retvec::Vector{T}
     dctplan::DCTPLAN
 end
 
-function plan_paduatransform{T}(::Type{T},N::Integer)
+PaduaTransformPlan{T,lex}(vals::Matrix{T},retvec::Vector{T},dctplan,::Type{Val{lex}}) =
+    PaduaTransformPlan{lex,typeof(dctplan),T}(vals,retvec,dctplan)
+
+function plan_paduatransform{T}(::Type{T},N::Integer,lex)
     n=Int(cld(-3+sqrt(1+8N),2))
     if N ≠ ((n+1)*(n+2))÷2
         error("Padua transforms can only be applied to vectors of length (n+1)*(n+2)/2.")
     end
     PaduaTransformPlan(Array{T}(n+2,n+1),Array{T}(N),
-        FFTW.plan_r2r!(Array{T}(n+2,n+1),FFTW.REDFT00))
+        FFTW.plan_r2r!(Array{T}(n+2,n+1),FFTW.REDFT00),lex)
 end
 
-plan_paduatransform{T}(v::AbstractVector{T}) = plan_paduatransform(eltype(v),length(v))
+plan_paduatransform{T}(::Type{T},N::Integer) = plan_paduatransform(T,N,Val{true})
+plan_paduatransform{T}(v::AbstractVector{T},lex...) = plan_paduatransform(eltype(v),length(v),lex...)
 
 """
 Padua Transform maps from interpolant values at the Padua points to the 2D Chebyshev coefficients.
@@ -112,7 +142,7 @@ function *{T}(P::PaduaTransformPlan,v::AbstractVector{T})
     return cfs
 end
 
-paduatransform(v::AbstractVector) = plan_paduatransform(v)*v
+paduatransform(v::AbstractVector,lex...) = plan_paduatransform(v,lex...)*v
 
 """
 Creates (n+2)x(n+1) matrix of interpolant values on the tensor grid at the (n+1)(n+2)/2 Padua points.
@@ -137,7 +167,7 @@ end
 """
 Creates length (n+1)(n+2)/2 vector from matrix of triangle Chebyshev coefficients.
 """
-function trianglecfsvec(P::PaduaTransformPlan,cfs::Matrix)
+function trianglecfsvec(P::PaduaTransformPlan{true},cfs::Matrix)
     m=size(cfs,2)
     l=1
     for d=1:m
@@ -150,6 +180,19 @@ function trianglecfsvec(P::PaduaTransformPlan,cfs::Matrix)
     return P.retvec
 end
 
+function trianglecfsvec(P::PaduaTransformPlan{false},cfs::Matrix)
+    m=size(cfs,2)
+    l=1
+    for d=1:m
+        @inbounds for k=d:-1:1
+            j=d-k+1
+            P.retvec[l]=cfs[k,j]
+            l+=1
+        end
+    end
+    return P.retvec
+end
+
 """
 Returns coordinates of the (n+1)(n+2)/2 Padua points.
 """

test/runtests.jl

@@ -223,37 +223,19 @@ println("Testing runtimes for (I)Padua Transforms")
 @time IPl*v
 
 println("Accuracy of 2d function interpolation at a point")
-function trianglecfsmat{T}(cfs::AbstractVector{T})
-    N=length(cfs)
-    n=Int(cld(-3+sqrt(1+8N),2))
-    @assert N==div((n+1)*(n+2),2)
-    cfsmat=Array(T,n+2,n+1)
-    cfsmat=fill!(cfsmat,0)
-    m=1
-    for d=1:n+1
-        @inbounds for k=1:d
-            j=d-k+1
-            cfsmat[k,j]=cfs[m]
-            if m==N
-                return cfsmat
-            else
-                m+=1
-            end
-        end
-    end
-    return cfsmat
-end
+
 """
 Interpolates a 2d function at a given point using 2d Chebyshev series.
 """
-function paduaeval(f::Function,x::AbstractFloat,y::AbstractFloat,m::Integer)
+function paduaeval(f::Function,x::AbstractFloat,y::AbstractFloat,m::Integer,lex)
     T=promote_type(typeof(x),typeof(y))
     M=div((m+1)*(m+2),2)
     pvals=Array(T,M)
     p=paduapoints(T,m)
     map!(f,pvals,p[:,1],p[:,2])
-    coeffs=paduatransform(pvals)
-    cfs_mat=trianglecfsmat(coeffs)
+    coeffs=paduatransform(pvals,lex)
+    plan=plan_ipaduatransform(pvals,lex)
+    cfs_mat=FastTransforms.trianglecfsmat(plan,coeffs)
     f_x=sum([cfs_mat[k,j]*cos((j-1)*acos(x))*cos((k-1)*acos(y)) for k=1:m+1, j=1:m+1])
     return f_x
 end
@@ -262,7 +244,13 @@ g_xy = (x,y) ->cos(exp(2*x+y))*sin(y)
 x=0.1;y=0.2
 m=130
 l=80
-f_m=paduaeval(f_xy,x,y,m)
-g_l=paduaeval(g_xy,x,y,l)
+f_m=paduaeval(f_xy,x,y,m,Val{true})
+g_l=paduaeval(g_xy,x,y,l,Val{true})
+@test_approx_eq f_xy(x,y) f_m
+@test_approx_eq g_xy(x,y) g_l
+
+
+f_m=paduaeval(f_xy,x,y,m,Val{false})
+g_l=paduaeval(g_xy,x,y,l,Val{false})
 @test_approx_eq f_xy(x,y) f_m
 @test_approx_eq g_xy(x,y) g_l