Skip to content

Commit 6397146

Browse files
authored
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
1 parent c3c0550 commit 6397146

File tree

8 files changed

+222
-51
lines changed

8 files changed

+222
-51
lines changed

Project.toml

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "ContinuumArrays"
22
uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
3-
version = "0.18.7"
3+
version = "0.19"
44

55
[deps]
66
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -31,14 +31,14 @@ ArrayLayouts = "1.0"
3131
BandedMatrices = "1"
3232
BlockArrays = "1"
3333
DomainSets = "0.6, 0.7"
34-
FastTransforms = "0.15, 0.16"
34+
FastTransforms = "0.15, 0.16, 0.17"
3535
FillArrays = "1.0"
36-
InfiniteArrays = "0.14, 0.15"
36+
InfiniteArrays = "0.15"
3737
Infinities = "0.1"
3838
IntervalSets = "0.7"
3939
LazyArrays = "2"
4040
Makie = "0.20, 0.21, 0.22"
41-
QuasiArrays = "0.11.8"
41+
QuasiArrays = "0.12"
4242
RecipesBase = "1.0"
4343
StaticArrays = "1.0"
4444
julia = "1.10"

docs/src/index.md

Lines changed: 24 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,12 +20,35 @@ product of the specified sizes, similar to `plan_fft`.
2020
5. `plotgrid(::MyBasis, n...)`: return `n` grid points suitable for plotting the basis. The default value for `n` is 10,000.
2121

2222

23-
## Routines
23+
## Differential Operators
2424

2525

2626
```@docs
2727
Derivative
2828
```
29+
30+
```@docs
31+
Laplacian
32+
```
33+
34+
```@docs
35+
AbsLaplacian
36+
```
37+
38+
```@docs
39+
laplacian
40+
```
41+
42+
```@docs
43+
abslaplacian
44+
```
45+
46+
```@docs
47+
weaklaplacian
48+
```
49+
50+
## Routines
51+
2952
```@docs
3053
transform
3154
```

src/ContinuumArrays.jl

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -19,11 +19,12 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
1919
ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
2020
LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize, _cutdim,
2121
AbstractQuasiFill, UnionDomain, sum_size, sum_layout, _cumsum, cumsum_layout, applylayout, equals_layout, layout_broadcasted, PolynomialLayout, dot_size,
22-
diff_layout, diff_size
22+
diff_layout, diff_size, AbstractQuasiVecOrMat
2323
import InfiniteArrays: Infinity, InfAxes
2424
import AbstractFFTs: Plan
2525

26-
export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform, basis, coefficients
26+
export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform, basis, coefficients,
27+
weaklaplacian, laplacian, Laplacian, AbsLaplacian, abslaplacian
2728

2829

2930

@@ -107,7 +108,7 @@ include("plotting.jl")
107108

108109
sum_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, dims) = _sum(expand(a), dims)
109110
dot_size(::InfiniteCardinal{1}, a, b) = dot(expand(a), expand(b))
110-
diff_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, dims) = diff(expand(a); dims=dims)
111+
diff_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, order...; dims...) = diff(expand(a), order...; dims...)
111112
function copy(d::Dot{<:ExpansionLayout,<:ExpansionLayout,<:AbstractQuasiArray,<:AbstractQuasiArray})
112113
a,b = d.A,d.B
113114
P,c = basis(a),coefficients(a)

src/bases/bases.jl

Lines changed: 79 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -421,6 +421,8 @@ basis_axes(ax, v) = error("Overload for $ax")
421421
coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any}}) where N = v.args[2]
422422
coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any,Vararg{Any}}}) where N = ApplyArray(*, tail(v.args)...)
423423

424+
coefficients(B::Basis{T}) where T = SquareEye{T}((axes(B,2),))
425+
424426

