Skip to content

Support order in Derivative #197

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 12 commits into from
Jan 29, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -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"
Expand Down Expand Up @@ -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"
Expand Down
25 changes: 24 additions & 1 deletion docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
```
Expand Down
7 changes: 4 additions & 3 deletions src/ContinuumArrays.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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



Expand Down Expand Up @@ -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)
Expand Down
89 changes: 79 additions & 10 deletions src/bases/bases.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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...)


####
Expand All @@ -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))
Expand Down
6 changes: 3 additions & 3 deletions src/bases/basisconcat.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))...)
106 changes: 79 additions & 27 deletions src/operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -108,51 +108,97 @@
# 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)

Check warning on line 155 in src/operators.jl

View check run for this annotation

Codecov / codecov/patch

src/operators.jl#L155

Added line #L155 was not covered by tests
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}
Expand All @@ -161,11 +207,17 @@
# 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)
2 changes: 1 addition & 1 deletion test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading
Loading