Skip to content

Commit 212481f

Browse files
authored
Support InfiniteArrays v0.9 (#69)
* Support InfniteArrays v0.9 * Improve mapping * move find to QuasiArrays.jl * Add QuadMap test * fix find, add in for map * Update runtests.jl
1 parent 322b156 commit 212481f

File tree

4 files changed

+126
-54
lines changed

4 files changed

+126
-54
lines changed

Project.toml

Lines changed: 2 additions & 2 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.3.6"
3+
version = "0.4.0"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -19,7 +19,7 @@ FillArrays = "0.9.3, 0.10, 0.11"
1919
InfiniteArrays = "0.8, 0.9"
2020
IntervalSets = "0.4, 0.5"
2121
LazyArrays = "0.19, 0.20"
22-
QuasiArrays = "0.3.5"
22+
QuasiArrays = "0.3.8"
2323
julia = "1.5"
2424

2525
[extras]

src/ContinuumArrays.jl

Lines changed: 41 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -90,16 +90,6 @@ function dot(x::Inclusion{T,<:AbstractInterval}, y::Inclusion{V,<:AbstractInterv
9090
end
9191

9292

93-
for find in (:findfirst, :findlast)
94-
@eval $find(f::Base.Fix2{typeof(isequal)}, d::Inclusion) = f.x in d.domain ? f.x : nothing
95-
end
96-
97-
function findall(f::Base.Fix2{typeof(isequal)}, d::Inclusion)
98-
r = findfirst(f,d)
99-
r === nothing ? eltype(d)[] : [r]
100-
end
101-
102-
10393
function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval})
10494
@_propagate_inbounds_meta
10595
isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
@@ -110,8 +100,38 @@ end
110100
# Maps
111101
###
112102

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+
113133
# Affine map represents A*x .+ b
114-
abstract type AbstractAffineQuasiVector{T,AA,X,B} <: AbstractQuasiVector{T} end
134+
abstract type AbstractAffineQuasiVector{T,AA,X,B} <: Map{T} end
115135

116136
summary(io::IO, a::AbstractAffineQuasiVector) = print(io, "$(a.A) * $(a.x) .+ ($(a.b))")
117137

@@ -130,7 +150,9 @@ AffineQuasiVector(x) = AffineQuasiVector(one(eltype(x)), x)
130150
AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*x.b .+ b)
131151

132152
axes(A::AbstractAffineQuasiVector) = axes(A.x)
153+
133154
affine_getindex(A, k) = A.A*A.x[k] .+ A.b
155+
Base.unsafe_getindex(A::AbstractAffineQuasiVector, k) = A.A*Base.unsafe_getindex(A.x,k) .+ A.b
134156
getindex(A::AbstractAffineQuasiVector, k::Number) = affine_getindex(A, k)
135157
function getindex(A::AbstractAffineQuasiVector, k::Inclusion)
136158
@boundscheck A.x[k] # throws bounds error if k ≠ x
@@ -167,17 +189,14 @@ function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::
167189
isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
168190
end
169191

170-
affine_igetindex(d, x) = d.A \ (x .- d.b)
171-
igetindex(d::AbstractAffineQuasiVector, x) = affine_igetindex(d, x)
172-
173-
for find in (:findfirst, :findlast, :findall)
174-
@eval $find(f::Base.Fix2{typeof(isequal)}, d::AbstractAffineQuasiVector) = $find(isequal(igetindex(d, f.x)), d.x)
175-
end
176-
177192
minimum(d::AbstractAffineQuasiVector) = signbit(d.A) ? last(d) : first(d)
178193
maximum(d::AbstractAffineQuasiVector) = signbit(d.A) ? first(d) : last(d)
179194

180195
union(d::AbstractAffineQuasiVector) = Inclusion(minimum(d)..maximum(d))
196+
invmap(d::AbstractAffineQuasiVector) = affine(union(d), axes(d,1))
197+
198+
199+
181200

