Support broadcasting mixed with operators (#37)

* Support broadcasting and operators

* Update test_fourier.jl

* Update Project.toml

* Drop Julia v1.5

* update Project

* Increase coverage

Author: dlfivefifty

.github/workflows/ci.yml

@@ -10,8 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
-          - '^1.6.0-0'
+          - '1.6'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.3.7"
+version = "0.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -25,25 +25,25 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "0.6.2, 0.7"
+ArrayLayouts = "0.7"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.15"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.7.3"
+ContinuumArrays = "0.8"
 DomainSets = "0.4, 0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
-FastTransforms = "0.11, 0.12"
+FastTransforms = "0.12"
 FillArrays = "0.11"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.10.3"
 InfiniteLinearAlgebra = "0.5.3"
 IntervalSets = "0.5"
 LazyArrays = "0.21"
 LazyBandedMatrices = "0.5.5"
-QuasiArrays = "0.5"
-SpecialFunctions = "0.10, 1"
-julia = "1.5"
+QuasiArrays = "0.6"
+SpecialFunctions = "1"
+julia = "1.6"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

src/ClassicalOrthogonalPolynomials.jl

@@ -15,8 +15,8 @@ import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupp
                 AbstractCachedVector, AbstractCachedMatrix
 import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
                         subdiagonaldata, diagonaldata, supdiagonaldata
-import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal
-import LinearAlgebra: pinv, factorize, qr, adjoint, transpose
+import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, AbstractLazyBandedLayout
+import LinearAlgebra: pinv, factorize, qr, adjoint, transpose, dot
 import BandedMatrices: AbstractBandedLayout, AbstractBandedMatrix, _BandedMatrix, bandeddata
 import FillArrays: AbstractFill, getindex_value
 
@@ -46,7 +46,7 @@ export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomia
             ∞, Derivative, .., Inclusion, 
             chebyshevt, chebyshevu, legendre, jacobi,
             legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh, normalizedjacobip,
-            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight
+            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted
 
 if VERSION < v"1.6-"
     oneto(n) = Base.OneTo(n)

src/clenshaw.jl

@@ -263,7 +263,7 @@ bandwidths(M::Clenshaw) = (length(M.c)-1,length(M.c)-1)
 Base.array_summary(io::IO, C::Clenshaw{T}, inds::Tuple{Vararg{OneToInf{Int}}}) where T =
     print(io, Base.dims2string(length.(inds)), " Clenshaw{$T} with $(length(C.c)) degree polynomial")
 
-struct ClenshawLayout <: AbstractBandedLayout end
+struct ClenshawLayout <: AbstractLazyBandedLayout end
 MemoryLayout(::Type{<:Clenshaw}) = ClenshawLayout()
 sublayout(::ClenshawLayout, ::Type{<:NTuple{2,AbstractUnitRange{Int}}}) = ClenshawLayout()
 sub_materialize(::ClenshawLayout, V) = BandedMatrix(V)

test/test_chebyshev.jl

@@ -130,9 +130,12 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
 
         @testset "Weighted" begin
             WT = Weighted(ChebyshevT())
+            @test WT == copy(WT)
+            @test WT \ WT == Eye(∞)
             @test wT[0.1,1:10] ≈ WT[0.1,1:10]
             @test WT \ (exp.(x) ./ sqrt.(1 .- x.^2)) ≈ wT \ (exp.(x) ./ sqrt.(1 .- x.^2))
             @test WT[:,1:20] \ (exp.(x) ./ sqrt.(1 .- x.^2)) ≈ (WT \ (exp.(x) ./ sqrt.(1 .- x.^2)))[1:20]
+            @test WT \ (x .* WT) == T \ (x .* T)
         end
 
         @testset "mapped" begin

test/test_fourier.jl

@@ -39,7 +39,6 @@ import QuasiArrays: MulQuasiArray
         @test blockisequal(axes(F'F), (axes(F,2),axes(F,2)))
         D = Derivative(axes(F,1))
         D̃ = (D*F).args[2]
-        @test BlockArrays.blockcolsupport(D̃,Block(3)) == Block.(3:3)
         @test (F\F)*D̃ isa BlockArray
         @test (F \ (D*F))[Block.(1:3),Block.(1:3)] == [0 0 0 0 0; 0 0.0 -1 0 0; 0 1 0 0 0; 0 0 0 0 -2; 0 0 0 2 0]
 

test/test_normalized.jl

@@ -1,4 +1,4 @@
-using ClassicalOrthogonalPolynomials, FillArrays, BandedMatrices, ContinuumArrays, ArrayLayouts, LazyArrays, Base64, Test
+using ClassicalOrthogonalPolynomials, FillArrays, BandedMatrices, ContinuumArrays, ArrayLayouts, LazyArrays, Base64, LinearAlgebra, Test
 import ClassicalOrthogonalPolynomials: NormalizedBasisLayout, recurrencecoefficients, Normalized, Clenshaw, weighted
 import LazyArrays: CachedVector, PaddedLayout
 import ContinuumArrays: MappedWeightedBasisLayout
@@ -80,6 +80,13 @@ import ContinuumArrays: MappedWeightedBasisLayout
         @testset "show" begin
             @test stringmime("text/plain", Normalized(Legendre())) == "Normalized(Legendre{Float64})"
         end
+
+        @testset "qr" begin
+            P = Legendre()
+            Q,R = qr(P)
+            @test Q == Normalized(Legendre())
+            @test R[1:10,1:10] == (P\Q)[1:10,1:10]
+        end
     end
 
     @testset "Chebyshev" begin

test/test_odes.jl

@@ -170,4 +170,14 @@ import SemiseparableMatrices: VcatAlmostBandedLayout
         u = L \ [airyai(-ε^(-2/3)); airyai(ε^(2/3)); zeros(∞)]
         @test T[-0.1,:]'u ≈ airyai(-0.1*ε^(-2/3))
     end
+
+    @testset "combo operators" begin
+        T = Chebyshev()
+        C = Ultraspherical(2)
+        x = axes(T,1)
+        D = Derivative(x)
+        L = x .* D + cos.(x) .* D^2
+        M = C \ (L * T)
+        @test C[0.1,:]' * (M * (T \ exp.(x))) ≈ (0.1 + cos(0.1))*exp(0.1)
+    end
 end