Support order in Derivative (#197)

* Support order in Derivative * Update diff and add weak/abs/laplacian * tests pass * Update bases.jl * unify operators def * Work on abslaplacian * Improve operators * Update Project.toml * add tests and docs * Update test_splines.jl * Increase coverage * add test
JuliaApproximation · Jan 29, 2025 · 6397146 · 6397146 · dlfivefifty · Jan 29, 2025
1 parent c3c0550
commit 6397146
Show file tree

Hide file tree

Showing 8 changed files with 222 additions and 51 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.18.7"
+version = "0.19"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -31,14 +31,14 @@ ArrayLayouts = "1.0"
 BandedMatrices = "1"
 BlockArrays = "1"
 DomainSets = "0.6, 0.7"
-FastTransforms = "0.15, 0.16"
+FastTransforms = "0.15, 0.16, 0.17"
 FillArrays = "1.0"
-InfiniteArrays = "0.14, 0.15"
+InfiniteArrays = "0.15"
 Infinities = "0.1"
 IntervalSets = "0.7"
 LazyArrays = "2"
 Makie = "0.20, 0.21, 0.22"
-QuasiArrays = "0.11.8"
+QuasiArrays = "0.12"
 RecipesBase = "1.0"
 StaticArrays = "1.0"
 julia = "1.10"

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -20,12 +20,35 @@ product of the specified sizes, similar to `plan_fft`.
 5. `plotgrid(::MyBasis, n...)`: return `n` grid points suitable for plotting the basis. The default value for `n` is 10,000. 
 
 
-## Routines
+## Differential Operators
 
 
 ```@docs
 Derivative
 ```
+
+```@docs
+Laplacian
+```
+
+```@docs
+AbsLaplacian
+```
+
+```@docs
+laplacian
+```
+
+```@docs
+abslaplacian
+```
+
+```@docs
+weaklaplacian
+```
+
+## Routines
+
 ```@docs
 transform
 ```

diff --git a/src/ContinuumArrays.jl b/src/ContinuumArrays.jl
@@ -19,11 +19,12 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize, _cutdim,
                     AbstractQuasiFill, UnionDomain, sum_size, sum_layout, _cumsum, cumsum_layout, applylayout, equals_layout, layout_broadcasted, PolynomialLayout, dot_size,
-                    diff_layout, diff_size
+                    diff_layout, diff_size, AbstractQuasiVecOrMat
 import InfiniteArrays: Infinity, InfAxes
 import AbstractFFTs: Plan
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform, basis, coefficients 
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform, basis, coefficients,
+        weaklaplacian, laplacian, Laplacian, AbsLaplacian, abslaplacian
 
 
 
@@ -107,7 +108,7 @@ include("plotting.jl")
 
 sum_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, dims) = _sum(expand(a), dims)
 dot_size(::InfiniteCardinal{1}, a, b) = dot(expand(a), expand(b))
-diff_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, dims) = diff(expand(a); dims=dims)
+diff_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, order...; dims...) = diff(expand(a), order...; dims...)
 function copy(d::Dot{<:ExpansionLayout,<:ExpansionLayout,<:AbstractQuasiArray,<:AbstractQuasiArray})
     a,b = d.A,d.B
     P,c = basis(a),coefficients(a)

diff --git a/src/bases/bases.jl b/src/bases/bases.jl
@@ -421,6 +421,8 @@ basis_axes(ax, v) = error("Overload for $ax")
 coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any}}) where N = v.args[2]
 coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any,Vararg{Any}}}) where N = ApplyArray(*, tail(v.args)...)
 
+coefficients(B::Basis{T}) where T = SquareEye{T}((axes(B,2),))
+
 
 function unweighted(lay::ExpansionLayout, a)
     wP,c = arguments(lay, a)