182201

183202
struct AffineMap{T,D,R} <: AbstractAffineQuasiVector{T,T,D,T}
@@ -191,9 +210,10 @@ AffineMap(domain::AbstractQuasiVector{T}, range::AbstractQuasiVector{V}) where {
191210
measure(x::Inclusion) = last(x)-first(x)
192211

193212
function getproperty(A::AffineMap, d::Symbol)
194-
d == :x && return A.domain
195-
d == :A && return measure(A.range)/measure(A.domain)
196-
d == :b && return (last(A.domain)*first(A.range) - first(A.domain)*last(A.range))/measure(A.domain)
213+
domain, range = getfield(A, :domain), getfield(A, :range)
214+
d == :x && return domain
215+
d == :A && return measure(range)/measure(domain)
216+
d == :b && return (last(domain)*first(range) - first(domain)*last(range))/measure(domain)
197217
getfield(A, d)
198218
end
199219

@@ -204,12 +224,6 @@ function getindex(A::AffineMap, k::Number)
204224
affine_getindex(A, k)
205225
end
206226

207-
function igetindex(A::AffineMap, k::Number)
208-
# ensure we exactly hit range
209-
k == first(A.range) && return first(A.domain)
210-
k == last(A.range) && return last(A.domain)
211-
affine_igetindex(A, k)
212-
end
213227

214228
first(A::AffineMap) = first(A.range)
215229
last(A::AffineMap) = last(A.range)

src/bases/bases.jl

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,8 @@ broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::MappedBasisLayouts) = Mappe
3434

3535
# A sub of a weight is still a weight
3636
sublayout(::WeightLayout, _) = WeightLayout()
37+
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{Map,AbstractUnitRange}}) = MappedBasisLayout()
38+
3739

3840
## Weighted basis interface
3941
unweightedbasis(P::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{AbstractQuasiVector,AbstractQuasiMatrix}}) = last(P.args)
@@ -112,7 +114,7 @@ copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:Any,<:AbstractQu
112114

113115
# expansion
114116
_grid(_, P) = error("Overload Grid")
115-
_grid(::MappedBasisLayout, P) = igetindex.(Ref(parentindices(P)[1]), grid(demap(P)))
117+
_grid(::MappedBasisLayout, P) = invmap(parentindices(P)[1])[grid(demap(P))]
116118
_grid(::SubBasisLayout, P) = grid(parent(P))
117119
_grid(::WeightedBasisLayouts, P) = grid(unweightedbasis(P))
118120
grid(P) = _grid(MemoryLayout(typeof(P)), P)
@@ -152,11 +154,22 @@ end
152154
\(a::ProjectionFactorization, b::AbstractVector) = (a.F \ b)[a.inds]
153155

154156
_factorize(::SubBasisLayout, L) = ProjectionFactorization(factorize(parent(L)), parentindices(L)[2])
155-
# function _factorize(::MappedBasisLayout, L)
156-
# kr, jr = parentindices(L)
157-
# P = parent(L)
158-
# ProjectionFactorization(factorize(view(P,:,jr)), parentindices(L)[2])
159-
# end
157+
158+
struct MappedFactorization{T, FAC<:Factorization{T}, MAP} <: Factorization{T}
159+
F::FAC
160+
map::MAP
161+
end
162+
163+
\(a::MappedFactorization, b::AbstractQuasiVector) = a.F \ view(b, a.map)
164+
\(a::MappedFactorization, b::AbstractVector) = a.F \ b
165+
166+
function invmap end
167+
168+
function _factorize(::MappedBasisLayout, L)
169+
kr, jr = parentindices(L)
170+
P = parent(L)
171+
MappedFactorization(factorize(view(P,:,jr)), invmap(parentindices(L)[1]))
172+
end
160173

161174
transform_ldiv(A, B, _) = factorize(A) \ B
162175
transform_ldiv(A, B) = transform_ldiv(A, B, size(A))

