Introduce plan_transform (#134)

* Introduce plan_transform

* v0.11

* Update bases.jl

* Update bases.jl

* Add transform and expand

* Derivative -> Diff

* Don't call using LazyArrays

* Update test_splines.jl

* Update bases.jl

* fix FactorizationPlan

* Update Project.toml

* Diff -> Derivative

* Diff -> Derivative

* Update ci.yml

* Update Project.toml

* Update Project.toml

* Revert "Update Project.toml"

This reverts commit c1fe83e.

* Revert "Update Project.toml"

This reverts commit f4ebd28.

* increase coverage

* Add plan_transform tests

* Increase coverage

* Update test_chebyshev.jl

* split plan_grid_transform and plan_transform

Author: dlfivefifty

.github/workflows/ci.yml

@@ -12,7 +12,6 @@ jobs:
         version:
           - '1.6'
           - '1'
-          - '^1.8.0-0'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,11 +1,12 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.10.3"
+version = "0.11"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
@@ -20,13 +21,14 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 ArrayLayouts = "0.7.7, 0.8"
 BandedMatrices = "0.16, 0.17"
 BlockArrays = "0.16"
-FastTransforms = "0.13"
+DomainSets = "0.5"
+FastTransforms = "0.13, 0.14"
 FillArrays = "0.12, 0.13"
 InfiniteArrays = "0.12"
 Infinities = "0.1"
 IntervalSets = "0.5, 0.6, 0.7"
 LazyArrays = "0.22"
-QuasiArrays = "0.9"
+QuasiArrays = "0.9.3"
 RecipesBase = "1.0"
 StaticArrays = "1.0"
 julia = "1.6"

README.md

@@ -69,11 +69,13 @@ julia> f[0.1]
 
 Creating a finite element method is possible using standard array terminology. 
 We always take the Lebesgue inner product associated with an axes, so in this
-case the mass matrix is just `L'L`. Combined with a derivative operator allows
+case the mass matrix is just `L'L`. Combined with a differentiation operator allows
 us to form the weak Laplacian.
 ```julia
 julia> B = L[:,2:end-1]; # drop boundary terms to impose zero Dirichlet
 
+julia> D = Derivative(L); # Differentiation operator
+
 julia> Δ = (D*B)'D*B # weak Laplacian
 4×4 BandedMatrices.BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
  10.0  -5.0    ⋅     ⋅ 

examples/heat.jl

@@ -9,9 +9,9 @@ x = axes(L,1)
 D = Derivative(x)
 M = L'L
 Δ = -((D*L)'D*L)
-u0 = copy(L \ exp.(x))
+u0 = L \ exp.(x)
 
-heat(u,(M,Δ),t) = M\(Δ*u)
+heat(u,(M,Δ),t) = M \ (Δ*u)
 prob = ODEProblem(heat,u0,(0.0,1.0),(M,Δ))
 sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
 

examples/poisson.jl

@@ -5,11 +5,10 @@ using ContinuumArrays, Plots
 ####
 
 L = LinearSpline(range(0,1; length=10_000))[:,2:end-1]
-x = axes(L,1)
-D = Derivative(x)
+D = Derivative(L)
 Δ = -((D*L)'D*L)
 M = L'L
-f = L \ exp.(x)
+f = expand(L, exp)
 u = L * (Δ \ (M*f))
 plot(u)
 

src/ContinuumArrays.jl

@@ -1,5 +1,5 @@
 module ContinuumArrays
-using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays, Infinities, InfiniteArrays, StaticArrays, BlockArrays, RecipesBase
+using IntervalSets, DomainSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays, Infinities, InfiniteArrays, StaticArrays, BlockArrays, RecipesBase
 import Base: @_inline_meta, @_propagate_inbounds_meta, axes, size, getindex, convert, prod, *, /, \, +, -, ==, ^,
                 IndexStyle, IndexLinear, ==, OneTo, _maybetail, tail, similar, copyto!, copy, diff,
                 first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum, sum, _sum,
@@ -21,14 +21,16 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     AbstractQuasiFill, UnionDomain, __sum, _cumsum, __cumsum, applylayout, _equals, layout_broadcasted, PolynomialLayout
 import InfiniteArrays: Infinity, InfAxes
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, ..
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform
 
 
 
 const QMul2{A,B} = Mul{<:Any, <:Any, <:A,<:B}
 const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}
 
 cardinality(::AbstractInterval) = ℵ₁
+cardinality(::Union{FullSpace{<:AbstractFloat},EuclideanDomain,DomainSets.RealNumbers,DomainSets.ComplexNumbers}) = ℵ₁
+cardinality(::Union{DomainSets.Integers,DomainSets.Rationals,DomainSets.NaturalNumbers}) = ℵ₀
 
 Inclusion(d::AbstractInterval{T}) where T = Inclusion{float(T)}(d)
 first(S::Inclusion{<:Any,<:AbstractInterval}) = leftendpoint(S.domain)
@@ -90,7 +92,6 @@ checkpoints(A::AbstractQuasiMatrix) = checkpoints(axes(A,1))
 
 include("operators.jl")
 include("bases/bases.jl")
-include("basisconcat.jl")
 
 include("plotting.jl")
 

src/bases/bases.jl

@@ -156,13 +156,12 @@ _grid(::WeightedBasisLayouts, P) = grid(unweighted(P))
 grid(P) = _grid(MemoryLayout(P), P)
 
 
-struct TransformFactorization{T,Grid,Plan,IPlan} <: Factorization{T}
+struct TransformFactorization{T,Grid,Plan} <: Factorization{T}
     grid::Grid
     plan::Plan
-    iplan::IPlan
 end
 
-TransformFactorization{T}(grid, plan) where T = TransformFactorization{T,typeof(grid),typeof(plan),Nothing}(grid, plan, nothing)
+TransformFactorization{T}(grid, plan) where T = TransformFactorization{T,typeof(grid),typeof(plan)}(grid, plan)
 
 """
     TransformFactorization(grid, plan)
@@ -173,44 +172,84 @@ associates a planned transform with a grid. That is, if `F` is a `TransformFacto
 TransformFactorization(grid, plan) = TransformFactorization{promote_type(eltype(eltype(grid)),eltype(plan))}(grid, plan)
 
 
-TransformFactorization{T}(grid, ::Nothing, iplan) where T = TransformFactorization{T,typeof(grid),Nothing,typeof(iplan)}(grid, nothing, iplan)
-
-"""
-    TransformFactorization(grid, nothing, iplan)
-
-associates a planned inverse transform with a grid. That is, if `F` is a `TransformFactorization`, then
-`F \\ f` is equivalent to `F.iplan \\ f[F.grid]`.
-"""
-TransformFactorization(grid, ::Nothing, iplan) = TransformFactorization{promote_type(eltype(eltype(grid)),eltype(iplan))}(grid, nothing, iplan)
 
 grid(T::TransformFactorization) = T.grid
 function size(T::TransformFactorization, k)
     @assert k == 2 # TODO: make consistent
     size(T.plan,1)
 end
 
-function size(T::TransformFactorization{<:Any,<:Any,Nothing}, k)
-    @assert k == 2 # TODO: make consistent
-    size(T.iplan,2)
+
+\(a::TransformFactorization, b::AbstractQuasiVector) = a.plan * convert(Array, b[a.grid])
+\(a::TransformFactorization, b::AbstractQuasiMatrix) = a.plan * convert(Array, b[a.grid,:])
+
+"""
+    InvPlan(factorization, dims)
+
+Takes a factorization and supports it applied to different dimensions.
+"""
+struct InvPlan{T, Fact, Dims} # <: Plan{T} We don't depend on AbstractFFTs
+    factorization::Fact
+    dims::Dims
 end
 
+InvPlan(fact, dims) = InvPlan{eltype(fact), typeof(fact), typeof(dims)}(fact, dims)
 
-\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiVector{T}) where T = a.iplan \  convert(Array{T}, b[a.grid])
-\(a::TransformFactorization, b::AbstractQuasiVector) = a.plan * convert(Array, b[a.grid])
-\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractVector) = a.iplan \  b
-\(a::TransformFactorization, b::AbstractVector) = a.plan * b
-ldiv!(ret::AbstractVecOrMat, a::TransformFactorization, b::AbstractVecOrMat) = mul!(ret, a.plan, b)
-ldiv!(ret::AbstractVecOrMat, a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractVecOrMat) = ldiv!(ret, a.iplan, b)
+size(F::InvPlan, k...) = size(F.factorization, k...)
+
+
+function *(P::InvPlan{<:Any,<:Any,Int}, x::AbstractVector)
+    @assert P.dims == 1
+    P.factorization \ x
+end
+
+function *(P::InvPlan{<:Any,<:Any,Int}, X::AbstractMatrix)
+    if P.dims == 1
+        P.factorization \ X
+    else
+        @assert P.dims == 2
+        permutedims(P.factorization \ permutedims(X))
+    end
+end
 
