Skip to content

Commit e67d37f

Browse files
dlfivefiftydlfivefifty
authored
Add AffineMap (#23)
* Add AffineMap * add find for isequal (which is vector-version of inverting a map) * Update runtests.jl * SubArray \, == for AffineMap * Change of variables works! * comine_mul * FillArrays 0.8, * Heat example * Update bases.jl * More general broadcast Layout * Update bases.jl * Update Project.toml * Update Project.toml * tests pass * Update ContinuumArrays.jl
1 parent 2b50508 commit e67d37f

File tree

7 files changed

+327
-61
lines changed

7 files changed

+327
-61
lines changed

‎Project.toml

+5-5
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.0.1"
3+
version = "0.0.2"
44

55
[deps]
66
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
@@ -11,11 +11,11 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
1111
QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
1212

1313
[compat]
14-
BandedMatrices = "0.12, 0.13"
15-
FillArrays = "0.7"
14+
BandedMatrices = "0.14"
15+
FillArrays = "0.8"
1616
IntervalSets = "0.3.1"
17-
LazyArrays = "0.12, 0.14"
18-
QuasiArrays = "0.0.2, 0.0.3, 0.0.4"
17+
LazyArrays = "0.14"
18+
QuasiArrays = "0.0.4"
1919
julia = "1.1"
2020

2121
[extras]

‎examples/heat.jl

+37
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
using ContinuumArrays, DifferentialEquations, Plots
2+
3+
##
4+
# natural boundary conditions
5+
##
6+
7+
L = LinearSpline(range(-1,1; length=20))
8+
x = axes(L,1)
9+
D = Derivative(x)
10+
M = L'L
11+
Δ = -((D*L)'D*L)
12+
u0 = copy(L \ exp.(x))
13+
14+
heat(u,(M,Δ),t) = M\(Δ*u)
15+
prob = ODEProblem(heat,u0,(0.0,1.0),(M,Δ))
16+
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
17+
18+
g = range(-1,1;length=1000)
19+
@gif for t in 0.0:0.01:1.0
20+
plot(g, L[g,:]*sol(t); ylims=(0,3))
21+
end
22+
23+
##
24+
# Dirichlet boundary conditions
25+
##
26+
27+
S = L[:,2:end-1]
28+
M = S'S
29+
Δ = -((D*S)'D*S)
30+
u0 = (L \ broadcast(x -> (1-x^2)*exp(x), x))[2:end-1]
31+
prob = ODEProblem(heat,u0,(0.0,1.0),(M,Δ))
32+
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
33+
34+
g = range(-1,1;length=1000)
35+
@gif for t in 0.0:0.01:1.0
36+
plot(g, S[g,:]*sol(t); ylims=(0,3))
37+
end

‎src/ContinuumArrays.jl

+74-22
Original file line numberDiff line numberDiff line change
@@ -1,20 +1,22 @@
11
module ContinuumArrays
22
using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays
3-
import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
3+
import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==,
44
IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
5-
first, last, show
6-
import Base.Broadcast: materialize, BroadcastStyle
7-
import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, adjointlayout, LdivApplyStyle
5+
first, last, show, isempty, findfirst, findlast, findall, Slice
6+
import Base.Broadcast: materialize, BroadcastStyle, broadcasted
7+
import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport,
8+
adjointlayout, LdivApplyStyle, arguments, broadcastlayout, lazy_getindex,
9+
sublayout, ApplyLayout, BroadcastLayout, combine_mul_styles
810
import LinearAlgebra: pinv
911
import BandedMatrices: AbstractBandedLayout, _BandedMatrix
1012
import FillArrays: AbstractFill, getindex_value
1113

1214
import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
13-
QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
15+
QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, quasimulapplystyle,
1416
ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
1517
LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle
1618

17-
export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis
19+
export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid
1820

1921
####
2022
# Interval indexing support
@@ -38,35 +40,85 @@ const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}
3840
cardinality(::AbstractInterval) = ℵ₁
3941
*(ℵ::AlephInfinity) = ℵ
4042

43+
Inclusion(d::AbstractInterval{T}) where T = Inclusion{float(T)}(d)
4144
first(S::Inclusion{<:Any,<:AbstractInterval}) = leftendpoint(S.domain)
4245
last(S::Inclusion{<:Any,<:AbstractInterval}) = rightendpoint(S.domain)
4346

47+
for find in (:findfirst, :findlast)
48+
@eval $find(f::Base.Fix2{typeof(isequal)}, d::Inclusion) = f.x in d.domain ? f.x : nothing
49+
end
4450

45-
checkindex(::Type{Bool}, inds::AbstractInterval, i::Number) = (leftendpoint(inds) <= i) & (i <= rightendpoint(inds))
46-
checkindex(::Type{Bool}, inds::AbstractInterval, i::Inclusion) = i.domain ⊆ inds
47-
function checkindex(::Type{Bool}, inds::AbstractInterval, I::AbstractArray)
48-
@_inline_meta
49-
b = true
50-
for i in I
51-
b &= checkindex(Bool, inds, i)
52-
end
53-
b
51+
function findall(f::Base.Fix2{typeof(isequal)}, d::Inclusion)
52+
r = findfirst(f,d)
53+
r === nothing ? eltype(d)[] : [r]
5454
end
5555

5656

57-
# we represent as a Mul with a banded matrix
58-
function materialize(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
59-
A = parent(V)
60-
_,jr = parentindices(V)
61-
first(jr) ≥ 1 || throw(BoundsError())
62-
P = _BandedMatrix(Ones{Int}(1,length(jr)), axes(A,2), first(jr)-1,1-first(jr))
63-
A*P
57+
function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval})
58+
@_propagate_inbounds_meta
59+
isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
6460
end
6561

62+
6663
BroadcastStyle(::Type{<:Inclusion{<:Any,<:AbstractInterval}}) = LazyQuasiArrayStyle{1}()
6764
BroadcastStyle(::Type{<:QuasiAdjoint{<:Any,<:Inclusion{<:Any,<:AbstractInterval}}}) = LazyQuasiArrayStyle{2}()
6865
BroadcastStyle(::Type{<:QuasiTranspose{<:Any,<:Inclusion{<:Any,<:AbstractInterval}}}) = LazyQuasiArrayStyle{2}()
6966

67+
68+
###
69+
# Maps
70+
###
71+
72+
# Affine map represents A*x .+ b
73+
struct AffineQuasiVector{T,AA,X,B} <: AbstractQuasiVector{T}
74+
A::AA
75+
x::X
76+
b::B
77+
end
78+
79+
AffineQuasiVector(A::AA, x::X, b::B) where {AA,X,B} =
80+
AffineQuasiVector{promote_type(eltype(AA), eltype(X), eltype(B)),AA,X,B}(A,x,b)
81+
82+
AffineQuasiVector(A, x) = AffineQuasiVector(A, x, zero(promote_type(eltype(A),eltype(x))))
83+
AffineQuasiVector(x) = AffineQuasiVector(one(eltype(x)), x)
84+
85+
AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*x.b .+ b)
86+
87+
axes(A::AffineQuasiVector) = axes(A.x)
88+
getindex(A::AffineQuasiVector, k::Number) = A.A*A.x[k] .+ A.b
89+
inbounds_getindex(A::AffineQuasiVector{<:Any,<:Any,<:Inclusion}, k::Number) = A.A*k .+ A.b
90+
isempty(A::AffineQuasiVector) = isempty(A.x)
91+
==(a::AffineQuasiVector, b::AffineQuasiVector) = a.A == b.A && a.x == b.x && a.b == b.b
92+
93+
BroadcastStyle(::Type{<:AffineQuasiVector}) = LazyQuasiArrayStyle{1}()
94+
95+
for op in(:*, :\, :+, :-)
96+
@eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), a::Number, x::Inclusion) = broadcast($op, a, AffineQuasiVector(x))
97+
end
98+
for op in(:/, :+, :-)
99+
@eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Inclusion, a::Number) = broadcast($op, AffineQuasiVector(x), a)
100+
end
101+
102+
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), a::Number, x::AffineQuasiVector) = AffineQuasiVector(a, x)
103+
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(\), a::Number, x::AffineQuasiVector) = AffineQuasiVector(inv(a), x)
104+
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(/), x::AffineQuasiVector, a::Number) = AffineQuasiVector(inv(a), x)
105+
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), a::Number, x::AffineQuasiVector) = AffineQuasiVector(one(eltype(x)), x, a)
106+
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), x::AffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, b)
107+
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), a::Number, x::AffineQuasiVector) = AffineQuasiVector(-one(eltype(x)), x, a)
108+
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), x::AffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, -b)
109+
110+
function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AffineQuasiVector{<:Real,<:Real,<:Inclusion{<:Real,<:AbstractInterval}})
111+
@_propagate_inbounds_meta
112+
isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
113+
end
114+
115+
for find in (:findfirst, :findlast, :findall)
116+
@eval $find(f::Base.Fix2{typeof(isequal)}, d::AffineQuasiVector) = $find(isequal(d.A \ (f.x .- d.b)), d.x)
117+
end
118+
119+
120+
121+
70122
include("operators.jl")
71123
include("bases/bases.jl")
72124

