Skip to content

Commit 5e8cf76

Browse files
authored
Support basis concatenation (#75)
* Create basisconcat.jl * move over block code, add Piecewise/VcatBlockDiagonal * Fix checkindex * Add Derivative for concat arrays * remove AbstractInclusion * ADd HvcatBasis * Update runtests.jl * Update Project.toml * increase coverge
1 parent d4dd514 commit 5e8cf76

File tree

9 files changed

+458
-247
lines changed

9 files changed

+458
-247
lines changed

Project.toml

Lines changed: 12 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,29 @@
11
name = "ContinuumArrays"
22
uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
3-
version = "0.4.2"
3+
version = "0.5.0"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
77
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
8+
BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
89
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
910
InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
1011
IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
1112
LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
1213
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
1314
QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
15+
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
1416

1517
[compat]
16-
ArrayLayouts = "0.4.10, 0.5"
17-
BandedMatrices = "0.15.17, 0.16"
18-
FillArrays = "0.9.3, 0.10, 0.11"
19-
InfiniteArrays = "0.8, 0.9"
20-
IntervalSets = "0.4, 0.5"
21-
LazyArrays = "0.19, 0.20"
22-
QuasiArrays = "0.3.8"
18+
ArrayLayouts = "0.5"
19+
BandedMatrices = "0.16"
20+
BlockArrays = "0.14"
21+
FillArrays = "0.11"
22+
InfiniteArrays = "0.9"
23+
IntervalSets = "0.5"
24+
LazyArrays = "0.20"
25+
QuasiArrays = "0.4.1"
26+
StaticArrays = "0.12, 1"
2327
julia = "1.5"
2428

2529
[extras]

src/ContinuumArrays.jl

Lines changed: 31 additions & 164 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,26 @@
11
module ContinuumArrays
2-
using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays, InfiniteArrays
2+
using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays, InfiniteArrays, StaticArrays, BlockArrays
33
import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==, ^,
4-
IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, diff,
4+
IndexStyle, IndexLinear, ==, OneTo, _maybetail, tail, similar, copyto!, copy, diff,
55
first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum, sum, _sum,
6-
getproperty, isone, iszero, zero, abs, <, , >, , string, summary
6+
getproperty, isone, iszero, zero, abs, <, , >, , string, summary, to_indices
77
import Base.Broadcast: materialize, BroadcastStyle, broadcasted
88
import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, most, combine_mul_styles, AbstractArrayApplyStyle,
99
adjointlayout, arguments, _mul_arguments, call, broadcastlayout, layout_getindex, UnknownLayout,
1010
sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles, applylayout,
1111
simplifiable, _simplify
1212
import LinearAlgebra: pinv, dot, norm2
1313
import BandedMatrices: AbstractBandedLayout, _BandedMatrix
14+
import BlockArrays: block, blockindex, unblock, blockedrange, _BlockedUnitRange, _BlockArray
1415
import FillArrays: AbstractFill, getindex_value, SquareEye
1516
import ArrayLayouts: mul
1617
import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
1718
QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
1819
ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
1920
LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize
20-
import InfiniteArrays: Infinity
21+
import InfiniteArrays: Infinity, InfAxes
2122

22-
export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform, affine
23+
export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform, affine, ..
2324

2425
####
2526
# Interval indexing support
@@ -90,177 +91,43 @@ function dot(x::Inclusion{T,<:AbstractInterval}, y::Inclusion{V,<:AbstractInterv
9091
end
9192

9293

93-
function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval})
94-
@_propagate_inbounds_meta
95-
isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
96-
end
97-
98-
99-
###
100-
# Maps
101-
###
102-
103-
"""
104-
A subtype of `Map` is used as a one-to-one map between two domains
105-
via `view`. The domain of the map `m` is `axes(m,1)` and the range
106-
is `union(m)`.
107-
108-
Maps must also overload `invmap` to give the inverse of the map, which
109-
is equivalent to `invmap(m)[x] == findfirst(isequal(x), m)`.
110-
"""
111-
112-
abstract type Map{T} <: AbstractQuasiVector{T} end
113-
114-
invmap(M::Map) = error("Overload invmap(::$(typeof(M)))")
115-
116-
117-
Base.in(x, m::Map) = x in union(m)
118-
Base.issubset(d::Map, b::IntervalSets.Domain) = union(d) b
119-
Base.union(d::Map) = axes(invmap(d),1)
120-
121-
for find in (:findfirst, :findlast)
122-
@eval function $find(f::Base.Fix2{typeof(isequal)}, d::Map)
123-
f.x in d || return nothing
124-
$find(isequal(invmap(d)[f.x]), union(d))
125-
end
126-
end
127-
128-
@eval function findall(f::Base.Fix2{typeof(isequal)}, d::Map)
129-
f.x in d || return eltype(axes(d,1))[]
130-
findall(isequal(invmap(d)[f.x]), union(d))
131-
end
132-
133-
function Base.getindex(d::Map, x::Inclusion)
134-
x == axes(d,1) || throw(BoundsError(d, x))
135-
d
136-
end
137-
138-
# Affine map represents A*x .+ b
139-
abstract type AbstractAffineQuasiVector{T,AA,X,B} <: Map{T} end
140-
141-
summary(io::IO, a::AbstractAffineQuasiVector) = print(io, "$(a.A) * $(a.x) .+ ($(a.b))")
142-
143-
struct AffineQuasiVector{T,AA,X,B} <: AbstractAffineQuasiVector{T,AA,X,B}
144-
A::AA
145-
x::X
146-
b::B
147-
end
148-
149-
AffineQuasiVector(A::AA, x::X, b::B) where {AA,X,B} =
150-
AffineQuasiVector{promote_type(eltype(AA), eltype(X), eltype(B)),AA,X,B}(A,x,b)
151-
152-
AffineQuasiVector(A, x) = AffineQuasiVector(A, x, zero(promote_type(eltype(A),eltype(x))))
153-
AffineQuasiVector(x) = AffineQuasiVector(one(eltype(x)), x)
154-
155-
AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*x.b .+ b)
156-
157-
axes(A::AbstractAffineQuasiVector) = axes(A.x)
158-
159-
affine_getindex(A, k) = A.A*A.x[k] .+ A.b
160-
Base.unsafe_getindex(A::AbstractAffineQuasiVector, k) = A.A*Base.unsafe_getindex(A.x,k) .+ A.b
161-
getindex(A::AbstractAffineQuasiVector, k::Number) = affine_getindex(A, k)
162-
function getindex(A::AbstractAffineQuasiVector, k::Inclusion)
163-
@boundscheck A.x[k] # throws bounds error if k ≠ x
164-
A
165-
end
166-
167-
getindex(A::AbstractAffineQuasiVector, ::Colon) = copy(A)
168-
169-
copy(A::AbstractAffineQuasiVector) = A
170-
171-
inbounds_getindex(A::AbstractAffineQuasiVector{<:Any,<:Any,<:Inclusion}, k::Number) = A.A*k .+ A.b
172-
isempty(A::AbstractAffineQuasiVector) = isempty(A.x)
173-
==(a::AbstractAffineQuasiVector, b::AbstractAffineQuasiVector) = a.A == b.A && a.x == b.x && a.b == b.b
174-
175-
BroadcastStyle(::Type{<:AbstractAffineQuasiVector}) = LazyQuasiArrayStyle{1}()
176-
177-
for op in(:*, :\, :+, :-)
178-
@eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), a::Number, x::Inclusion) = broadcast($op, a, AffineQuasiVector(x))
179-
end
180-
for op in(:/, :+, :-)
181-
@eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Inclusion, a::Number) = broadcast($op, AffineQuasiVector(x), a)
182-
end
183-
184-
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(a, x)
185-
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(\), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(inv(a), x)
186-
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(/), x::AbstractAffineQuasiVector, a::Number) = AffineQuasiVector(inv(a), x)
187-
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(one(eltype(x)), x, a)
188-
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, b)
189-
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(-one(eltype(x)), x, a)
190-
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, -b)
191-
192-
function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AbstractAffineQuasiVector)
193-
@_propagate_inbounds_meta
194-
isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
195-
end
196-
197-
minimum(d::AbstractAffineQuasiVector) = signbit(d.A) ? last(d) : first(d)
198-
maximum(d::AbstractAffineQuasiVector) = signbit(d.A) ? first(d) : last(d)
199-
200-
union(d::AbstractAffineQuasiVector) = Inclusion(minimum(d)..maximum(d))
201-
invmap(d::AbstractAffineQuasiVector) = affine(union(d), axes(d,1))
202-
203-
204-
94+
include("maps.jl")
20595