-\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiMatrix{T}) where T = a \  convert(Array{T}, b[a.grid,:])
-\(a::TransformFactorization, b::AbstractQuasiMatrix) = a \ convert(Array, b[a.grid,:])
-\(a::TransformFactorization, b::AbstractMatrix) = ldiv!(Array{promote_type(eltype(a),eltype(b))}(undef,size(a,2),size(b,2)), a, b)
+function *(P::InvPlan{<:Any,<:Any,Int}, X::AbstractArray{<:Any,3})
+    Y = similar(X)
+    if P.dims == 1
+        for j in axes(X,3)
+            Y[:,:,j] = P.factorization \ X[:,:,j]
+        end
+    elseif P.dims == 2
+        for k in axes(X,1)
+            Y[k,:,:] = P.factorization \ X[k,:,:]
+        end
+    else
+        @assert P.dims == 3
+        for k in axes(X,1), j in axes(X,2)
+            Y[k,j,:] = P.factorization \ X[k,j,:]
+        end
+    end
+    Y
+end
 
-function _factorize(::AbstractBasisLayout, L, dims...; kws...)
+function *(P::InvPlan, X::AbstractArray)
+    for d in P.dims
+        X = InvPlan(P.factorization, d) * X
+    end
+    X
+end
+
+
+function plan_grid_transform(L, arr, dims=1:ndims(arr))
     p = grid(L)