‎src/bases/bases.jl

+126-7
Original file line numberDiff line numberDiff line change
@@ -1,27 +1,146 @@
11
abstract type Basis{T} <: LazyQuasiMatrix{T} end
2+
abstract type Weight{T} <: LazyQuasiVector{T} end
23

34

45
const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}
56

6-
MemoryLayout(::Type{<:Basis}) = LazyLayout()
7+
struct WeightLayout <: MemoryLayout end
8+
struct BasisLayout <: MemoryLayout end
9+
struct AdjointBasisLayout <: MemoryLayout end
10+
11+
MemoryLayout(::Type{<:Basis}) = BasisLayout()
12+
MemoryLayout(::Type{<:Weight}) = WeightLayout()
13+
14+
adjointlayout(::Type, ::BasisLayout) = AdjointBasisLayout()
15+
transposelayout(::Type{<:Real}, ::BasisLayout) = AdjointBasisLayout()
16+
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = BasisLayout()
17+
18+
combine_mul_styles(::BasisLayout) = LazyQuasiArrayApplyStyle()
19+
combine_mul_styles(::AdjointBasisLayout) = LazyQuasiArrayApplyStyle()
20+
721
ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
822
pinv(J::Basis) = apply(pinv,J)
923

1024
_multup(a::Tuple) = Mul(a...)
1125
_multup(a) = a
1226

