Skip to content

Commit f26a3d3

Browse files
dlfivefiftydlfivefifty
authored
Make Legendre transform type stable (#152)
* Make Legendre transform type stable * transform_ldiv -> transform_ldiv_size * tests pass * Update Project.toml * increase cov * Update test_normalized.jl * Update ci.yml
1 parent 6f1ce10 commit f26a3d3

File tree

13 files changed

+91
-34
lines changed

13 files changed

+91
-34
lines changed

‎.github/workflows/ci.yml‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,17 @@
11
name: CI
22
on:
3-
- push
4-
- pull_request
3+
push:
4+
branches:
5+
- main
6+
paths-ignore:
7+
- 'LICENSE'
8+
- 'README.md'
9+
- '.github/workflows/TagBot.yml'
10+
pull_request:
11+
paths-ignore:
12+
- 'LICENSE'
13+
- 'README.md'
14+
- '.github/workflows/TagBot.yml'
515
jobs:
616
test:
717
name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}

‎Project.toml‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -25,11 +25,11 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
2525
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
2626

2727
[compat]
28-
ArrayLayouts = "1.0.1"
28+
ArrayLayouts = "1.3.1"
2929
BandedMatrices = "0.17.17"
3030
BlockArrays = "0.16.9"
3131
BlockBandedMatrices = "0.12"
32-
ContinuumArrays = "0.15.2"
32+
ContinuumArrays = "0.16"
3333
DomainSets = "0.6"
3434
FFTW = "1.1"
3535
FastGaussQuadrature = "0.5"

‎src/ClassicalOrthogonalPolynomials.jl‎

Lines changed: 22 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,18 @@ MemoryLayout(::Type{<:OrthogonalPolynomial}) = OPLayout()
7575

7676
sublayout(::AbstractOPLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:Slice}}) = MappedOPLayout()
7777

78+
"""
79+
MappedOPLayout
80+
81+
represents an OP that is (usually affine) mapped OP
82+
"""
7883
struct MappedOPLayout <: AbstractOPLayout end
84+
85+
"""
86+
WeightedOPLayout
87+
88+
represents an OP multiplied by its orthogonality weight.
89+
"""
7990
struct WeightedOPLayout{Lay<:AbstractOPLayout} <: AbstractWeightedBasisLayout end
8091

8192
isorthogonalityweighted(::WeightedOPLayout, _) = true
@@ -96,10 +107,6 @@ broadcastbasis_layout(::typeof(+), ::MappedOPLayout, M::MappedBasisLayout, P, Q)
96107
broadcastbasis_layout(::typeof(+), L::MappedBasisLayout, ::MappedOPLayout, P, Q) = broadcastbasis_layout(+, L, MappedBasisLayout(), P, Q)
97108
sum_layout(::MappedOPLayout, A, dims) = sum_layout(MappedBasisLayout(), A, dims)
98109

99-
# demap to avoid Golub-Welsch fallback
100-
ContinuumArrays.transform_ldiv_if_columns(L::Ldiv{MappedOPLayout,Lay}, ax::OneTo) where Lay = ContinuumArrays.transform_ldiv_if_columns(Ldiv{MappedBasisLayout,Lay}(L.A,L.B), ax)
101-
ContinuumArrays.transform_ldiv_if_columns(L::Ldiv{MappedOPLayout,ApplyLayout{typeof(hcat)}}, ax::OneTo) = ContinuumArrays.transform_ldiv_if_columns(Ldiv{MappedBasisLayout,UnknownLayout}(L.A,L.B), ax)
102-
103110
equals_layout(::AbstractOPLayout, ::AbstractWeightedBasisLayout, _, _) = false # Weighted-Legendre doesn't exist
104111
equals_layout(::AbstractWeightedBasisLayout, ::AbstractOPLayout, _, _) = false # Weighted-Legendre doesn't exist
105112