test/runtests.jl

Lines changed: 64 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, FastTransforms, InfiniteArrays, Test, Base64
22
import ContinuumArrays: ℵ₁, materialize, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, ℵ₁,
3-
MappedBasisLayout, AdjointMappedBasisLayout, MappedWeightedBasisLayout, igetindex, TransformFactorization, Weight, WeightedBasisLayout, SubWeightedBasisLayout, WeightLayout,
4-
Expansion, basis
3+
MappedBasisLayout, AdjointMappedBasisLayout, MappedWeightedBasisLayout, TransformFactorization, Weight, WeightedBasisLayout, SubWeightedBasisLayout, WeightLayout,
4+
Expansion, basis, invmap, Map
55
import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
66
import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout, LdivStyle, MulStyle
77

@@ -24,8 +24,9 @@ end
2424

2525
@testset "Inclusion" begin
2626
x = Inclusion(-1..1)
27-
@test x[0.1] === 0.1
28-
@test x[0.0] === 0.0
27+
@test eltype(x) == Float64
28+
@test x[0.1] 0.1
29+
@test x[0] x[0.0] 0.0
2930

3031
x = Inclusion(-1.0..1)
3132
X = QuasiDiagonal(x)
@@ -100,9 +101,9 @@ end
100101
@test_throws BoundsError a[-3]
101102
@test a[-2] == first(a) == -1
102103
@test a[3] == last(a) == 1
103-
@test igetindex(a,-0.16) 0.1
104-
@test igetindex(a,1) == 3
105-
@test igetindex(a,-1) == -2
104+
@test invmap(a)[-0.16] 0.1
105+
@test invmap(a)[1] == 3
106+
@test invmap(a)[-1] == -2
106107
@test union(a) == Inclusion(-1..1)
107108

108109
@test affine(0..1, -1..1) == y
@@ -116,11 +117,11 @@ end
116117

117118
@testset "DiracDelta" begin
118119
δ = DiracDelta(-1..3)
119-
@test axes(δ) === (axes(δ,1),) === (Inclusion(-1..3),)
120-
@test size(δ) === (length(δ),) === (ℵ₁,)
121-
@test δ[1.1] === 0.0
122-
@test δ[0.0] === Inf
123-
@test Base.IndexStyle(δ) === Base.IndexLinear()
120+
@test axes(δ) (axes(δ,1),) (Inclusion(-1..3),)
121+
@test size(δ) (length(δ),) (ℵ₁,)
122+
@test δ[1.1] 0.0
123+
@test δ[0.0] Inf
124+
@test Base.IndexStyle(δ) Base.IndexLinear()
124125

125126
@test stringmime("text/plain", δ) == "δ at 0.0 over Inclusion(-1..3)"
126127
x = Inclusion(-1..1)
@@ -131,14 +132,14 @@ end
131132
@testset "HeavisideSpline" begin
132133
H = HeavisideSpline([1,2,3])
133134

134-
@test axes(H) === (axes(H,1),axes(H,2)) === (Inclusion(1..3), Base.OneTo(2))
135-
@test size(H) === (size(H,1),size(H,2)) === (ℵ₁, 2)
135+
@test axes(H) (axes(H,1),axes(H,2)) (Inclusion(1..3), Base.OneTo(2))
136+
@test size(H) (size(H,1),size(H,2)) (ℵ₁, 2)
136137

