ContinuumArrays v0.15

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.11.1"
+version = "0.11.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,7 +29,7 @@ ArrayLayouts = "1.0.1"
 BandedMatrices = "0.17.17"
 BlockArrays = "0.16.9"
 BlockBandedMatrices = "0.12"
-ContinuumArrays = "0.14"
+ContinuumArrays = "0.15"
 DomainSets = "0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.5"

src/ClassicalOrthogonalPolynomials.jl

@@ -30,14 +30,14 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
                     _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector,
-                    AbstractQuasiFill, _equals, QuasiArrayLayout, PolynomialLayout, diff_layout
+                    AbstractQuasiFill, equals_layout, QuasiArrayLayout, PolynomialLayout, diff_layout
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import InfiniteLinearAlgebra: chop!, chop, pad, choplength, compatible_resize!
 import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
-    grid, plotgrid, _plotgrid, _grid, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
+    grid, plotgrid, plotgrid_layout, plotvalues_layout, grid_layout, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
     AffineQuasiVector, AffineMap, AbstractWeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
-    checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, _broadcastbasis,
+    checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, broadcastbasis_layout,
     plan_transform, plan_grid_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg
@@ -87,26 +87,26 @@ end
 isorthogonalityweighted(wS) = isorthogonalityweighted(MemoryLayout(wS), wS)
 
 