425427
function unweighted(lay::ExpansionLayout, a)
426428
wP,c = arguments(lay, a)
@@ -660,29 +662,45 @@ cumsum_layout(::ExpansionLayout, A, dims) = cumsum_layout(ApplyLayout{typeof(*)}
660662
###
661663
# diff
662664
###
665+
diff_layout(::AbstractBasisLayout, Vm, order...; dims...) = error("Overload diff(::$(typeof(Vm)))")
666+
function diff_layout(::AbstractBasisLayout, a, order::Int; dims...)
667+
order < 0 && throw(ArgumentError("order must be non-negative"))
668+
order == 0 && return a
669+
isone(order) ? diff(a) : diff(diff(a), order-1)
670+
end
663671

664-
diff_layout(::AbstractBasisLayout, Vm, dims...) = error("Overload diff(::$(typeof(Vm)))")
672+
function diff_layout(::SubBasisLayout, Vm, order::Int; dims::Integer=1)
673+
dims == 1 || error("not implemented")
674+
diff(parent(Vm), order)[:,parentindices(Vm)[2]]
675+
end
665676

666-
function diff_layout(::SubBasisLayout, Vm, dims::Integer=1)
677+
function diff_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
667678
dims == 1 || error("not implemented")
668-
diff(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
679+
diff(parent(Vm), order...)[:,parentindices(Vm)[2]]
669680
end
670681

671-
function diff_layout(::WeightedBasisLayout{SubBasisLayout}, Vm, dims::Integer=1)
682+
683+
function diff_layout(::WeightedBasisLayout{SubBasisLayout}, Vm, order...; dims::Integer=1)
672684
dims == 1 || error("not implemented")
673685
w = weight(Vm)
674686
V = unweighted(Vm)
675-
view(diff(w .* parent(V)), parentindices(V)...)
687+
view(diff(w .* parent(V), order...), parentindices(V)...)
676688
end
677689

678-
function diff_layout(::MappedBasisLayouts, V, dims::Integer=1)
679-
kr = basismap(V)
680-
@assert kr isa AbstractAffineQuasiVector
681-
D = diff(demap(V); dims=dims)
690+
diff_layout(::MappedBasisLayouts, V, order::Int; dims...) = diff_mapped(basismap(V), V, order::Int; dims...)
691+
diff_layout(::MappedBasisLayouts, V, order...; dims...) = diff_mapped(basismap(V), V, order...; dims...)
692+
693+
function diff_mapped(kr::AbstractAffineQuasiVector, V; dims...)
694+
D = diff(demap(V); dims...)
682695
view(basis(D), kr, :) * (kr.A*coefficients(D))
683696
end
684697

685-
diff_layout(::ExpansionLayout, A, dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, dims...)
698+
function diff_mapped(kr::AbstractAffineQuasiVector, V, order::Int; dims...)
699+
D = diff(demap(V), order; dims...)
700+
view(basis(D), kr, :) * (kr.A^order*coefficients(D))
701+
end
702+
703+
diff_layout(::ExpansionLayout, A, order...; dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
686704

687705

688706
####
@@ -708,6 +726,57 @@ function grammatrix_layout(::MappedBasisLayouts, P)
708726
grammatrix(Q)/kr.A
709727
end
710728

729+
"""
730+
abslaplacian(A, α=1)
731+
732+
computes ``|Δ|^α * A``.
733+
"""
734+
abslaplacian(A, order...; dims...) = abslaplacian_layout(MemoryLayout(A), A, order...; dims...)
735+
abslaplacian_layout(layout, A, order...; dims...) = abslaplacian_axis(axes(A,1), A, order...; dims...)
736+
abslaplacian_axis(::Inclusion{<:Number}, A, order=1; dims...) = -diff(A, 2order; dims...)
737+
738+
"""
739+
abslaplacian(A, k=1)
740+
741+
computes ``Δ^k * A``.
742+
"""
743+
laplacian(A, order...; dims...) = laplacian_layout(MemoryLayout(A), A, order...; dims...)
744+
laplacian_layout(layout, A, order...; dims...) = laplacian_axis(axes(A,1), A, order...; dims...)
745+
laplacian_axis(::Inclusion{<:Number}, A, order=1; dims...) = diff(A, 2order; dims...)
746+
747+
748+
laplacian_layout(::ExpansionLayout, A, order...; dims...) = laplacian_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
749+
abslaplacian_layout(::ExpansionLayout, A, order...; dims...) = abslaplacian_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
750+
751+
function abslaplacian_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
752+
dims == 1 || error("not implemented")
753+
abslaplacian(parent(Vm), order...)[:,parentindices(Vm)[2]]
754+
end
755+
756+
function laplacian_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
757+
dims == 1 || error("not implemented")
758+
laplacian(parent(Vm), order...)[:,parentindices(Vm)[2]]
759+
end
760+
761+
function laplacian_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, order...; dims=1)
762+
a = arguments(LAY, V)
763+
dims == 1 || throw(ArgumentError("cannot take laplacian a vector along dimension $dims"))
764+
*(laplacian(a[1], order...), tail(a)...)
765+
end
766+
767+
function abslaplacian_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, order...; dims=1)
768+
a = arguments(LAY, V)
769+
dims == 1 || throw(ArgumentError("cannot take abslaplacian a vector along dimension $dims"))
770+
*(abslaplacian(a[1], order...), tail(a)...)
771+
end
772+
773+
774+
775+
"""
776+
weaklaplacian(A)
777+
778+
represents the weak Laplacian.
779+
"""
711780
weaklaplacian(A) = weaklaplacian_layout(MemoryLayout(A), A)
712781
weaklaplacian_layout(_, A) = weaklaplacian_axis(axes(A,1), A)
713782
weaklaplacian_axis(::Inclusion{<:Number}, A) = -(diff(A)'diff(A))

src/bases/basisconcat.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8,9 +8,9 @@ abstract type AbstractConcatBasis{T} <: Basis{T} end
88

99
copy(S::AbstractConcatBasis) = S
1010

11-
function diff(S::AbstractConcatBasis; dims::Integer)
11+
function diff(S::AbstractConcatBasis, order...; dims::Integer=1)
1212
dims == 1 || error("not implemented")
13-
args = arguments.(Ref(ApplyLayout{typeof(*)}()), diff.(S.args; dims=dims))
13+
args = arguments.(Ref(ApplyLayout{typeof(*)}()), diff.(S.args, order...; dims=dims))
1414
all(length.(args) .== 2) || error("Not implemented")
1515
concatbasis(S, map(first, args)...) * mortar(Diagonal([map(last, args)...]))
1616
end
@@ -112,4 +112,4 @@ function QuasiArrays._getindex(::Type{IND}, A::HvcatBasis{T}, (x,j)::IND) where
112112
end
113113

114114

115-
diff(H::ApplyQuasiMatrix{<:Any,typeof(hcat)}; dims::Integer=1) = hcat((diff.(H.args; dims=dims))...)
115+
diff(H::ApplyQuasiMatrix{<:Any,typeof(hcat)}, order...; dims::Integer=1) = hcat((diff.(H.args, order...; dims=dims))...)

src/operators.jl

Lines changed: 79 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -108,51 +108,97 @@ show(io::IO, ::MIME"text/plain", δ::DiracDelta) = show(io, δ)
108108
# Differentiation
109109
#########
110110

111+
abstract type AbstractDifferentialQuasiMatrix{T} <: LazyQuasiMatrix{T} end
112+
111113
"""
112114
Derivative(axis)
113115
114116
represents the differentiation (or finite-differences) operator on the
115117
specified axis.
116118
"""
117-
struct Derivative{T,D} <: LazyQuasiMatrix{T}
118-
axis::Inclusion{T,D}
119+
struct Derivative{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
120+
axis::D
121+
order::Order
119122
end
120123

121-
Derivative{T}(axis::Inclusion{<:Any,D}) where {T,D} = Derivative{T,D}(axis)
122-
Derivative{T}(domain) where T = Derivative{T}(Inclusion(domain))
124+
"""
125+
Laplacian(axis)
126+
127+
represents the laplacian operator `Δ` on the
128+
specified axis.
129+
"""
130+
struct Laplacian{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
131+
axis::D
132+
order::Order
133+
end
123134

124-
Derivative(L::AbstractQuasiMatrix) = Derivative(axes(L,1))
135+
"""
136+
AbsLaplacian(axis)
125137
126-
show(io::IO, a::Derivative) = summary(io, a)
127-
function summary(io::IO, D::Derivative)
128-
print(io, "Derivative(")
129-
summary(io,D.axis)
130-
print(io,")")
138+
represents the positive-definite/negative/absolute-value laplacian operator `|Δ| ≡ -Δ` on the
139+
specified axis.
140+
"""
141+
struct AbsLaplacian{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
142+
axis::D
143+
order::Order
131144
end
132145

133-
axes(D::Derivative) = (D.axis, D.axis)
134-
==(a::Derivative, b::Derivative) = a.axis == b.axis
135-
copy(D::Derivative) = Derivative(copy(D.axis))
146+
operatororder(D) = something(D.order, 1)
147+
148+
show(io::IO, a::AbstractDifferentialQuasiMatrix) = summary(io, a)
149+
axes(D::AbstractDifferentialQuasiMatrix) = (D.axis, D.axis)
150+
==(a::AbstractDifferentialQuasiMatrix, b::AbstractDifferentialQuasiMatrix) = a.axis == b.axis && operatororder(a) == operatororder(b)
151+
copy(D::AbstractDifferentialQuasiMatrix) = D
152+
153+
136154

137-
@simplify function *(D::Derivative, B::AbstractQuasiMatrix)
155+
@simplify function *(D::AbstractDifferentialQuasiMatrix, B::AbstractQuasiVecOrMat)
138156
T = typeof(zero(eltype(D)) * zero(eltype(B)))
139-
diff(convert(AbstractQuasiMatrix{T}, B); dims=1)
157+
operatorcall(D, convert(AbstractQuasiArray{T}, B), D.order)
140158
end
141159

142-
@simplify function *(D::Derivative, B::AbstractQuasiVector)
143-
T = typeof(zero(eltype(D)) * zero(eltype(B)))
144-
diff(convert(AbstractQuasiVector{T}, B))
160+
^(D::AbstractDifferentialQuasiMatrix{T}, k::Integer) where T = similaroperator(D, D.axis, operatororder(D) .* k)
161+
162+
function view(D::AbstractDifferentialQuasiMatrix, kr::Inclusion, jr::Inclusion)
163+
@boundscheck axes(D,1) == kr == jr || throw(BoundsError(D,(kr,jr)))
164+
D
145165
end
146166

167+
operatorcall(D::AbstractDifferentialQuasiMatrix, B, order) = operatorcall(D)(B, order)
168+
operatorcall(D::AbstractDifferentialQuasiMatrix, B, ::Nothing) = operatorcall(D)(B)
147169

148170

171+
operatorcall(::Derivative) = diff
172+
operatorcall(::Laplacian) = laplacian
173+
operatorcall(::AbsLaplacian) = abslaplacian
149174

150-
^(D::Derivative, k::Integer) = ApplyQuasiArray(^, D, k)
151175

176+
for Op in (:Derivative, :Laplacian, :AbsLaplacian)
177+
@eval begin
178+
$Op{T, D}(axis::D, order) where {T,D<:Inclusion} = $Op{T,D,typeof(order)}(axis, order)
179+
$Op{T, D}(axis::D) where {T,D<:Inclusion} = $Op{T,D,Nothing}(axis, nothing)
180+
$Op{T}(axis::D, order...) where {T,D<:Inclusion} = $Op{T,D}(axis, order...)
181+
$Op{T}(domain, order...) where T = $Op{T}(Inclusion(domain), order...)
152182

153-
function view(D::Derivative, kr::Inclusion, jr::Inclusion)
154-
@boundscheck axes(D,1) == kr == jr || throw(BoundsError(D,(kr,jr)))
155-
D
183+
$Op(ax::AbstractQuasiVector{T}, order...) where T = $Op{eltype(eltype(ax))}(ax, order...)
184+
$Op(L::AbstractQuasiMatrix, order...) = $Op(axes(L,1), order...)
185+
186+
similaroperator(D::$Op, ax, ord) = $Op{eltype(D)}(ax, ord)
187+
188+
simplifiable(::typeof(*), A::$Op, B::$Op) = Val(true)
189+
*(D1::$Op, D2::$Op) = similaroperator(convert(AbstractQuasiMatrix{promote_type(eltype(D1),eltype(D2))}, D1), D1.axis, operatororder(D1)+operatororder(D2))
190+
191+
192+
function summary(io::IO, D::$Op)
193+
print(io, "$($Op)(")
194+
summary(io, D.axis)
195+
if !isnothing(D.order)
196+
print(io, ", ")
197+
print(io, D.order)
198+
end
199+
print(io,")")
200+
end
201+
end
156202
end
157203

158204
# struct Multiplication{T,F,A} <: AbstractQuasiMatrix{T}
@@ -161,11 +207,17 @@ end
161207
# end
162208

163209

164-
const Identity{T,D} = QuasiDiagonal{T,Inclusion{T,D}}
165-
166-
Identity(d::Inclusion) = QuasiDiagonal(d)
167-
168210
struct OperatorLayout <: AbstractLazyLayout end
169-
MemoryLayout(::Type{<:Derivative}) = OperatorLayout()
211+
MemoryLayout(::Type{<:AbstractDifferentialQuasiMatrix}) = OperatorLayout()
170212
# copy(M::Mul{OperatorLayout, <:ExpansionLayout}) = simplify(M)
171213
# copy(M::Mul{OperatorLayout, <:AbstractLazyLayout}) = M.A * expand(M.B)
214+
215+
216+
# Laplacian
217+
218+
abs::Laplacian{T}) where T = AbsLaplacian{T}(axes(Δ,1), Δ.order)
219+
-::Laplacian{<:Any,<:Any,Nothing}) = abs(Δ)
220+
-::AbsLaplacian{T,<:Any,Nothing}) where T = Laplacian{T}.axis)
221+
222+
^::AbsLaplacian{T}, k::Real) where T = AbsLaplacian{T}.axis, operatororder(Δ)*k)
223+
^::AbsLaplacian{T}, k::Integer) where T = AbsLaplacian{T}.axis, operatororder(Δ)*k)

test/runtests.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, App
3939
D = Derivative(x)
4040
@test D == Derivative{Float64}(x) == Derivative{Float64}(-1..1)
4141
@test D*x QuasiOnes(x)
42-
@test D^2 * x QuasiZeros(x)
42+
@test D^2 * x == zeros(x)
4343
@test D*[x D*x] == [D*x D^2*x]
4444
@test stringmime("text/plain", D) == "Derivative(Inclusion($(-1..1)))"
4545
@test_throws DimensionMismatch Derivative(Inclusion(0..1)) * x

0 commit comments

Comments
 (0)