-    TransformFactorization(p, nothing, factorize(L[p,:]))
+    p, InvPlan(factorize(L[p,:]), dims)
 end
 
+plan_transform(P, arr, dims...) = plan_grid_transform(P, arr, dims...)[2]
+
+_factorize(::AbstractBasisLayout, L, dims...; kws...) =
+    TransformFactorization(plan_grid_transform(L, Array{eltype(L)}(undef, size(L,2), dims...), 1)...)
+
+
 
 """
     ProjectionFactorization(F, inds)
@@ -225,9 +264,22 @@ end
 
 \(a::ProjectionFactorization, b::AbstractQuasiVector) = (a.F \ b)[a.inds]
 \(a::ProjectionFactorization, b::AbstractQuasiMatrix) = (a.F \ b)[a.inds,:]
-\(a::ProjectionFactorization, b::AbstractVector) = (a.F \ b)[a.inds]
 
-_factorize(::SubBasisLayout, L, dims...; kws...) = ProjectionFactorization(factorize(parent(L), dims...; kws...), parentindices(L)[2])
+
+
+# if parent is finite dimensional default to its transform and project down
+_sub_factorize(::Tuple{Any,Int}, (kr,jr), L, dims...; kws...) = ProjectionFactorization(factorize(parent(L), dims...; kws...), jr)
+_sub_factorize(::Tuple{Any,Int}, (kr,jr)::Tuple{Any,OneTo}, L, dims...; kws...) = ProjectionFactorization(factorize(parent(L), dims...; kws...), jr)
+
+# ∞-dimensional parents need to use transforms. For now we assume the size of the transform is equal to the size of the truncation
+_sub_factorize(::Tuple{Any,Any}, (kr,jr)::Tuple{Any,OneTo}, L, dims...; kws...) =
+    TransformFactorization(plan_grid_transform(parent(L), Array{eltype(L)}(undef, last(jr), dims...), 1)...)
+
+# If jr is not OneTo we project
+_sub_factorize(::Tuple{Any,Any}, (kr,jr), L, dims...; kws...) =
+    ProjectionFactorization(factorize(parent(L)[:,OneTo(maximum(jr))]), jr)
+
+_factorize(::SubBasisLayout, L, dims...; kws...) = _sub_factorize(size(parent(L)), parentindices(L), L, dims...; kws...)
 
 
 """
@@ -260,6 +312,29 @@ plan_ldiv(A, B::AbstractQuasiMatrix) = factorize(A, size(B,2))
 transform_ldiv(A::AbstractQuasiArray{T}, B::AbstractQuasiArray{V}, _) where {T,V} = plan_ldiv(A, B) \ B
 transform_ldiv(A, B) = transform_ldiv(A, B, size(A))
 