-_equals(::MappedOPLayout, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
-_equals(::MappedOPLayout, ::MappedBasisLayouts, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
-_equals(::MappedBasisLayouts, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
+equals_layout(::MappedOPLayout, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
+equals_layout(::MappedOPLayout, ::MappedBasisLayouts, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
+equals_layout(::MappedBasisLayouts, ::MappedOPLayout, P, Q) = demap(P) == demap(Q) && basismap(P) == basismap(Q)
 
-_broadcastbasis(::typeof(+), ::MappedOPLayout, ::MappedOPLayout, P, Q) = _broadcastbasis(+, MappedBasisLayout(), MappedBasisLayout(), P, Q)
-_broadcastbasis(::typeof(+), ::MappedOPLayout, M::MappedBasisLayout, P, Q) = _broadcastbasis(+, MappedBasisLayout(), M, P, Q)
-_broadcastbasis(::typeof(+), L::MappedBasisLayout, ::MappedOPLayout, P, Q) = _broadcastbasis(+, L, MappedBasisLayout(), P, Q)
+broadcastbasis_layout(::typeof(+), ::MappedOPLayout, ::MappedOPLayout, P, Q) = broadcastbasis_layout(+, MappedBasisLayout(), MappedBasisLayout(), P, Q)
+broadcastbasis_layout(::typeof(+), ::MappedOPLayout, M::MappedBasisLayout, P, Q) = broadcastbasis_layout(+, MappedBasisLayout(), M, P, Q)
+broadcastbasis_layout(::typeof(+), L::MappedBasisLayout, ::MappedOPLayout, P, Q) = broadcastbasis_layout(+, L, MappedBasisLayout(), P, Q)
 sum_layout(::MappedOPLayout, A, dims) = sum_layout(MappedBasisLayout(), A, dims)
 
 # demap to avoid Golub-Welsch fallback
 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)
 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)
 
-_equals(::AbstractOPLayout, ::AbstractWeightedBasisLayout, _, _) = false # Weighted-Legendre doesn't exist
-_equals(::AbstractWeightedBasisLayout, ::AbstractOPLayout, _, _) = false # Weighted-Legendre doesn't exist
+equals_layout(::AbstractOPLayout, ::AbstractWeightedBasisLayout, _, _) = false # Weighted-Legendre doesn't exist
+equals_layout(::AbstractWeightedBasisLayout, ::AbstractOPLayout, _, _) = false # Weighted-Legendre doesn't exist
 
-_equals(::WeightedOPLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP) == unweighted(wQ)
-_equals(::WeightedOPLayout, ::WeightedBasisLayout, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
-_equals(::WeightedBasisLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
-_equals(::WeightedBasisLayout{<:AbstractOPLayout}, ::WeightedBasisLayout{<:AbstractOPLayout}, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
+equals_layout(::WeightedOPLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP) == unweighted(wQ)
+equals_layout(::WeightedOPLayout, ::WeightedBasisLayout, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
+equals_layout(::WeightedBasisLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
+equals_layout(::WeightedBasisLayout{<:AbstractOPLayout}, ::WeightedBasisLayout{<:AbstractOPLayout}, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
 
 
 copy(L::Ldiv{MappedOPLayout,Lay}) where Lay<:MappedBasisLayouts = copy(Ldiv{MappedBasisLayout,Lay}(L.A,L.B))
@@ -255,10 +255,11 @@ _tritrunc(X, n) = _tritrunc(MemoryLayout(X), X, n)
 jacobimatrix(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Inclusion,OneTo}}) =
     _tritrunc(jacobimatrix(parent(V)), maximum(parentindices(V)[2]))
 
-_grid(::AbstractOPLayout, P, n::Integer) = eigvals(symtridiagonalize(jacobimatrix(P[:,OneTo(n)])))
-_grid(::MappedOPLayout, P, n::Integer) = _grid(MappedBasisLayout(), P, n)
-_plotgrid(::AbstractOPLayout, P, n::Integer) = grid(P, min(40n, MAX_PLOT_POINTS))
-_plotgrid(::MappedOPLayout, P, n::Integer) = _plotgrid(MappedBasisLayout(), P, n)
+grid_layout(::AbstractOPLayout, P, n::Integer) = eigvals(symtridiagonalize(jacobimatrix(P[:,OneTo(n)])))
+grid_layout(::MappedOPLayout, P, n::Integer) = grid_layout(MappedBasisLayout(), P, n)
+plotgrid_layout(::AbstractOPLayout, P, n::Integer) = grid(P, min(40n, MAX_PLOT_POINTS))
+plotgrid_layout(::MappedOPLayout, P, n::Integer) = plotgrid_layout(MappedBasisLayout(), P, n)
+plotvalues_layout(::ExpansionLayout{MappedOPLayout}, f, x...) = plotvalues_layout(ExpansionLayout{MappedBasisLayout}(), f, x...)
 
 function golubwelsch(X)
     D, V = eigen(symtridiagonalize(X))  # Eigenvalue decomposition

src/normalized.jl

@@ -55,9 +55,9 @@ Base.iterate(S::QuasiQR, ::Val{:done}) = nothing
 
 axes(Q::Normalized) = axes(Q.P)
 
-_equals(::AbstractNormalizedOPLayout, ::AbstractNormalizedOPLayout, P, Q) = orthogonalityweight(P) == orthogonalityweight(Q)
-_equals(::AbstractNormalizedOPLayout, ::AbstractOPLayout, P, Q) = isnormalized(Q) && orthogonalityweight(P) == orthogonalityweight(Q)
-_equals(::AbstractOPLayout, ::AbstractNormalizedOPLayout, P, Q) = isnormalized(P) && orthogonalityweight(P) == orthogonalityweight(Q)
+equals_layout(::AbstractNormalizedOPLayout, ::AbstractNormalizedOPLayout, P, Q) = orthogonalityweight(P) == orthogonalityweight(Q)
+equals_layout(::AbstractNormalizedOPLayout, ::AbstractOPLayout, P, Q) = isnormalized(Q) && orthogonalityweight(P) == orthogonalityweight(Q)
+equals_layout(::AbstractOPLayout, ::AbstractNormalizedOPLayout, P, Q) = isnormalized(P) && orthogonalityweight(P) == orthogonalityweight(Q)
 
 
 

test/test_chebyshev.jl

@@ -152,6 +152,7 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
         @test chebyshevt(0..1) == chebyshevt(Inclusion(0..1)) == chebyshevt(T) == T
         @test (T \ [x exp.(x)])[1:20,1] ≈ (T\x)[1:20]
         @test (T \ [x exp.(x)])[1:20,2] ≈ (T\exp.(x))[1:20]
+        @test grid(T, 5) == (grid(Chebyshev(), 5) .+ 1)/2
 
         Tn = Chebyshev()[2x .- 1, [1,3,4]]
         @test (axes(Tn,1) .* Tn).args[2][1:5,:] ≈ (axes(T,1) .* T).args[2][1:5,[1,3,4]]
@@ -169,7 +170,14 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
         end
 
         @testset "broadcast" begin
+            x = Inclusion(0..1)
+            T = Chebyshev()[2x .- 1,:]
             @test (x.^2 .* T)[0.1,1:10] ≈ 0.1^2 * T[0.1,1:10]
+
+            f = T * [1:3; zeros(∞)]
+            g = LinearSpline(-1:1)[parentindices(T)[1],:]  * [1:3;]
+            @test_throws ErrorException (f + g)[0.1]
+            @test_throws ErrorException (g + f)[0.1]
         end
     end
 
@@ -228,11 +236,15 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
             x = Inclusion(0..1)
             wT̃ = wT[2x .- 1, :]
             @test MemoryLayout(wT̃) isa MappedWeightedBasisLayout
+            @test wT̃ == wT̃
+            @test wT̃ ≠ T[2x .- 1, :]
+            @test T[2x .- 1, :] ≠ wT̃
             v = wT̃ * (wT̃ \ @.(exp(x)/(sqrt(x)*sqrt(1-x))))
             @test v[0.1] ≈ let x = 0.1; exp(x)/(sqrt(x)*sqrt(1-x)) end
 
             WT̃ = w[2x .- 1] .* T[2x .- 1, :]
             @test MemoryLayout(WT̃) isa WeightedBasisLayout{MappedOPLayout}
+            @test WT̃ ≠ T[2x .- 1, :]
             v = WT̃ * (WT̃ \ @.(exp(x)/(sqrt(x)*sqrt(1-x))))
             @test v[0.1] ≈ let x = 0.1; exp(x)/(sqrt(x)*sqrt(1-x)) end
 
@@ -254,6 +266,18 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
                 c = T \ sqrt.(x.^2 .- 1)
                 @test [T[begin,:]'; L] \ [sqrt(2^2-1); zeros(∞)] ≈ c
             end
+
+            @testset "algebra" begin
+                T = chebyshevt(2..3)
+                U = chebyshevu(2..3)
+                wT = Weighted(ChebyshevT())[parentindices(T)...]
+                f = T * [1:3; zeros(∞)]
+                g = U * [1:3; zeros(∞)]
+                h = wT * [1:3; zeros(∞)]
+
+                @test (f + g)[2.1] ≈ f[2.1]+g[2.1]
+                @test (f - g)[2.1] ≈ f[2.1]-g[2.1]
+            end
         end
     end
 
@@ -542,6 +566,10 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
     un2 = Tn2 \ (2 * x .^2 .- 1)
     u = T \ (2 * x .^2 .- 1)
     @test un ≈ un2 ≈ u[1:10]
+
+    f = T * [1:3; zeros(∞)]
+    g = chebyshevt(0..1) * [1:3; zeros(∞)]
+    @test_broken (f + g)[0.1] ≈ f[0.1] + g[0.1] # ContinuumArrays needs to check maps are equal
 end
 
 @testset "block structure" begin