expand support for broadcastbasis with Jacobi and Ultraspherical OPs (#271)

* expand support for broadcastbasis with Jacobi and Ultraspherical OPs

* fix default broadcastbasis

* Update legendre.jl

* Update test_jacobi.jl

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
-version = "0.15.17"
+version = "0.15.18"
 authors = ["Sheehan Olver <solver@mac.com>"]
 
 [deps]

src/classical/chebyshev.jl

@@ -327,5 +327,6 @@ end
 # algebra
 ####
 
-broadcastbasis(::typeof(+), ::ChebyshevT, U::ChebyshevU) = U
-broadcastbasis(::typeof(+), U::ChebyshevU, ::ChebyshevT) = U
+broadcastbasis(::typeof(+), ::Chebyshev{kind,T}, ::Chebyshev{kind,V}) where {kind,T,V} = Chebyshev{kind,promote_type(T,V)}()
+broadcastbasis(::typeof(+), ::ChebyshevT{T}, ::ChebyshevU{V}) where {T,V} = ChebyshevU{promote_type(T,V)}()
+broadcastbasis(::typeof(+), ::ChebyshevU{T}, ::ChebyshevT{V}) where {T,V} = ChebyshevU{promote_type(T,V)}()

src/classical/jacobi.jl

@@ -159,6 +159,9 @@ Jacobi{T}(a::V, b::V) where {T,V} = Jacobi{T,V}(a, b)
 Jacobi{T}(a, b) where T = Jacobi{T}(promote(a,b)...)
 Jacobi(a::V, b::T) where {T,V} = Jacobi{float(promote_type(T,V))}(a, b)
 
+Jacobi(P::Jacobi) = P
+Jacobi{T}(P::Jacobi) where T = Jacobi{T}(P.a, P.b)
+
 AbstractQuasiArray{T}(w::Jacobi) where T = Jacobi{T}(w.a, w.b)
 AbstractQuasiMatrix{T}(w::Jacobi) where T = Jacobi{T}(w.a, w.b)
 
@@ -464,6 +467,8 @@ function \(w_A::WeightedJacobi, B::AbstractJacobi)
 end
 
 
+broadcastbasis(::typeof(+),  A::AbstractJacobi, B::AbstractJacobi) = Jacobi(max(Jacobi(A).a,Jacobi(B).a), max(Jacobi(A).b,Jacobi(B).b))
+
 function broadcastbasis(::typeof(+), w_A::WeightedJacobi, w_B::WeightedJacobi)
     wA,A = w_A.args
     wB,B = w_B.args

src/classical/legendre.jl

@@ -263,3 +263,6 @@ function _sum(P::Legendre{T}, dims::Int) where T
 end
 
 _sum(p::SubQuasiArray{T,1,Legendre{T},<:Tuple{Inclusion,Int}}, ::Colon) where T = parentindices(p)[2] == 1 ? convert(T, 2) : zero(T)
+
+
+broadcastbasis(::typeof(+), ::Legendre{T}, ::Legendre{V}) where {T,V} = Legendre{promote_type(T,V)}()

src/classical/ultraspherical.jl

@@ -234,12 +234,12 @@ end
 
 
 
-function \(U::Ultraspherical{<:Any,<:Integer}, C::ChebyshevT)
+function \(U::Ultraspherical, C::ChebyshevT)
     T = promote_type(eltype(U), eltype(C))
     (U\Ultraspherical{T}(1)) * (ChebyshevU{T}()\C)
 end
 
-function \(U::Ultraspherical{<:Any,<:Integer}, C::ChebyshevU)
+function \(U::Ultraspherical, C::ChebyshevU)
     T = promote_type(eltype(U), eltype(C))
     U\Ultraspherical(C)
 end
@@ -326,3 +326,10 @@ end
 \(A::Ultraspherical, w_B::WeightedUltraspherical) = (UltrasphericalWeight(one(A.λ)/2) .* A) \ w_B
 \(A::Legendre, wB::WeightedUltraspherical) = Ultraspherical(A) \ wB
 
+
+
+broadcastbasis(::typeof(+),  A::Ultraspherical{T}, B::Ultraspherical{V}) where {T,V} = Ultraspherical{promote_type(T,V)}(max(A.λ,B.λ))
+broadcastbasis(::typeof(+),  A::ChebyshevT, B::Ultraspherical) = Ultraspherical(max(zero(real(eltype(A))),B.λ)) # we assume B.λ > 0
+broadcastbasis(::typeof(+),  A::Ultraspherical, B::ChebyshevT) = Ultraspherical(max(A.λ,zero(real(eltype(B))))) # we assume A.λ > 0
+broadcastbasis(::typeof(+),  A::Legendre, B::Ultraspherical) = Ultraspherical(max(-one(real(eltype(A)))/2,B.λ))
+broadcastbasis(::typeof(+),  A::Ultraspherical, B::Legendre) = Ultraspherical(max(A.λ,-one(real(eltype(B)))/2))

test/test_chebyshev.jl

@@ -602,6 +602,10 @@ import BandedMatrices: isbanded
         @test W \ (a .* W) isa Clenshaw
         @test (a .* W)[0.1,1:5] ≈ a[0.1] * W[0.1,1:5]
     end
+
+    @testset "broadcastbasis" begin
+        @test (expand(ChebyshevT(), exp) + expand(ChebyshevU(), cos))[0.1] ≈ exp(0.1)+cos(0.1)
+    end
 end
 
 struct QuadraticMap{T} <: Map{T} end

test/test_jacobi.jl

@@ -28,7 +28,7 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
         @test Jacobi(1,2) .* Legendre() == Jacobi(1,2) .* Jacobi(0,0)
 
         J = Jacobi(1,2)
-        @test AbstractQuasiArray{Float32}(J) ≡ AbstractQuasiMatrix{Float32}(J) ≡ Jacobi{Float32}(1,2)
+        @test AbstractQuasiArray{Float32}(J) ≡ AbstractQuasiMatrix{Float32}(J) ≡ Jacobi{Float32}(1,2) ≡ Jacobi{Float32}(J)
     end
 
     @testset "basis" begin
@@ -561,4 +561,11 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
     @testset "ChebyshevInterval constructior" begin
         @test jacobi(1,2,ChebyshevInterval()) ≡ Jacobi(1,2)
     end
+
+    @testset "broadcastbasis" begin
+        @test (expand(Jacobi(1,2), exp) + expand(Jacobi(2,1),cos))[0.1] ≈ exp(0.1)+cos(0.1)
+        @test (expand(Legendre(), exp) + expand(Jacobi(2,1),cos))[0.1] ≈ exp(0.1)+cos(0.1)
+        @test (expand(ChebyshevT(), exp) + expand(Jacobi(3/2,1/2),cos))[0.1] ≈ exp(0.1)+cos(0.1)
+        @test (expand(Ultraspherical(2), exp) + expand(Jacobi(3/2,1/2),cos))[0.1] ≈ exp(0.1)+cos(0.1)
+    end
 end

test/test_ultraspherical.jl

@@ -227,4 +227,10 @@ using ClassicalOrthogonalPolynomials: grammatrix
         op = (s .+ im) .* T
         @test (C \ op)[1,1] ≈ 0.5 + im
     end
+
+    @testset "broadcastbasis" begin
+        @test (expand(Ultraspherical(1), exp) + expand(Ultraspherical(2),cos))[0.1] ≈ exp(0.1)+cos(0.1)
+        @test (expand(Legendre(), exp) + expand(Ultraspherical(3/2), cos))[0.1] ≈ (expand(Ultraspherical(3/2), exp) + expand(Legendre(), cos))[0.1] ≈ exp(0.1)+cos(0.1)
+        @test (expand(ChebyshevT(), exp) + expand(Ultraspherical(2), cos))[0.1] ≈ (expand(Ultraspherical(2), exp) + expand(ChebyshevT(), cos))[0.1] ≈ exp(0.1)+cos(0.1)
+    end
 end