+
+"""
+    transform(A, f)
+
+finds the coefficients of a function `f` expanded in a basis defined as the columns of a quasi matrix `A`.
+It is equivalent to
+```
+A \\ f.(axes(A,1))
+```
+"""
+transform(A, f) = A \ f.(axes(A,1))
+
+"""
+    expand(A, f)
+
+expands a function `f` im a basis defined as the columns of a quasi matrix `A`.
+It is equivalent to
+```
+A / A \\ f.(axes(A,1))
+```
+"""
+expand(A, f) = A * transform(A, f)
+
 copy(L::Ldiv{<:AbstractBasisLayout}) = transform_ldiv(L.A, L.B)
 # TODO: redesign to use simplifiable(\, A, B)
 copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = transform_ldiv(L.A, L.B)
@@ -573,4 +648,6 @@ end
 __sum(::ExpansionLayout, A, dims) = __sum(ApplyLayout{typeof(*)}(), A, dims)
 __cumsum(::ExpansionLayout, A, dims) = __cumsum(ApplyLayout{typeof(*)}(), A, dims)
 
+include("basisconcat.jl")
+include("basiskron.jl")
 include("splines.jl")

src/basisconcat.jl renamed to src/bases/basisconcat.jl

@@ -0,0 +1,3 @@
+struct KronBasisLayout <: AbstractBasisLayout end
+
+QuasiArrays.kronlayout(::AbstractBasisLayout...) = KronBasisLayout()

src/bases/basiskron.jl

@@ -94,7 +94,7 @@ end
 end
 
 
-## Derivative
+## Differentiation
 function copyto!(dest::MulQuasiMatrix{<:Any,<:Tuple{<:HeavisideSpline,<:Any}},
                  M::QMul2{<:Derivative,<:LinearSpline})
     D, L = M.A, M.B

src/bases/splines.jl

@@ -92,17 +92,24 @@ show(io::IO, δ::DiracDelta) = print(io, "δ at $(δ.x) over $(axes(δ,1))")
 show(io::IO, ::MIME"text/plain", δ::DiracDelta) = show(io, δ)
 
 #########
-# Derivative
+# Differentiation
 #########
 
+"""
+Derivative(axis)
 
+represents the differentiation (or finite-differences) operator on the
+specified axis.
+"""
 struct Derivative{T,D} <: LazyQuasiMatrix{T}
     axis::Inclusion{T,D}
 end
 
 Derivative{T}(axis::Inclusion{<:Any,D}) where {T,D} = Derivative{T,D}(axis)
 Derivative{T}(domain) where T = Derivative{T}(Inclusion(domain))
 
+Derivative(L::AbstractQuasiMatrix) = Derivative(axes(L,1))
+
 function summary(io::IO, D::Derivative)
     print(io, "Derivative(")
     summary(io,D.axis)

src/operators.jl

@@ -1,4 +1,4 @@
-using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, InfiniteArrays, Test, Base64, RecipesBase
+using ContinuumArrays, QuasiArrays, IntervalSets, DomainSets, FillArrays, LinearAlgebra, BandedMatrices, InfiniteArrays, Test, Base64, RecipesBase
 import ContinuumArrays: ℵ₁, materialize, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, ℵ₁,
                         MappedBasisLayout, AdjointMappedBasisLayout, MappedWeightedBasisLayout, TransformFactorization, Weight, WeightedBasisLayout, SubWeightedBasisLayout, WeightLayout,
                         basis, invmap, Map, checkpoints, _plotgrid, mul
@@ -12,6 +12,9 @@ import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, App
         @test eltype(x) == Float64
         @test x[0.1] ≡ 0.1
         @test x[0] ≡ x[0.0] ≡ 0.0
+
+        @test size(Inclusion(ℝ),1) == ℵ₁
+        @test size(Inclusion(ℤ),1) == ℵ₀
     end
 
     @testset "broadcast" begin

test/runtests.jl

@@ -0,0 +1,11 @@
+using ContinuumArrays, QuasiArrays, StaticArrays, Test
+
+@testset "Basis Kron" begin
+    L = LinearSpline(range(0,1; length=4))
+    K = QuasiKron(L, L)
+
+    xy = axes(K,1)
+    f = xy -> ((x,y) = xy; exp(x*cos(y)))
+    K \ f.(xy)
+
+end

test/test_basiskron.jl

@@ -52,7 +52,7 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
     @testset "basics" begin
         F = factorize(T)
         g = grid(F)
-        @test T \ exp.(x) == F \ exp.(x) == F \ exp.(g) == chebyshevtransform(exp.(g), Val(1))
+        @test T \ exp.(x) == F \ exp.(x) == chebyshevtransform(exp.(g), Val(1))
         @test all(checkpoints(T) .∈ Ref(axes(T,1)))
 
         @test T \ [exp.(x) cos.(x)] ≈ [F \ exp.(x) F \ cos.(x)]