137138
@test_throws BoundsError H[0.1, 1]
138-
@test H[1.1,1] === H'[1,1.1] === transpose(H)[1,1.1] === 1.0
139-
@test H[2.1,1] === H'[1,2.1] === transpose(H)[1,2.1] === 0.0
140-
@test H[1.1,2] === 0.0
141-
@test H[2.1,2] === 1.0
139+
@test H[1.1,1] H'[1,1.1] transpose(H)[1,1.1] 1.0
140+
@test H[2.1,1] H'[1,2.1] transpose(H)[1,2.1] 0.0
141+
@test H[1.1,2] 0.0
142+
@test H[2.1,2] 1.0
142143
@test_throws BoundsError H[2.1,3]
143144
@test_throws BoundsError H'[3,2.1]
144145
@test_throws BoundsError transpose(H)[3,2.1]
@@ -483,6 +484,12 @@ end
483484
@test_throws BoundsError K[Inclusion(0..0.5), Inclusion(0..0.5)][1,1]
484485
end
485486

487+
"""
488+
This is a simple implementation of Chebyshev for testing. Use OrthogonalPolynomialsQuasi
489+
for the real implementation.
490+
"""
491+
492+
486493
struct Chebyshev <: Basis{Float64}
487494
n::Int
488495
end
@@ -503,6 +510,20 @@ LinearAlgebra.factorize(L::Chebyshev) =
503510
# This is wrong but just for tests
504511
Base.broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:Chebyshev}, b::Chebyshev) = b * Matrix(I, 5, 5)
505512

513+
struct QuadraticMap{T} <: Map{T} end
514+
struct InvQuadraticMap{T} <: Map{T} end
515+
516+
QuadraticMap() = QuadraticMap{Float64}()
517+
InvQuadraticMap() = InvQuadraticMap{Float64}()
518+
519+
Base.getindex(::QuadraticMap, r::Number) = 2r^2-1
520+
Base.axes(::QuadraticMap{T}) where T = (Inclusion(0..1),)
521+
Base.axes(::InvQuadraticMap{T}) where T = (Inclusion(-1..1),)
522+
Base.getindex(d::InvQuadraticMap, x::Number) = sqrt((x+1)/2)
523+
ContinuumArrays.invmap(::QuadraticMap{T}) where T = InvQuadraticMap{T}()
524+
ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
525+
526+
506527
@testset "Chebyshev" begin
507528
T = Chebyshev(5)
508529
w = ChebyshevWeight()
@@ -536,9 +557,33 @@ Base.broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:Che
536557
@test wT[y,:][[0.1,0.2],1:5] == (w[y] .* T[y,:])[[0.1,0.2],1:5] == (w .* T[:,1:5])[y,:][[0.1,0.2],:]
537558
@test MemoryLayout(wT[y,1:3]) isa MappedWeightedBasisLayout
538559
@test wT[y,1:3][[0.1,0.2],1:2] == wT[y[[0.1,0.2]],1:2]
560+
561+
@testset "QuadraticMap" begin
562+
m = QuadraticMap()
563+
mi = InvQuadraticMap()
564+
@test 0.1 m
565+
@test -0.1 m
566+
@test 2 m
567+
@test 0.1 mi
568+
@test -0.1 mi
569+
570+
@test m[findfirst(isequal(0.1), m)] 0.1
571+
@test m[findlast(isequal(0.1), m)] 0.1
572+
@test m[findall(isequal(0.1), m)] [0.1]
573+
574+
T = Chebyshev(5)
575+
M = T[m,:]
576+
@test MemoryLayout(M) isa MappedBasisLayout
577+
@test M[0.1,:] T[2*0.1^2-1,:]
578+
x = axes(M,1)
579+
@test x == Inclusion(0..1)
580+
@test M \ exp.(x) T \ exp.(sqrt.((axes(T,1) .+ 1)/2))
581+
end
539582
end
540583

541584
@testset "Broadcasted" begin
585+
T = Chebyshev(5)
586+
x = axes(T,1)
542587
a = 1 .+ x .+ x.^2
543588
# The following are wrong, just testing dispatch
544589
@test T \ (a .* T) == I
@@ -549,4 +594,4 @@ Base.broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:Che
549594
= T * (T \ a)
550595
@test T \ (ã .* ã) [1.5,1,0.5,0,0]
551596
end
552-
end
597+
end

0 commit comments

Comments
 (0)