Add plan_transform (#106)

* Add plan_transform

* Work on plan_transform

* Update ClassicalOrthogonalPolynomials.jl

* Move cardinality

* Add MulPlan

* Update interlace.jl

* plan_grid_transform

* Update ci.yml

* add transform tests

* Update test_normalized.jl

Author: dlfivefifty

.github/workflows/ci.yml

@@ -11,7 +11,7 @@ jobs:
       matrix:
         version:
           - '1.7'
-          - '^1.8.0-0'
+          - '1'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.6.5"
+version = "0.6.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,11 +29,11 @@ ArrayLayouts = "0.8"
 BandedMatrices = "0.17"
 BlockArrays = "0.16.9"
 BlockBandedMatrices = "0.11.6"
-ContinuumArrays = "0.10"
+ContinuumArrays = "0.11"
 DomainSets = "0.5.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
-FastTransforms = "0.13, 0.14"
+FastTransforms = "0.14.4"
 FillArrays = "0.13"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.3"

src/ClassicalOrthogonalPolynomials.jl

@@ -36,7 +36,8 @@ import InfiniteLinearAlgebra: chop!, chop, choplength, compatible_resize!
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
     inbounds_getindex, grid, plotgrid, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
     AffineQuasiVector, AffineMap, WeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
-    checkpoints, weight, unweighted, MappedBasisLayouts, __sum, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, _broadcastbasis
+    checkpoints, weight, unweighted, MappedBasisLayouts, __sum, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, _broadcastbasis,
+    plan_transform, plan_grid_transform
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints, Plan
 
@@ -52,20 +53,13 @@ export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomia
             ∞, Derivative, .., Inclusion,
             chebyshevt, chebyshevu, legendre, jacobi, ultraspherical,
             legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh, normalizedjacobip,
-            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted, PiecewiseInterlace
+            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted, PiecewiseInterlace, plan_transform
 
 
 import Base: oneto
 
 
 include("interlace.jl")
-
-
-cardinality(::FullSpace{<:AbstractFloat}) = ℵ₁
-cardinality(::EuclideanDomain) = ℵ₁
-cardinality(::Union{DomainSets.RealNumbers,DomainSets.ComplexNumbers}) = ℵ₁
-cardinality(::Union{DomainSets.Integers,DomainSets.Rationals,DomainSets.NaturalNumbers}) = ℵ₀
-
 include("standardchop.jl")
 include("adaptivetransform.jl")
 
@@ -95,9 +89,9 @@ _equals(::MappedOPLayout, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) && basi
 _equals(::MappedOPLayout, ::MappedBasisLayouts, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
 _equals(::MappedBasisLayouts, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
 
-_broadcastbasis(::typeof(+), ::MappedOPLayout, ::MappedOPLayout, P, Q) where {L,M} = _broadcastbasis(+, MappedBasisLayout(), MappedBasisLayout(), P, Q)
-_broadcastbasis(::typeof(+), ::MappedOPLayout, M::MappedBasisLayout, P, Q) where L = _broadcastbasis(+, MappedBasisLayout(), M, P, Q)
-_broadcastbasis(::typeof(+), L::MappedBasisLayout, ::MappedOPLayout, P, Q) where M = _broadcastbasis(+, L, MappedBasisLayout(), P, Q)
+_broadcastbasis(::typeof(+), ::MappedOPLayout, ::MappedOPLayout, P, Q) = _broadcastbasis(+, MappedBasisLayout(), MappedBasisLayout(), P, Q)
+_broadcastbasis(::typeof(+), ::MappedOPLayout, M::MappedBasisLayout, P, Q) = _broadcastbasis(+, MappedBasisLayout(), M, P, Q)
+_broadcastbasis(::typeof(+), L::MappedBasisLayout, ::MappedOPLayout, P, Q) = _broadcastbasis(+, L, MappedBasisLayout(), P, Q)
 __sum(::MappedOPLayout, A, dims) = __sum(MappedBasisLayout(), A, dims)
 
 # demap to avoid Golub-Welsch fallback
@@ -231,11 +225,6 @@ function recurrencecoefficients(C::SubQuasiArray{T,2,<:Any,<:Tuple{AbstractAffin
     A * kr.A, A*kr.b + B, C
 end
 
-
-_vec(a) = vec(a)
-_vec(a::InfiniteArrays.ReshapedArray) = _vec(parent(a))
-_vec(a::Adjoint{<:Any,<:AbstractVector}) = a'
-
 include("clenshaw.jl")
 include("ratios.jl")
 include("normalized.jl")
@@ -270,24 +259,76 @@ function golubwelsch(V::SubQuasiArray)
     x,w
 end
 
-function factorize(L::SubQuasiArray{T,2,<:Normalized,<:Tuple{Inclusion,OneTo}}, dims...; kws...) where T
-    x,w = golubwelsch(L)
-    TransformFactorization(x, L[x,:]'*Diagonal(w))
+"""
+    MulPlan(matrix, dims)
+
+Takes a matrix and supports it applied to different dimensions.
+"""
+struct MulPlan{T, Fact, Dims} # <: Plan{T} We don't depend on AbstractFFTs
+    matrix::Fact
+    dims::Dims
 end
 
+MulPlan(fact, dims) = MulPlan{eltype(fact), typeof(fact), typeof(dims)}(fact, dims)
+
+function *(P::MulPlan{<:Any,<:Any,Int}, x::AbstractVector)
+    @assert P.dims == 1
+    P.matrix * x
+end
 
-function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,OneTo}}, dims...; kws...) where T
-    Q = Normalized(parent(L))[parentindices(L)...]
-    D = L \ Q
-    F = factorize(Q, dims...; kws...)
-    TransformFactorization(F.grid, D*F.plan)
+function *(P::MulPlan{<:Any,<:Any,Int}, X::AbstractMatrix)
+    if P.dims == 1
+        P.matrix * X
+    else
+        @assert P.dims == 2
+        permutedims(P.matrix * permutedims(X))
+    end
 end
 
-function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{<:Inclusion,<:AbstractUnitRange}}, dims...; kws...) where T
-    _,jr = parentindices(L)
-    ProjectionFactorization(factorize(parent(L)[:,oneto(maximum(jr))], dims...; kws...), jr)
+function *(P::MulPlan{<:Any,<:Any,Int}, X::AbstractArray{<:Any,3})
+    Y = similar(X)
+    if P.dims == 1
+        for j in axes(X,3)
+            Y[:,:,j] = P.matrix * X[:,:,j]
+        end
+    elseif P.dims == 2
+        for k in axes(X,1)
+            Y[k,:,:] = P.matrix * X[k,:,:]
+        end
+    else
+        @assert P.dims == 3
+        for k in axes(X,1), j in axes(X,2)
+            Y[k,j,:] = P.matrix * X[k,j,:]
+        end
+    end
+    Y
 end
 
+function *(P::MulPlan, X::AbstractArray)
+    for d in P.dims
+        X = MulPlan(P.matrix, d) * X
+    end
+    X
+end
+
+*(A::AbstractMatrix, P::MulPlan) = MulPlan(A*P.matrix, P.dims)
+
+
+function plan_grid_transform(Q::Normalized, arr, dims=1:ndims(arr))
+    L = Q[:,OneTo(size(arr,1))]
+    x,w = golubwelsch(L)
+    x, MulPlan(L[x,:]'*Diagonal(w), dims)
+end
+
+function plan_grid_transform(P::OrthogonalPolynomial, arr, dims...)
+    Q = Normalized(P)
+    x, A = plan_grid_transform(Q, arr, dims...)
+    n = size(arr,1)
+    D = (P \ Q)[1:n, 1:n]
+    x, D * A
+end
+
+
 function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial})
     axes(A,1) == axes(B,1) || throw(DimensionMismatch())
     _,jA = parentindices(A)

src/classical/chebyshev.jl

@@ -103,15 +103,16 @@ Jacobi(C::ChebyshevU{T}) where T = Jacobi(one(T)/2,one(T)/2)
 #######
 
 
-factorize(L::SubQuasiArray{T,2,<:ChebyshevT,<:Tuple{<:Inclusion,<:OneTo}}) where T =
-    TransformFactorization(grid(L), plan_chebyshevtransform(Array{T}(undef, size(L,2))))
-
-factorize(L::SubQuasiArray{T,2,<:ChebyshevT,<:Tuple{<:Inclusion,<:OneTo}}, m) where T =
-    TransformFactorization(grid(L), plan_chebyshevtransform(Array{T}(undef, size(L,2), m),1))
-
-# TODO: extend plan_chebyshevutransform
-factorize(L::SubQuasiArray{T,2,<:ChebyshevU,<:Tuple{<:Inclusion,<:OneTo}}) where T<:FastTransforms.fftwNumber =
-    TransformFactorization(grid(L), plan_chebyshevutransform(Array{T}(undef, size(L,2))))
+function plan_grid_transform(T::ChebyshevT, arr, dims...)
+    n = size(arr,1)
+    x = grid(T[:,oneto(n)])
+    x, plan_chebyshevtransform(arr, dims...)
+end
+function plan_grid_transform(U::ChebyshevU{<:FastTransforms.fftwNumber}, arr, dims...)
+    n = size(arr,1)
+    x = grid(U[:,oneto(n)])
+    x, plan_chebyshevutransform(arr, dims...)
+end
 
 
 ########

src/interlace.jl

@@ -241,7 +241,7 @@ function \(F::SetindexFactorization{T,<:AbstractFill}, v::AbstractQuasiVector) w
     for (k,x) in enumerate(F̃.grid)
         data[k,:] = v[x]
     end
-    blockvec(Matrix(transpose(F̃ \ data))) # call Matrix to avoid ReshapedArray
+    blockvec(Matrix(transpose(F̃.plan * data))) # call Matrix to avoid ReshapedArray
 end
 
 

test/test_chebyshev.jl

@@ -1,4 +1,4 @@
-using ClassicalOrthogonalPolynomials, QuasiArrays, ContinuumArrays, BandedMatrices, LazyArrays, 
+using ClassicalOrthogonalPolynomials, QuasiArrays, ContinuumArrays, BandedMatrices, LazyArrays,
         FastTransforms, ArrayLayouts, Test, FillArrays, Base64, BlockArrays, LazyBandedMatrices, ForwardDiff
 import ClassicalOrthogonalPolynomials: Clenshaw, recurrencecoefficients, clenshaw, paddeddata, jacobimatrix, oneto, Weighted, MappedOPLayout
 import LazyArrays: ApplyStyle
@@ -120,6 +120,13 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
             @time U = T / T \ cos.(x .* (0:1000)');
             @test U[0.1,:] ≈ cos.(0.1 * (0:1000))
             @test U[[0.1,0.2],:] ≈ [cos.(0.1 * (0:1000)'); cos.(0.2 * (0:1000)')]
+
+            # support tensors but for grids
+            X = randn(150, 2, 2)
+            F = plan_transform(T, X, 1)
+            F_m = plan_transform(T, X[:,:,1], 1)
+            @test (F * X)[:,:,1] ≈ F_m * X[:,:,1]
+            @test (F * X)[:,:,2] ≈ F_m * X[:,:,2]
         end
     end
 
@@ -183,11 +190,11 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
             @test sum(wT; dims=1)[:,1:10] ≈ [π zeros(1,9)]
             @test sum(wT[:,1]) ≈ π
             @test iszero(sum(wT[:,2]))
-  
+
             @test (wT \ wU)[1:10,1:10] ≈ inv(wU \ wT)[1:10,1:10]
             @test_skip (wU \ WT)[1,1] == 2
         end
-        
+
         @testset "Derivative" begin
             x = axes(wT,1)
             D = Derivative(x)
@@ -445,7 +452,7 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
         T = ChebyshevT(); U = ChebyshevU()
         D = Derivative(axes(T,1))
         V = view(T,:,[1,3,4])
-        @test (U\(D*V))[1:5,:] == (U \ (V'D')')[1:5,:] == (U\(D*T))[1:5,[1,3,4]] 
+        @test (U\(D*V))[1:5,:] == (U \ (V'D')')[1:5,:] == (U\(D*T))[1:5,[1,3,4]]
     end
 
     @testset "plot" begin

test/test_normalized.jl

@@ -209,4 +209,47 @@ import ContinuumArrays: MappedWeightedBasisLayout
         @test grid(Q[:,1:5]) == grid(Q[:,collect(1:5)]) == grid(P[:,1:5])
         @test plotgrid(Q[:,1:5]) == plotgrid(Q[:,collect(1:5)]) == plotgrid(P[:,1:5])
     end
+
+    @testset "Transform" begin
+        Q = Normalized(Hermite())
+        n = 20
+        Qₙ = Q[:,Base.OneTo(n)]
+        x = axes(Q,1)
+        g = grid(Qₙ)
+        v = exp.(g)
+        P = plan_transform(Q, v)
+        @test P * v ≈ Qₙ[g,:] \ exp.(g) ≈ transform(Qₙ, exp)
+
+        V = cos.(g .* (1:3)')
+        P = plan_transform(Q, V, 1)
+        @test P * V ≈ Qₙ \ cos.(x .* (1:3)')
+
+        X = randn(n, n)
+        P₂ = plan_transform(Q, X, 2)
+
+        P = plan_transform(Q, X)
+
+        PX = P * X
+        for k = 1:n
+            X[:, k] = Qₙ[g,:] \ X[:, k]
+        end
+        for k = 1:n
+            X[k, :] = Qₙ[g,:] \ X[k, :]
+        end
+        @test PX ≈ X
+
+        X = randn(n, n, n)
+        P = plan_transform(Q, X)
+        PX = P * X
+        for k = 1:n, j = 1:n
+            X[:, k, j] = Qₙ[g,:] \ X[:, k, j]
+        end
+        for k = 1:n, j = 1:n
+            X[k, :, j] = Qₙ[g,:] \ X[k, :, j]
+        end
+        for k = 1:n, j = 1:n
+            X[k, j, :] = Qₙ[g,:] \ X[k, j, :]
+        end
+        @test PX ≈ X
+    end
 end