@@ -109,7 +116,14 @@ equals_layout(::WeightedBasisLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP
109116
equals_layout(::WeightedBasisLayout{<:AbstractOPLayout}, ::WeightedBasisLayout{<:AbstractOPLayout}, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
110117

111118

112-
copy(L::Ldiv{MappedOPLayout,Lay}) where Lay<:MappedBasisLayouts = copy(Ldiv{MappedBasisLayout,Lay}(L.A,L.B))
119+
copy(L::Ldiv{MappedOPLayout,Lay}) where Lay = copy(Ldiv{MappedBasisLayout,Lay}(L.A,L.B))
120+
copy(L::Ldiv{MappedOPLayout,Lay}) where Lay<:AbstractLazyLayout = copy(Ldiv{MappedBasisLayout,Lay}(L.A,L.B))
121+
copy(L::Ldiv{MappedOPLayout,Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{MappedBasisLayout,Lay}(L.A,L.B))
122+
copy(L::Ldiv{MappedOPLayout,BroadcastLayout{typeof(-)}}) = copy(Ldiv{MappedBasisLayout,BroadcastLayout{typeof(-)}}(L.A,L.B))
123+
copy(L::Ldiv{MappedOPLayout,BroadcastLayout{typeof(+)}}) = copy(Ldiv{MappedBasisLayout,BroadcastLayout{typeof(+)}}(L.A,L.B))
124+
copy(L::Ldiv{MappedOPLayout,BroadcastLayout{typeof(*)}}) = copy(Ldiv{MappedBasisLayout,BroadcastLayout{typeof(*)}}(L.A,L.B))
125+
copy(L::Ldiv{MappedOPLayout,ApplyLayout{typeof(hcat)}}) = copy(Ldiv{MappedBasisLayout,ApplyLayout{typeof(hcat)}}(L.A,L.B))
126+
copy(L::Ldiv{MappedOPLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{MappedBasisLayout,ApplyLayout{typeof(*)}}(L.A,L.B))
113127

114128
# OPs are immutable
115129
copy(a::OrthogonalPolynomial) = a
@@ -168,6 +182,7 @@ singularities(w) = singularities(MemoryLayout(w), w)
168182
singularities(::ExpansionLayout, f) = singularities(basis(f))
169183

170184
singularitiesview(w, ::Inclusion) = w # for now just assume it doesn't change
185+
singularitiesview(w, ind) = view(w, ind)
171186
singularities(S::SubQuasiArray) = singularitiesview(singularities(parent(S)), parentindices(S)[1])
172187

173188
basis_axes(::Inclusion{<:Any,<:AbstractInterval}, v) = convert(AbstractQuasiMatrix{eltype(v)}, basis_singularities(singularities(v)))
@@ -291,14 +306,14 @@ end
291306
plan_grid_transform(::MappedOPLayout, L, szs::NTuple{N,Int}, dims=1:N) where N =
292307
plan_grid_transform(MappedBasisLayout(), L, szs, dims)
293308

294-
function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial})
309+
@simplify function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial})
295310
axes(A,1) == axes(B,1) || throw(DimensionMismatch())
296311
_,jA = parentindices(A)
297312
_,jB = parentindices(B)
298313
(parent(A) \ parent(B))[jA, jB]
299314
end
300315

301-
function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}})
316+
@simplify function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}})
302317
axes(A,1) == axes(B,1) || throw(DimensionMismatch())
303318
parent(A) \ parent(B)
304319
end

‎src/adaptivetransform.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11

2-
transform_ldiv(A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}, ::Tuple{<:Any,InfiniteCardinal{0}}) where {T,V} =
2+
ContinuumArrays.transform_ldiv_size(::Tuple{<:Any,InfiniteCardinal{0}}, A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}) where {T,V} =
33
adaptivetransform_ldiv(A, f)
44

55

@@ -31,7 +31,7 @@ function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiVector
3131

3232
for jr in increasingtruncations(ax)
3333
An = A[:,jr]
34-
cfs = An \ f
34+
cfs = transform_ldiv(An, f)
3535
maxabsc = maximum(abs, cfs)
3636
if maxabsc ≤ tol && maxabsfr ≤ tol # probably zero
3737
return pad(similar(cfs,0), ax)

