ContinuumArray 0.19 (#227)

* ContinuumArray 0.19

* fix diff_layout

* Update ultraspherical.jl

* JacobiWeight template

* Update Project.toml

* Update normalized.jl

* add tests

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.14.4"
+version = "0.15"
 
 
 [deps]
@@ -38,23 +38,23 @@ ArrayLayouts = "1.3.1"
 BandedMatrices = "1"
 BlockArrays = "1"
 BlockBandedMatrices = "0.13"
-ContinuumArrays = "0.18.3"
+ContinuumArrays = "0.19"
 DomainSets = "0.6, 0.7"
 DynamicPolynomials = "0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "1"
 FastTransforms = "0.16.6, 0.17"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
-InfiniteArrays = " 0.14, 0.15"
-InfiniteLinearAlgebra = "0.8, 0.9"
+InfiniteArrays = " 0.15"
+InfiniteLinearAlgebra = "0.9"
 IntervalSets = "0.7"
 LazyArrays = "2.2"
-LazyBandedMatrices = "0.10, 0.11"
+LazyBandedMatrices = "0.11"
 MutableArithmetics = "1"
-QuasiArrays = "0.11"
+QuasiArrays = "0.12"
 RecurrenceRelationshipArrays = "0.1.2"
-RecurrenceRelationships = "0.1.1, 0.2"
+RecurrenceRelationships = "0.2"
 SpecialFunctions = "1.0, 2"
 julia = "1.10"
 

src/ClassicalOrthogonalPolynomials.jl

@@ -34,7 +34,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import InfiniteLinearAlgebra: chop!, chop, pad, choplength, compatible_resize!, partialcholesky!
-import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
+import ContinuumArrays: Basis, Weight, basis_axes, @simplify, AbstractAffineQuasiVector, ProjectionFactorization,
     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_layout,
@@ -90,6 +90,8 @@ represents an OP multiplied by its orthogonality weight.
 """
 struct WeightedOPLayout{Lay<:AbstractOPLayout} <: AbstractWeightedBasisLayout end
 
+grid_layout(::WeightedOPLayout, P, n) = grid(unweighted(P), n)
+
 isorthogonalityweighted(::WeightedOPLayout, _) = true
 function isorthogonalityweighted(::AbstractWeightedBasisLayout, wS)
     w,S = arguments(wS)
@@ -246,7 +248,7 @@ grammatrix_layout(::WeightedOPLayout{MappedOPLayout}, P) = grammatrix_layout(Map
 
 OrthogonalPolynomial(w::Weight) =error("Override for $(typeof(w))")
 
-@simplify *(B::Identity, C::OrthogonalPolynomial) = ApplyQuasiMatrix(*, C, jacobimatrix(C))
+@simplify *(B::QuasiDiagonal{<:Any,<:Inclusion}, C::OrthogonalPolynomial) = ApplyQuasiMatrix(*, C, jacobimatrix(C))
 
 function layout_broadcasted(::Tuple{PolynomialLayout,AbstractOPLayout}, ::typeof(*), x::Inclusion, C)
     x == axes(C,1) || throw(DimensionMismatch())

src/classical/jacobi.jl

@@ -38,13 +38,16 @@ julia> axes(J)
 (Inclusion(-1.0 .. 1.0 (Chebyshev)),)
 ```
 """
-struct JacobiWeight{T} <: AbstractJacobiWeight{T}
-    a::T
-    b::T
-    JacobiWeight{T}(a, b) where T = new{T}(convert(T,a), convert(T,b))
+struct JacobiWeight{T,V} <: AbstractJacobiWeight{T}
+    a::V
+    b::V
+    JacobiWeight{T,V}(a, b) where {T,V} = new{T,V}(convert(V,a), convert(V,b))
 end
 
-JacobiWeight(a::V, b::T) where {T,V} = JacobiWeight{promote_type(T,V)}(a,b)
+JacobiWeight{T}(a::V, b::V) where {T,V} = JacobiWeight{T,V}(a, b)
+JacobiWeight{T}(a, b) where T = JacobiWeight{T}(promote(a,b)...)
+JacobiWeight(a::V, b::T) where {T,V} = JacobiWeight{float(promote_type(T,V))}(a, b)
+
 
 """
     jacobiweight(a,b, d::AbstractInterval)
@@ -71,9 +74,9 @@ AbstractQuasiVector{T}(w::JacobiWeight) where T = JacobiWeight{T}(w.a, w.b)
 
 ==(A::JacobiWeight, B::JacobiWeight) = A.b == B.b && A.a == B.a
 
-function getindex(w::JacobiWeight, x::Number)
+function getindex(w::JacobiWeight{T}, x::Number) where T
     x ∈ axes(w,1) || throw(BoundsError())
-    (1-x)^w.a * (1+x)^w.b
+    convert(T, (1-x)^w.a * (1+x)^w.b)
 end
 
 show(io::IO, P::JacobiWeight) = summary(io, P)
@@ -479,6 +482,8 @@ end
 broadcastbasis(::typeof(+), w_A::Weighted{<:Any,<:Jacobi}, w_B::Weighted{<:Any,<:Jacobi}) = broadcastbasis(+, convert(WeightedBasis,w_A), convert(WeightedBasis,w_B))
 broadcastbasis(::typeof(+), w_A::Weighted{<:Any,<:Jacobi}, w_B::WeightedJacobi) = broadcastbasis(+, convert(WeightedBasis,w_A), w_B)
 broadcastbasis(::typeof(+), w_A::WeightedJacobi, w_B::Weighted{<:Any,<:Jacobi}) = broadcastbasis(+, w_A, convert(WeightedBasis,w_B))
+broadcastbasis(::typeof(+), A::Jacobi, B::Weighted{<:Any,<:Jacobi{<:Any,<:Integer}}) = A
+broadcastbasis(::typeof(+), A::Weighted{<:Any,<:Jacobi{<:Any,<:Integer}}, B::Jacobi) = B
 
 function \(w_A::WeightedJacobi, w_B::WeightedJacobi)
     wA,A = w_A.args

src/classical/ultraspherical.jl

@@ -139,7 +139,9 @@ end
 
 # higher order 
 
-function diff(::ChebyshevT{T}, m::Integer; dims=1) where T
+function diff(S::ChebyshevT{T}, m::Integer; dims=1) where T
+    iszero(m) && return S
+    isone(m) && return diff(S)
     μ = pochhammer(one(T),m-1)*convert(T,2)^(m-1)
     D = _BandedMatrix((μ * (0:∞))', ℵ₀, -m, m)
     ApplyQuasiMatrix(*, Ultraspherical{T}(m), D)

src/lanczos.jl

@@ -233,4 +233,4 @@ broadcastbasis(::typeof(+), P::Union{Normalized,LanczosPolynomial}, Q::Union{Nor
 broadcastbasis(::typeof(+), P::Union{Normalized,LanczosPolynomial}, Q) = broadcastbasis(+, P.P, Q)
 broadcastbasis(::typeof(+), P, Q::Union{Normalized,LanczosPolynomial}) = broadcastbasis(+, P, Q.P)
 
-diff_layout(::LanczosLayout, P::AbstractQuasiMatrix, dims...) = diff_layout(ApplyLayout{typeof(*)}(), P, dims...)
+diff_layout(::LanczosLayout, P::AbstractQuasiMatrix, order...; dims...) = diff_layout(ApplyLayout{typeof(*)}(), P, order...; dims...)

src/normalized.jl

@@ -299,6 +299,7 @@ _sum(p::SubQuasiArray{T,1,<:Weighted,<:Tuple{Inclusion,Int}}, ::Colon) where T =
 demap(W::Weighted) = Weighted(demap(W.P))
 basismap(W::Weighted) = basismap(W.P)
 const MappedOPLayouts = Union{MappedOPLayout,WeightedOPLayout{MappedOPLayout}}
-diff_layout(::MappedOPLayouts, A, dims...) = diff_layout(MappedBasisLayout(), A, dims...)
+diff_layout(::MappedOPLayouts, A, order::Int; dims...) = diff_layout(MappedBasisLayout(), A, order; dims...)
+diff_layout(::MappedOPLayouts, A, order...; dims...) = diff_layout(MappedBasisLayout(), A, order...; dims...)
 
-diff_layout(::AbstractNormalizedOPLayout, A, dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, dims...)
+diff_layout(::AbstractNormalizedOPLayout, A, order...; dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)

test/test_jacobi.jl

@@ -537,4 +537,19 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
     @testset "conversion not implemented" begin 
         @test_throws ArgumentError Jacobi(0,0) \ Jacobi(1.1,2.1)
     end
+
+    @testset "broadcastbasis" begin
+        a = Jacobi(1,1) * [1; zeros(∞)]
+        b = Weighted(Jacobi(1,1)) * [1; zeros(∞)]
+        c = Weighted(Jacobi(2,3)) * [1; zeros(∞)]
+
+        @test basis(a+b) == basis(b+a) == basis(a+c) == basis(c+a) == Jacobi(1,1)
+        @test a[0.1]+b[0.1] ≈ (a+b)[0.1] ≈ (b+a)[0.1]
+    end
+
+    @testset "Weighted expand" begin
+        W = Weighted(Jacobi(1,1))
+        @test expand(W, x -> (1-x^2)*exp(x))[0.1] ≈ (1-0.1^2)*exp(0.1)
+        @test grid(W, 5) == grid(W.P, 5)
+    end
 end

test/test_lanczos.jl

@@ -303,16 +303,23 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedColumns, or
         @test jacobimatrix(Q)[1,1] ≈ 1/3
         @test Q[0.5,1:3] ≈ [1, 1.369306393762913, 0.6469364618834543]
     end
-end
-
-@testset "#197" begin
-    @test _emptymaximum(1:5) == 5 
-    @test _emptymaximum(1:0) == 0
-    x = Inclusion(ChebyshevInterval())
-    f = exp.(x)
-    QQ = LanczosPolynomial(f)
-    R = LanczosConversion(QQ.data)
-    v = cache(Zeros(∞))
-    @test (R \ v)[1:500] == zeros(500)
-    @test (R * v)[1:500] == zeros(500)
+
+    @testset "#197" begin
+        @test _emptymaximum(1:5) == 5 
+        @test _emptymaximum(1:0) == 0
+        x = Inclusion(ChebyshevInterval())
+        f = exp.(x)
+        QQ = LanczosPolynomial(f)
+        R = LanczosConversion(QQ.data)
+        v = cache(Zeros(∞))
+        @test (R \ v)[1:500] == zeros(500)
+        @test (R * v)[1:500] == zeros(500)
+    end
+
+    @testset "diff" begin
+        P = Normalized(Legendre())
+        x = axes(P,1)
+        Q = LanczosPolynomial(exp.(x))
+        @test diff(Q)[0.1,3] ≈ -0.1637907411174539
+    end
 end