@@ -660,29 +662,45 @@ cumsum_layout(::ExpansionLayout, A, dims) = cumsum_layout(ApplyLayout{typeof(*)}
 ###
 # diff
 ###
+diff_layout(::AbstractBasisLayout, Vm, order...; dims...) = error("Overload diff(::$(typeof(Vm)))")
+function diff_layout(::AbstractBasisLayout, a, order::Int; dims...)
+    order < 0 && throw(ArgumentError("order must be non-negative"))
+    order == 0 && return a
+    isone(order) ? diff(a) : diff(diff(a), order-1)
+end
 
-diff_layout(::AbstractBasisLayout, Vm, dims...) = error("Overload diff(::$(typeof(Vm)))")
+function diff_layout(::SubBasisLayout, Vm, order::Int; dims::Integer=1)
+    dims == 1 || error("not implemented")
+    diff(parent(Vm), order)[:,parentindices(Vm)[2]]
+end
 
-function diff_layout(::SubBasisLayout, Vm, dims::Integer=1)
+function diff_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
     dims == 1 || error("not implemented")
-    diff(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
+    diff(parent(Vm), order...)[:,parentindices(Vm)[2]]
 end
 
-function diff_layout(::WeightedBasisLayout{SubBasisLayout}, Vm, dims::Integer=1)
+
+function diff_layout(::WeightedBasisLayout{SubBasisLayout}, Vm, order...; dims::Integer=1)
     dims == 1 || error("not implemented")
     w = weight(Vm)
     V = unweighted(Vm)
-    view(diff(w .* parent(V)), parentindices(V)...)
+    view(diff(w .* parent(V), order...), parentindices(V)...)
 end
 
-function diff_layout(::MappedBasisLayouts, V, dims::Integer=1)
-    kr = basismap(V)
-    @assert kr isa AbstractAffineQuasiVector
-    D = diff(demap(V); dims=dims)
+diff_layout(::MappedBasisLayouts, V, order::Int; dims...) = diff_mapped(basismap(V), V, order::Int; dims...)
+diff_layout(::MappedBasisLayouts, V, order...; dims...) = diff_mapped(basismap(V), V, order...; dims...)
+
+function diff_mapped(kr::AbstractAffineQuasiVector, V; dims...)
+    D = diff(demap(V); dims...)
     view(basis(D), kr, :) * (kr.A*coefficients(D))
 end
 
-diff_layout(::ExpansionLayout, A, dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, dims...)
+function diff_mapped(kr::AbstractAffineQuasiVector, V, order::Int; dims...)
+    D = diff(demap(V), order; dims...)
+    view(basis(D), kr, :) * (kr.A^order*coefficients(D))
+end
+
+diff_layout(::ExpansionLayout, A, order...; dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
 
 
 ####
@@ -708,6 +726,57 @@ function grammatrix_layout(::MappedBasisLayouts, P)
     grammatrix(Q)/kr.A
 end
 
+"""
+    abslaplacian(A, α=1)
+
+computes ``|Δ|^α * A``. 
+"""
+abslaplacian(A, order...; dims...) = abslaplacian_layout(MemoryLayout(A), A, order...; dims...)
+abslaplacian_layout(layout, A, order...; dims...) = abslaplacian_axis(axes(A,1), A, order...; dims...)
+abslaplacian_axis(::Inclusion{<:Number}, A, order=1; dims...) = -diff(A, 2order; dims...)
+
+"""
+    abslaplacian(A, k=1)
+
+computes ``Δ^k * A``. 
+"""
+laplacian(A, order...; dims...) = laplacian_layout(MemoryLayout(A), A, order...; dims...)
+laplacian_layout(layout, A, order...; dims...) = laplacian_axis(axes(A,1), A, order...; dims...)
+laplacian_axis(::Inclusion{<:Number}, A, order=1; dims...) = diff(A, 2order; dims...)
+
+
+laplacian_layout(::ExpansionLayout, A, order...; dims...) = laplacian_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
+abslaplacian_layout(::ExpansionLayout, A, order...; dims...) = abslaplacian_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
+
+function abslaplacian_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
+    dims == 1 || error("not implemented")
+    abslaplacian(parent(Vm), order...)[:,parentindices(Vm)[2]]
+end
+
+function laplacian_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
+    dims == 1 || error("not implemented")
+    laplacian(parent(Vm), order...)[:,parentindices(Vm)[2]]
+end
+
+function laplacian_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, order...; dims=1)
+    a = arguments(LAY, V)
+    dims == 1 || throw(ArgumentError("cannot take laplacian a vector along dimension $dims"))
+    *(laplacian(a[1], order...), tail(a)...)
+end
+
+function abslaplacian_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, order...; dims=1)
+    a = arguments(LAY, V)
+    dims == 1 || throw(ArgumentError("cannot take abslaplacian a vector along dimension $dims"))
+    *(abslaplacian(a[1], order...), tail(a)...)
+end
+
+
+
+"""
+    weaklaplacian(A)
+
+represents the weak Laplacian.
+"""
 weaklaplacian(A) = weaklaplacian_layout(MemoryLayout(A), A)
 weaklaplacian_layout(_, A) = weaklaplacian_axis(axes(A,1), A)
 weaklaplacian_axis(::Inclusion{<:Number}, A) = -(diff(A)'diff(A))

diff --git a/src/bases/basisconcat.jl b/src/bases/basisconcat.jl
@@ -8,9 +8,9 @@ abstract type AbstractConcatBasis{T} <: Basis{T} end
 
 copy(S::AbstractConcatBasis) = S
 
-function diff(S::AbstractConcatBasis; dims::Integer)
+function diff(S::AbstractConcatBasis, order...; dims::Integer=1)
     dims == 1 || error("not implemented")
-    args = arguments.(Ref(ApplyLayout{typeof(*)}()), diff.(S.args; dims=dims))
+    args = arguments.(Ref(ApplyLayout{typeof(*)}()), diff.(S.args, order...; dims=dims))
     all(length.(args) .== 2) || error("Not implemented")
     concatbasis(S, map(first, args)...) * mortar(Diagonal([map(last, args)...]))
 end
@@ -112,4 +112,4 @@ function QuasiArrays._getindex(::Type{IND}, A::HvcatBasis{T}, (x,j)::IND) where
 end
 
 
-diff(H::ApplyQuasiMatrix{<:Any,typeof(hcat)}; dims::Integer=1) = hcat((diff.(H.args; dims=dims))...)
+diff(H::ApplyQuasiMatrix{<:Any,typeof(hcat)}, order...; dims::Integer=1) = hcat((diff.(H.args, order...; dims=dims))...)
diff --git a/src/operators.jl b/src/operators.jl
@@ -108,51 +108,97 @@ show(io::IO, ::MIME"text/plain", δ::DiracDelta) = show(io, δ)
 # Differentiation
 #########
 
+abstract type AbstractDifferentialQuasiMatrix{T} <: LazyQuasiMatrix{T} end
+
 """
 Derivative(axis)
 
 represents the differentiation (or finite-differences) operator on the
 specified axis.
 """
-struct Derivative{T,D} <: LazyQuasiMatrix{T}
-    axis::Inclusion{T,D}
+struct Derivative{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
+    axis::D
+    order::Order
 end
 
-Derivative{T}(axis::Inclusion{<:Any,D}) where {T,D} = Derivative{T,D}(axis)
-Derivative{T}(domain) where T = Derivative{T}(Inclusion(domain))
+"""
+Laplacian(axis)
+
+represents the laplacian operator `Δ` on the
+specified axis.
+"""
+struct Laplacian{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
+    axis::D
+    order::Order
+end
 
-Derivative(L::AbstractQuasiMatrix) = Derivative(axes(L,1))
+"""
+AbsLaplacian(axis)
 
-show(io::IO, a::Derivative) = summary(io, a)
-function summary(io::IO, D::Derivative)
-    print(io, "Derivative(")
-    summary(io,D.axis)
-    print(io,")")
+represents the positive-definite/negative/absolute-value laplacian operator `|Δ| ≡ -Δ` on the
+specified axis.
+"""
+struct AbsLaplacian{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
+    axis::D
+    order::Order
 end
 
-axes(D::Derivative) = (D.axis, D.axis)
-==(a::Derivative, b::Derivative) = a.axis == b.axis
-copy(D::Derivative) = Derivative(copy(D.axis))
+operatororder(D) = something(D.order, 1)
+
+show(io::IO, a::AbstractDifferentialQuasiMatrix) = summary(io, a)
+axes(D::AbstractDifferentialQuasiMatrix) = (D.axis, D.axis)
+==(a::AbstractDifferentialQuasiMatrix, b::AbstractDifferentialQuasiMatrix) = a.axis == b.axis && operatororder(a) == operatororder(b)
+copy(D::AbstractDifferentialQuasiMatrix) = D
+
+
 
-@simplify function *(D::Derivative, B::AbstractQuasiMatrix)
+@simplify function *(D::AbstractDifferentialQuasiMatrix, B::AbstractQuasiVecOrMat)
     T = typeof(zero(eltype(D)) * zero(eltype(B)))
-    diff(convert(AbstractQuasiMatrix{T}, B); dims=1)
+    operatorcall(D, convert(AbstractQuasiArray{T}, B), D.order)
 end
 
-@simplify function *(D::Derivative, B::AbstractQuasiVector)
-    T = typeof(zero(eltype(D)) * zero(eltype(B)))
-    diff(convert(AbstractQuasiVector{T}, B))
+^(D::AbstractDifferentialQuasiMatrix{T}, k::Integer) where T = similaroperator(D, D.axis, operatororder(D) .* k)
+
+function view(D::AbstractDifferentialQuasiMatrix, kr::Inclusion, jr::Inclusion)
+    @boundscheck axes(D,1) == kr == jr || throw(BoundsError(D,(kr,jr)))
+    D
 end
 
+operatorcall(D::AbstractDifferentialQuasiMatrix, B, order) = operatorcall(D)(B, order)
+operatorcall(D::AbstractDifferentialQuasiMatrix, B, ::Nothing) = operatorcall(D)(B)
 
 
+operatorcall(::Derivative) = diff
+operatorcall(::Laplacian) = laplacian
+operatorcall(::AbsLaplacian) = abslaplacian
 
-^(D::Derivative, k::Integer) = ApplyQuasiArray(^, D, k)
 
+for Op in (:Derivative, :Laplacian, :AbsLaplacian)
+    @eval begin
+        $Op{T, D}(axis::D, order) where {T,D<:Inclusion} = $Op{T,D,typeof(order)}(axis, order)
+        $Op{T, D}(axis::D) where {T,D<:Inclusion} = $Op{T,D,Nothing}(axis, nothing)
+        $Op{T}(axis::D, order...) where {T,D<:Inclusion} = $Op{T,D}(axis, order...)
+        $Op{T}(domain, order...) where T = $Op{T}(Inclusion(domain), order...)
 
-function view(D::Derivative, kr::Inclusion, jr::Inclusion)
-    @boundscheck axes(D,1) == kr == jr || throw(BoundsError(D,(kr,jr)))
-    D
+        $Op(ax::AbstractQuasiVector{T}, order...) where T = $Op{eltype(eltype(ax))}(ax, order...)
+        $Op(L::AbstractQuasiMatrix, order...) = $Op(axes(L,1), order...)
+
+        similaroperator(D::$Op, ax, ord) = $Op{eltype(D)}(ax, ord)
+
+        simplifiable(::typeof(*), A::$Op, B::$Op) = Val(true)
+        *(D1::$Op, D2::$Op) = similaroperator(convert(AbstractQuasiMatrix{promote_type(eltype(D1),eltype(D2))}, D1), D1.axis, operatororder(D1)+operatororder(D2))
+
+
+        function summary(io::IO, D::$Op)
+            print(io, "$($Op)(")
+            summary(io, D.axis)
+            if !isnothing(D.order)
+                print(io, ", ")
+                print(io, D.order)
+            end
+            print(io,")")
+        end
+    end
 end
 
 # struct Multiplication{T,F,A} <: AbstractQuasiMatrix{T}
@@ -161,11 +207,17 @@ end
 # end
 
 
-const Identity{T,D} = QuasiDiagonal{T,Inclusion{T,D}}
-
-Identity(d::Inclusion) = QuasiDiagonal(d)
-
 struct OperatorLayout <: AbstractLazyLayout end
-MemoryLayout(::Type{<:Derivative}) = OperatorLayout()
+MemoryLayout(::Type{<:AbstractDifferentialQuasiMatrix}) = OperatorLayout()
 # copy(M::Mul{OperatorLayout, <:ExpansionLayout}) = simplify(M)
 # copy(M::Mul{OperatorLayout, <:AbstractLazyLayout}) = M.A * expand(M.B)
+
+
+# Laplacian
+
+abs(Δ::Laplacian{T}) where T = AbsLaplacian{T}(axes(Δ,1), Δ.order)
+-(Δ::Laplacian{<:Any,<:Any,Nothing}) = abs(Δ)
+-(Δ::AbsLaplacian{T,<:Any,Nothing}) where T = Laplacian{T}(Δ.axis)
+
+^(Δ::AbsLaplacian{T}, k::Real) where T = AbsLaplacian{T}(Δ.axis, operatororder(Δ)*k)
+^(Δ::AbsLaplacian{T}, k::Integer) where T = AbsLaplacian{T}(Δ.axis, operatororder(Δ)*k)
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -39,7 +39,7 @@ import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, App
         D = Derivative(x)
         @test D == Derivative{Float64}(x) == Derivative{Float64}(-1..1)
         @test D*x ≡ QuasiOnes(x)
-        @test D^2 * x ≡ QuasiZeros(x)
+        @test D^2 * x == zeros(x)
         @test D*[x D*x] == [D*x D^2*x]
         @test stringmime("text/plain", D) == "Derivative(Inclusion($(-1..1)))"
         @test_throws DimensionMismatch Derivative(Inclusion(0..1)) * x