‎src/classical/jacobi.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -177,9 +177,9 @@ plotgrid(P::AbstractJacobi{T}, n::Integer) where T = ChebyshevGrid{2,T}(min(40n,
177177

178178
ldiv(P::Jacobi{V}, f::Inclusion{T}) where {T,V} = _op_ldiv(P, f)
179179
ldiv(P::Jacobi{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)
180-
function ldiv(P::Jacobi{V}, f::AbstractQuasiVector) where V
180+
function transform_ldiv(P::Jacobi{V}, f::AbstractQuasiArray) where V
181181
T = ChebyshevT{V}()
182-
[cheb2jac(paddeddata(T \ f), P.a, P.b); zeros(V,∞)]
182+
pad(cheb2jac(paddeddata(T \ f), P.a, P.b), axes(P,2), tail(axes(f))...)
183183
end
184184

185185

‎src/classical/legendre.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -94,7 +94,7 @@ end
9494
ldiv(P::Legendre{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)
9595
function transform_ldiv(::Legendre{V}, f::Union{AbstractQuasiVector,AbstractQuasiMatrix}) where V
9696
T = ChebyshevT{V}()
97-
dat = T \ f
97+
dat = transform_ldiv(T, f)
9898
pad(th_cheb2leg(paddeddata(dat)), axes(dat)...)
9999
end
100100

‎src/lanczos.jl‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -200,16 +200,16 @@ Base.array_summary(io::IO, C::LanczosJacobiBand{T}, inds::Tuple{Vararg{OneToInf{
200200
QuasiArrays.ApplyQuasiArray(Q::LanczosPolynomial) = ApplyQuasiArray(*, arguments(ApplyLayout{typeof(*)}(), Q)...)
201201

202202

203-
function \(A::LanczosPolynomial{T}, B::LanczosPolynomial{V}) where {T,V}
204-
A == B && return Eye{promote_type(T,V)}(∞)
203+
@simplify function \(A::LanczosPolynomial, B::LanczosPolynomial)
204+
A == B && return Eye{promote_type(eltype(A),eltype(B))}(∞)
205205
inv(LanczosConversion(A.data)) * (A.P \ B.P) * LanczosConversion(B.data)
206206
end
207-
\(A::OrthogonalPolynomial, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)
208-
\(A::Normalized, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)
209-
\(Q::LanczosPolynomial, A::OrthogonalPolynomial) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
210-
\(Q::LanczosPolynomial, A::Normalized) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
211-
\(Q::LanczosPolynomial, A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
212-
\(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)
207+
@simplify \(A::OrthogonalPolynomial, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)
208+
@simplify \(A::Normalized, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)
209+
@simplify \(Q::LanczosPolynomial, A::OrthogonalPolynomial) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
210+
@simplify \(Q::LanczosPolynomial, A::Normalized) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
211+
@simplify \(Q::LanczosPolynomial, A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
212+
@simplify \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)
213213

214214
ArrayLayouts.mul(Q::LanczosPolynomial, C::AbstractArray) = ApplyQuasiArray(*, Q, C)
215215
function ldiv(Qn::SubQuasiArray{<:Any,2,<:LanczosPolynomial,<:Tuple{<:Inclusion,<:Any}}, C::AbstractQuasiArray)

‎src/normalized.jl‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -116,17 +116,20 @@ copy(M::Mul{<:AdjointBasisLayout{<:NormalizedOPLayout},Blay}) where Blay<:Abstra
116116
# table stable identity if A.P == B.P
117117
@inline _normalized_ldiv(An, C, Bn) = An \ (C * Bn)
118118
@inline _normalized_ldiv(An, C::Eye{T}, Bn) where T = FillArrays.SquareEye{promote_type(eltype(An),T,eltype(Bn))}(ℵ₀)
119+
120+
simplifiable(::Ldiv{<:NormalizedOPLayout,<:NormalizedOPLayout}) = Val(true)
119121
@inline copy(L::Ldiv{<:NormalizedOPLayout,<:NormalizedOPLayout}) = _normalized_ldiv(Diagonal(L.A.scaling), L.A.P \ L.B.P, Diagonal(L.B.scaling))
120122
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,<:AbstractNormalizedOPLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A, L.B))
121123
@inline copy(L::Ldiv{Lay,<:AbstractNormalizedOPLayout}) where Lay = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
122124
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,Lay}) where Lay = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
123-
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,Lay,<:Any,<:AbstractQuasiVector}) where Lay = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
124125
@inline copy(L::Ldiv{Lay,<:AbstractNormalizedOPLayout}) where Lay<:AbstractBasisLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
125126
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
126127
@inline copy(L::Ldiv{Lay,<:AbstractNormalizedOPLayout}) where Lay<:AbstractLazyLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
127128
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,Lay}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
128-
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,Lay,<:Any,<:AbstractQuasiVector}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
129+
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,<:ExpansionLayout}) = copy(Ldiv{BasisLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
129130
@inline copy(L::Ldiv{ApplyLayout{typeof(*)},<:AbstractNormalizedOPLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A, L.B))
131+
@inline copy(L::Ldiv{MappedOPLayout,<:AbstractNormalizedOPLayout}) = copy(Ldiv{MappedOPLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
132+
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,ApplyLayout{typeof(hcat)}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(hcat)}}(L.A, L.B))
130133
for Lay in (:(ApplyLayout{typeof(*)}),:(BroadcastLayout{typeof(+)}),:(BroadcastLayout{typeof(-)}))
131134
@eval begin
132135
@inline copy(L::Ldiv{<:AbstractNormalizedOPLayout,$Lay}) = copy(Ldiv{ApplyLayout{typeof(*)},$Lay}(L.A, L.B))
@@ -146,7 +149,9 @@ function _norm_expand_ldiv(A, w_B)
146149
B̃,D = arguments(ApplyLayout{typeof(*)}(), B)
147150
(A \ (w .* B̃)) * D
148151
end
152+
simplifiable(::Ldiv{<:AbstractNormalizedOPLayout,<:WeightedBasisLayout{<:AbstractNormalizedOPLayout}}) = Val(true)
149153
copy(L::Ldiv{<:AbstractNormalizedOPLayout,<:WeightedBasisLayout{<:AbstractNormalizedOPLayout}}) = _norm_expand_ldiv(L.A, L.B)
154+
simplifiable(::Ldiv{OPLayout,<:WeightedBasisLayout{<:AbstractNormalizedOPLayout}}) = Val(true)
150155
copy(L::Ldiv{OPLayout,<:WeightedBasisLayout{<:AbstractNormalizedOPLayout}}) = _norm_expand_ldiv(L.A, L.B)
151156

152157
###

‎test/test_chebyshev.jl‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -179,6 +179,13 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
179179
@test_throws ErrorException (f + g)[0.1]
180180
@test_throws ErrorException (g + f)[0.1]
181181
end
182+
183+
@testset "broadcast +" begin
184+
x = Inclusion(0..1)
185+
T = Chebyshev()[2x .- 1,:]
186+
@test T \ (exp.(x) .+ cos.(x)) ≈ transform(T, x -> exp(x)+cos(x))
187+
@test T \ (exp.(x) .- cos.(x)) ≈ transform(T, x -> exp(x)-cos(x))
188+
end
182189
end
183190

184191
@testset "weighted" begin
@@ -570,6 +577,8 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
570577
f = T * [1:3; zeros(∞)]
571578
g = chebyshevt(0..1) * [1:3; zeros(∞)]
572579
@test_broken (f + g)[0.1] ≈ f[0.1] + g[0.1] # ContinuumArrays needs to check maps are equal
580+
581+
@test ClassicalOrthogonalPolynomials.singularities(T) === LegendreWeight()[QuadraticMap()]
573582
end
574583

575584
@testset "block structure" begin

‎test/test_jacobi.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -346,10 +346,10 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
346346
@test Jacobi(0.0,0.0) \ Legendre() == Eye(∞)
347347
@test ((Ultraspherical(3/2) \ Jacobi(1,1))*(Jacobi(1,1) \ Ultraspherical(3/2)))[1:10,1:10] ≈ Eye(10)
348348
f = Jacobi(0.0,0.0)*[[1,2,3]; zeros(∞)]
349-
g = (Legendre() \ f) - f.args[2]
350-
@test_skip norm(g) ≤ 1E-15
351-
@test_broken (Legendre() \ f) == f.args[2]
352-
@test (Legendre() \ f)[1:10] ≈ f.args[2][1:10]
349+
g = (Legendre() \ f) - coefficients(f)
350+
@test_skip norm(g) ≤ 1E-15
351+
@test (Legendre() \ f) == coefficients(f)
352+
@test (Legendre() \ f)[1:10] ≈ coefficients(f)[1:10]
353353
f = Jacobi(1.0,1.0)*[[1,2,3]; zeros(∞)]
354354
g = Ultraspherical(3/2)*(Ultraspherical(3/2)\f)
355355
@test f[0.1] ≈ g[0.1]

0 commit comments

Comments
 (0)