96+
const QInfAxes = Union{Inclusion,AbstractAffineQuasiVector}
20697

207-
struct AffineMap{T,D,R} <: AbstractAffineQuasiVector{T,T,D,T}
208-
domain::D
209-
range::R
210-
end
21198

212-
AffineMap(domain::AbstractQuasiVector{T}, range::AbstractQuasiVector{V}) where {T,V} =
213-
AffineMap{promote_type(T,V), typeof(domain),typeof(range)}(domain,range)
99+
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes}) = V
100+
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,QInfAxes}) = V
101+
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{Any,QInfAxes}) = V
102+
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,Any}) = V
214103

215-
measure(x::Inclusion) = last(x)-first(x)
104+
# ambiguity error
105+
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{InfAxes,QInfAxes}) = V
106+
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,InfAxes}) = V
216107

217-
function getproperty(A::AffineMap, d::Symbol)
218-
domain, range = getfield(A, :domain), getfield(A, :range)
219-
d == :x && return domain
220-
d == :A && return measure(range)/measure(domain)
221-
d == :b && return (last(domain)*first(range) - first(domain)*last(range))/measure(domain)
222-
getfield(A, d)
223-
end
108+
#
109+
# BlockQuasiArrays
224110

225-
function getindex(A::AffineMap, k::Number)
226-
# ensure we exactly hit range
227-
k == first(A.domain) && return first(A.range)
228-
k == last(A.domain) && return last(A.range)
229-
affine_getindex(A, k)
111+
BlockArrays.blockaxes(::Inclusion) = blockaxes(Base.OneTo(1)) # just use 1 block
112+
function BlockArrays.blockaxes(A::AbstractQuasiArray{T,N}, d) where {T,N}
113+
@_inline_meta
114+
d::Integer <= N ? blockaxes(A)[d] : Base.OneTo(1)
230115
end
231116