1327

14-
==(A::Basis, B::Basis) = axes(A) ≠ axes(B) ||
15-
throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
28+
function ==(A::Basis, B::Basis)
29+
axes(A) == axes(B) && throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
30+
false
31+
end
32+
1633

1734
ApplyStyle(::typeof(\), ::Type{<:Basis}, ::Type{<:AbstractQuasiMatrix}) = LdivApplyStyle()
1835
ApplyStyle(::typeof(\), ::Type{<:Basis}, ::Type{<:AbstractQuasiVector}) = LdivApplyStyle()
36+
ApplyStyle(::typeof(\), ::Type{<:SubQuasiArray{<:Any,2,<:Basis}}, ::Type{<:AbstractQuasiMatrix}) = LdivApplyStyle()
37+
ApplyStyle(::typeof(\), ::Type{<:SubQuasiArray{<:Any,2,<:Basis}}, ::Type{<:AbstractQuasiVector}) = LdivApplyStyle()
1938

20-
function copy(P::Ldiv{<:Any,<:Any,<:Basis,<:Basis})
21-
A, B = P.A, P.B
22-
A == B || throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
23-
Eye(size(A,2))
39+
copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
40+
function copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(-)}})
41+
a,b = arguments(L.B)
42+
(L.A\a)-(L.A\b)
43+
end
44+
45+
for Bas1 in (:Basis, :WeightedBasis), Bas2 in (:Basis, :WeightedBasis)
46+
@eval begin
47+
function copy(P::Ldiv{<:Any,<:Any,<:$Bas1,<:$Bas2})
48+
A, B = P.A, P.B
49+
A == B || throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
50+
Eye(size(A,2))
51+
end
52+
function copy(P::Ldiv{<:Any,<:Any,<:SubQuasiArray{<:Any,2,<:$Bas1},<:SubQuasiArray{<:Any,2,<:$Bas2}})
53+
A, B = P.A, P.B
54+
(parent(A) == parent(B) && parentindices(A) == parentindices(B)) ||
55+
throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
56+
Eye(size(A,2))
57+
end
58+
59+
function copy(P::Ldiv{<:Any,<:Any,<:SubQuasiArray{<:Any,2,<:$Bas1,<:Tuple{<:AffineQuasiVector,<:Slice}},
60+
<:SubQuasiArray{<:Any,2,<:$Bas2,<:Tuple{<:AffineQuasiVector,<:Slice}}})
61+
A, B = P.A, P.B
62+
parent(A)\parent(B)
63+
end
64+
function copy(P::Ldiv{<:Any,<:Any,<:SubQuasiArray{<:Any,2,<:$Bas1,<:Tuple{<:AffineQuasiVector,<:Slice}},
65+
<:SubQuasiArray{<:Any,2,<:$Bas2,<:Tuple{<:AffineQuasiVector,<:Any}}})
66+
A, B = P.A, P.B
67+
# use lazy_getindex to avoid sparse arrays
68+
lazy_getindex(parent(A)\parent(B),:,parentindices(B)[2])
69+
end
70+
71+
function ==(A::SubQuasiArray{<:Any,2,<:$Bas1}, B::SubQuasiArray{<:Any,2,<:$Bas2})
72+
all(parentindices(A) == parentindices(B)) && parent(A) == parent(B)
73+
end
74+
end
75+
end
76+
77+
78+
# expansion
79+
grid(P) = error("Overload Grid")
80+
function transform(L)
81+
p = grid(L)
82+
p,L[p,:]
83+
end
84+
85+
function copy(L::Ldiv{BasisLayout,<:Any,<:Any,<:AbstractQuasiVector})
86+
p,T = transform(L.A)
87+
T \ L.B[p]
88+
end
89+
90+
copy(L::Ldiv{BasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
91+
copy(Ldiv{LazyLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
92+
93+
function copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix,<:AbstractQuasiVector})
94+
p,T = transform(L.A)
95+
T \ L.B[p]
96+
end
97+
98+
## materialize views
99+
100+
# materialize(S::SubQuasiArray{<:Any,2,<:ApplyQuasiArray{<:Any,2,typeof(*),<:Tuple{<:Basis,<:Any}}}) =
101+
# *(arguments(S)...)
102+
103+
104+
# Differentiation of sub-arrays
105+
function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}})
106+
A, B = M.args
107+
P = parent(B)
108+
(Derivative(axes(P,1))*P)[parentindices(B)...]
109+
end
110+
111+
function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AffineQuasiVector,<:Any}}})
112+
A, B = M.args
113+
P = parent(B)
114+
kr,jr = parentindices(B)
115+
(Derivative(axes(P,1))*P*kr.A)[kr,jr]
116+
end
117+
118+
function copy(L::Ldiv{BasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix})
119+
args = arguments(L.B)
120+
# this is a temporary hack
121+
if args isa Tuple{AbstractQuasiMatrix,Number}
122+
(L.A \ first(args))*last(args)
123+
elseif args isa Tuple{Number,AbstractQuasiMatrix}
124+
first(args)*(L.A \ last(args))
125+
else
126+
error("Not implemented")
127+
end
128+
end
129+
130+
131+
# we represent as a Mul with a banded matrix
132+
sublayout(::BasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = ApplyLayout{typeof(*)}()
133+
sublayout(::BasisLayout, ::Type{<:Tuple{<:AffineQuasiVector,<:AbstractUnitRange}}) = BasisLayout()
134+
function arguments(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
135+
A = parent(V)
136+
_,jr = parentindices(V)
137+
first(jr) ≥ 1 || throw(BoundsError())
138+
P = _BandedMatrix(Ones{Int}(1,length(jr)), axes(A,2), first(jr)-1,1-first(jr))
139+
A,P
24140
end
25141

26142

27143
include("splines.jl")
144+
145+
146+

‎src/bases/splines.jl

+8-4
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,10 @@
1-
struct Spline{order,T} <: Basis{T}
2-
points::Vector{T}
1+
struct Spline{order,T,P<:AbstractVector} <: Basis{T}
2+
points::P
33
end
4+
Spline{o,T}(pts::P) where {o,T,P} = Spline{o,T,P}(pts)
45

5-
const LinearSpline{T} = Spline{1,T}
6-
const HeavisideSpline{T} = Spline{0,T}
6+
const LinearSpline = Spline{1}
7+
const HeavisideSpline = Spline{0}
78

89
Spline{o}(pts::AbstractVector{T}) where {o,T} = Spline{o,float(T)}(pts)
910

@@ -36,6 +37,9 @@ function getindex(B::HeavisideSpline{T}, x::Number, k::Int) where T
3637
return zero(T)
3738
end
3839

40+
grid(L::LinearSpline) = L.points
41+
transform(L::LinearSpline{T}) where T = grid(L),Eye{T}(size(L,2))
42+
3943
## Sub-bases
4044

4145

‎src/operators.jl

-6
Original file line numberDiff line numberDiff line change
@@ -154,12 +154,6 @@ axes(D::Derivative) = (D.axis, D.axis)
154154
copy(D::Derivative) = Derivative(copy(D.axis))
155155

156156

157-
function copy(M::QMul2{<:Derivative,<:SubQuasiArray})
158-
A, B = M.args
159-
P = parent(B)
160-
(Derivative(axes(P,1))*P)[parentindices(B)...]
161-
end
162-
163157

164158
# struct Multiplication{T,F,A} <: AbstractQuasiMatrix{T}
165159
# f::F

0 commit comments

Comments
 (0)