117+
@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{Block{1}, Vararg{Any}}) =
118+
(unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
119+
@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{BlockRange{1,R}, Vararg{Any}}) where R =
120+
(unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
121+
@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{BlockIndex{1}, Vararg{Any}}) =
122+
(inds[1][I[1]], to_indices(A, _maybetail(inds), tail(I))...)
232123

233-
first(A::AffineMap) = first(A.range)
234-
last(A::AffineMap) = last(A.range)
235-
236-
affine(a::AbstractQuasiVector, b::AbstractQuasiVector) = AffineMap(a, b)
237-
affine(a, b::AbstractQuasiVector) = affine(Inclusion(a), b)
238-
affine(a::AbstractQuasiVector, b) = affine(a, Inclusion(b))
239-
affine(a, b) = affine(Inclusion(a), Inclusion(b))
240-
241-
242-
# mapped vectors
243-
const AffineMappedQuasiVector = SubQuasiArray{<:Any, 1, <:Any, <:Tuple{AbstractAffineQuasiVector}}
244-
const AffineMappedQuasiMatrix = SubQuasiArray{<:Any, 2, <:Any, <:Tuple{AbstractAffineQuasiVector,Slice}}
245-
246-
==(a::AffineMappedQuasiVector, b::AffineMappedQuasiVector) = parentindices(a) == parentindices(b) && parent(a) == parent(b)
247-
248-
_sum(V::AffineMappedQuasiVector, ::Colon) = parentindices(V)[1].A \ sum(parent(V))
249-
250-
# pretty print for bases
251-
summary(io::IO, P::AffineMappedQuasiMatrix) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
252-
summary(io::IO, P::AffineMappedQuasiVector) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
253-
254-
const QInfAxes = Union{Inclusion,AbstractAffineQuasiVector}
255-
256-
257-
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes}) = V
258-
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,QInfAxes}) = V
259-
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{Any,QInfAxes}) = V
260-
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,Any}) = V
124+
checkpoints(d::AbstractInterval{T}) where T = width(d) .* SVector{3,float(T)}(0.823972,0.01,0.3273484) .+ leftendpoint(d)
125+
checkpoints(x::Inclusion) = checkpoints(x.domain)
126+
checkpoints(A::AbstractQuasiMatrix) = checkpoints(axes(A,1))
261127

262128

263129
include("operators.jl")
264130
include("bases/bases.jl")
131+
include("basisconcat.jl")
265132

266133
end

src/bases/bases.jl

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -69,13 +69,13 @@ end
6969
@inline function copy(P::Ldiv{<:AbstractBasisLayout,<:AbstractBasisLayout})
7070
A, B = P.A, P.B
7171
A == B || throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
72-
SquareEye{eltype(P)}((axes(A,2),))
72+
SquareEye{eltype(eltype(P))}((axes(A,2),)) # use double eltype for array-valued
7373
end
7474
@inline function copy(P::Ldiv{<:SubBasisLayouts,<:SubBasisLayouts})
7575
A, B = P.A, P.B
7676
parent(A) == parent(B) ||
7777
throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
78-
Eye{eltype(P)}((axes(A,2),axes(B,2)))
78+
Eye{eltype(eltype(P))}((axes(A,2),axes(B,2)))
7979
end
8080

8181
@inline function copy(P::Ldiv{<:MappedBasisLayouts,<:MappedBasisLayouts})
@@ -119,8 +119,6 @@ _grid(::SubBasisLayout, P) = grid(parent(P))
119119
_grid(::WeightedBasisLayouts, P) = grid(unweightedbasis(P))
120120
grid(P) = _grid(MemoryLayout(typeof(P)), P)
121121

122-
# TODO: Move over from OrthogonalPolynomialsQuasi
123-
function checkpoints end
124122

125123
struct TransformFactorization{T,Grid,Plan,IPlan} <: Factorization{T}
126124
grid::Grid

0 commit comments

Comments
 (0)