simplify jacobimatrix setup

Author: dlfivefifty

src/ContinuumArrays.jl

@@ -1,8 +1,9 @@
 module ContinuumArrays
 using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, InfiniteArrays, DomainSets, InfiniteLinearAlgebra, QuasiArrays
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
-                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy
-import Base.Broadcast: materialize
+                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
+                first, last
+import Base.Broadcast: materialize, BroadcastStyle
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, adjointlayout, LdivApplyStyle
 import LinearAlgebra: pinv
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
@@ -11,7 +12,7 @@ import FillArrays: AbstractFill, getindex_value
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
-                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout
+                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle
 
 export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, JacobiWeight, Jacobi, Legendre, Chebyshev, Ultraspherical,
             fullmaterialize
@@ -21,6 +22,9 @@ export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, JacobiWeig
 ####
 struct AlephInfinity{N} <: Integer end
 
+==(::AlephInfinity, ::Int) = false
+==(::Int, ::AlephInfinity) = false
+
 const ℵ₁ = AlephInfinity{1}()
 
 
@@ -30,6 +34,9 @@ const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}
 cardinality(::AbstractInterval) = ℵ₁
 *(ℵ::AlephInfinity) = ℵ
 
+first(S::Inclusion{<:Any,<:AbstractInterval}) = leftendpoint(S.domain)
+last(S::Inclusion{<:Any,<:AbstractInterval}) = rightendpoint(S.domain)
+
 
 checkindex(::Type{Bool}, inds::AbstractInterval, i::Number) = (leftendpoint(inds) <= i) & (i <= rightendpoint(inds))
 checkindex(::Type{Bool}, inds::AbstractInterval, i::Inclusion) = i.domain ⊆ inds
@@ -52,7 +59,7 @@ function materialize(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:Abstra
     A*P
 end
 
-
+BroadcastStyle(::Type{<:Inclusion{<:Any,<:AbstractInterval}}) = LazyQuasiArrayStyle{1}()
 
 include("operators.jl")
 include("bases/bases.jl")

src/bases/jacobi.jl

@@ -45,8 +45,17 @@ axes(::AbstractJacobi) = (Inclusion(ChebyshevInterval()), OneTo(∞))
 # Jacobi Matrix
 ########
 
-@simplify \(A::Legendre, *(B::Identity, C::Legendre)) =
-    _BandedMatrix(Vcat(((0:∞)./(1:2:∞))', Zeros(1,∞), ((1:∞)./(1:2:∞))'), ∞, 1,1)
+jacobimatrix(::Legendre) = _BandedMatrix(Vcat(((0:∞)./(1:2:∞))', Zeros(1,∞), ((1:∞)./(1:2:∞))'), ∞, 1,1)
+
+function jacobimatrix(J::Jacobi) 
+    b,a = J.b,J.a
+    n = 0:∞
+    B = @. 2*(n+1)*(n+a+b+1) / ((2n+a+b+1)*(2n+a+b+2))
+    A = (a^2-b^2) ./ ((2n.+a.+b).*(2n.+a.+b.+2))
+    C = @. 2*(n+a)*(n+b) / ((2n+a+b)*(2n+a+b+1))
+
+    _BandedMatrix(Vcat(C',A',B'), ∞, 1,1)
+end
 
 
 @simplify *(X::Identity, P::Legendre) = ApplyQuasiMatrix(*, P, P\(X*P))

src/bases/orthogonalpolynomials.jl

@@ -1,7 +1,14 @@
 abstract type OrthogonalPolynomial{T} <: Basis{T} end
 
-@inline jacobioperator(P::OrthogonalPolynomial) =
-    materialize(applied(\, P, applied(*, Diagonal(axes(P,1)), P)))
+
+
+@simplify function \(A::OrthogonalPolynomial, *(B::Identity, C::OrthogonalPolynomial))
+    A === C || error("Override  $(typeof(A)) \\ (x * $(typeof(C)))")
+    jacobimatrix(A)
+end
+
+
+@simplify *(X::Identity, P::OrthogonalPolynomial) = ApplyQuasiMatrix(*, P, apply(\, P, applied(*, X, P)))
 
 
 function forwardrecurrence!(v::AbstractVector{T}, b::AbstractVector, a::AbstractVector, c::AbstractVector, x) where T
@@ -60,13 +67,13 @@ _vec(a::Adjoint{<:Any,<:AbstractVector}) = a'
 bands(J) = _vec.(J.data.args)
 
 function getindex(P::OrthogonalPolynomial{T}, x::Number, n::OneTo) where T
-    J = jacobioperator(P)
+    J = jacobimatrix(P)
     b,a,c = bands(J)
     forwardrecurrence!(similar(n,T),b,a,c,x)
 end
 
 function getindex(P::OrthogonalPolynomial{T}, x::AbstractVector, n::OneTo) where T
-    J = jacobioperator(P)
+    J = jacobimatrix(P)
     b,a,c = bands(J)
     V = Matrix{T}(undef,length(x),length(n))
     for k = eachindex(x)

src/bases/ultraspherical.jl

@@ -32,14 +32,9 @@ Ultraspherical(λ::Λ) where Λ = Ultraspherical{Float64,Λ}(λ)
 # Jacobi Matrix
 ########
 
-@simplify function \(C::Chebyshev, *(X::Identity, C2::Chebyshev))
-    T = promote_type(eltype(C), eltype(X), eltype(C2))
-    _BandedMatrix(Vcat(Fill(one(T)/2,1,∞), Zeros(1,∞), Hcat(one(T), Fill(one(T)/2,1,∞))), ∞, 1, 1)
-end
-
-@simplify *(X::Identity, P::Chebyshev) = ApplyQuasiMatrix(*, P, apply(\, P, applied(*, X, P)))
+jacobimatrix(C::Chebyshev{T}) where T = _BandedMatrix(Vcat(Fill(one(T)/2,1,∞), Zeros(1,∞), Hcat(one(T), Fill(one(T)/2,1,∞))), ∞, 1, 1)
 
-@simplify function \(P::Ultraspherical, *(X::Identity, U::Ultraspherical))
+function jacobimatrix(P::Ultraspherical)
     T = promote_type(eltype(P), eltype(X), eltype(U))
     λ = P.λ
     λ == U.λ || throw(ArgumentError())
@@ -48,8 +43,6 @@ end
                         ((one(T):∞) ./ (2 .*((zero(T):∞) .+ λ)))'), ∞, 1, 1)
 end
 
-@simplify *(X::Identity, P::Ultraspherical) = ApplyQuasiMatrix(*, P, apply(\, P, applied(*, X, P)))
-
 
 ##########
 # Derivatives

test/runtests.jl

@@ -1,15 +1,20 @@
 using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, Test, InfiniteArrays, ForwardDiff
 import ContinuumArrays: ℵ₁, materialize, Chebyshev, Ultraspherical, jacobioperator, SimplifyStyle
-import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LmaterializeApplyStyle
+import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle, LmaterializeApplyStyle
 import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport
 import ForwardDiff: Dual
 
 
 @testset "Inclusion" begin
-    @test_throws InexactError Inclusion(-1..1)[0.1]
-    @test Inclusion(-1..1)[0.0] === 0
-    X = QuasiDiagonal(Inclusion(-1.0..1))
+    x = Inclusion(-1..1)
+    @test_throws InexactError x[0.1]
+    @test x[0.0] === 0
+    x = Inclusion(-1.0..1)
+    X = QuasiDiagonal(x)
     @test X[-1:0.1:1,-1:0.1:1] == Diagonal(-1:0.1:1)
+    @test Base.BroadcastStyle(typeof(x)) == LazyQuasiArrayStyle{1}()
+    @test x .* x isa BroadcastQuasiArray
+    @test (x.*x)[0.1] == 0.1^2
 end
 
 @testset "DiracDelta" begin
@@ -200,199 +205,4 @@ end
     @test u[0.1] ≈ 0.00012678835289369413
 end
 
-@testset "Ultraspherical" begin
-    T = Chebyshev()
-    U = Ultraspherical(1)
-    C = Ultraspherical(2)
-    D = Derivative(axes(T,1))
-
-    @test T\T === pinv(T)*T === Eye(∞)
-    @test U\U === pinv(U)*U === Eye(∞)
-    @test C\C === pinv(C)*C === Eye(∞)
-
-    @test ApplyStyle(\,typeof(U),typeof(applied(*,D,T))) == SimplifyStyle()
-    @test materialize(@~ U\(D*T)) isa BandedMatrix
-    D₀ = U\(D*T)
-    @test_broken D₀ isa BandedMatrix
-    @test D₀[1:10,1:10] isa BandedMatrix{Float64}
-    @test D₀[1:10,1:10] == diagm(1 => 1:9)
-    @test colsupport(D₀,1) == 1:0
-
-    D₁ = C\(D*U)
-    @test D₁ isa BandedMatrix
-    @test apply(*,D₁,D₀.args...)[1:10,1:10] == diagm(2 => 4:2:18)
-    @test (D₁*D₀)[1:10,1:10] == diagm(2 => 4:2:18)
-
-    S₀ = (U\T)[1:10,1:10]
-    @test S₀ isa BandedMatrix{Float64}
-    @test S₀ == diagm(0 => [1.0; fill(0.5,9)], 2=> fill(-0.5,8))
-
-    S₁ = (C\U)[1:10,1:10]
-    @test S₁ isa BandedMatrix{Float64}
-    @test S₁ == diagm(0 => 1 ./ (1:10), 2=> -(1 ./ (3:10)))
-end
-
-@testset "Jacobi" begin
-    S = Jacobi(true,true)
-    W = Diagonal(JacobiWeight(true,true))
-    D = Derivative(axes(W,1))
-    P = Legendre()
-
-    @test pinv(pinv(S)) === S
-    @test P\P === pinv(P)*P === Eye(∞)
-
-    Bi = pinv(Jacobi(2,2))
-    @test Bi isa QuasiArrays.PInvQuasiMatrix
-
-    A = apply(\, Jacobi(2,2), applied(*, D, S))
-    @test A isa BandedMatrix
-    A = Jacobi(2,2) \ (D*S)
-    @test typeof(A) == typeof(pinv(Jacobi(2,2))*(D*S))
-
-    @test A isa MulMatrix
-    @test isbanded(A)
-    @test bandwidths(A) == (-1,1)
-    @test size(A) == (∞,∞)
-    @test A[1:10,1:10] == diagm(1 => 1:0.5:5)
-
-    @test_broken @inferred(D*S)
-    M = D*S
-    @test M isa MulQuasiMatrix
-    @test M.args[1] == Jacobi(2,2)
-    @test M.args[2][1:10,1:10] == A[1:10,1:10]
-
-    L = Diagonal(JacobiWeight(true,false))
-    @test apply(\, Jacobi(false,true), applied(*,L,S)) isa BandedMatrix
-    @test_broken @inferred(Jacobi(false,true)\(L*S))
-    A = Jacobi(false,true)\(L*S)
-    @test A isa BandedMatrix
-    @test size(A) == (∞,∞)
-
-    L = Diagonal(JacobiWeight(false,true))
-    @test_broken @inferred(Jacobi(true,false)\(L*S))
-    A = Jacobi(true,false)\(L*S)
-    @test A isa BandedMatrix
-    @test size(A) == (∞,∞)
-
-    A,B = (P'P),P\(W*S)
-
-    M = Mul(A,B)
-    @test M[1,1] == 4/3
-
-    M = ApplyMatrix{Float64}(*,A,B)
-    M̃ = M[1:10,1:10]
-    @test M̃ isa BandedMatrix
-    @test bandwidths(M̃) == (2,0)
-
-    @test A*B isa MulArray
-
-    A,B,C = (P\(W*S))',(P'P),P\(W*S)
-    M = ApplyArray(*,A,B,C)
-    @test bandwidths(M) == (2,2)
-    @test M[1,1] ≈  1+1/15
-    @test typeof(M) == typeof(A*B*C)
-    M = A*B*C
-    @test bandwidths(M) == (2,2)
-    @test M[1,1] ≈  1+1/15
-end
-
-@testset "P-FEM" begin
-    S = Jacobi(true,true)
-    W = Diagonal(JacobiWeight(true,true))
-    D = Derivative(axes(W,1))
-    P = Legendre()
-    
-    M = P\(D*W*S)
-    @test M isa ApplyArray
-    @test M[1:10,1:10] == diagm(-1 => -2.0:-2:-18.0)
-
-    N = 10
-    A = D*W*S[:,1:N]
-    @test A.args[1] == P    
-    @test P\((D*W)*S[:,1:N]) isa AbstractMatrix
-    @test P\(D*W*S[:,1:N]) isa AbstractMatrix
-
-    L = D*W*S
-    Δ = L'L
-    @test Δ isa MulMatrix
-    @test Δ[1:3,1:3] isa BandedMatrix
-    @test bandwidths(Δ) == (0,0)
-
-    L = D*W*S[:,1:N]
-
-    A  = apply(*, (L').args..., L.args...)
-    @test A isa MulQuasiMatrix
-
-    A  = *((L').args..., L.args...)
-    @test A isa MulQuasiMatrix
-
-    @test apply(*,L',L) isa QuasiArrays.ApplyQuasiArray
-
-    Δ = L'L
-    @test Δ isa MulMatrix
-    @test bandwidths(Δ) == (0,0)
-end
-
-@testset "Chebyshev evaluation" begin
-    P = Chebyshev()
-    @test @inferred(P[0.1,Base.OneTo(0)]) == Float64[]
-    @test @inferred(P[0.1,Base.OneTo(1)]) == [1.0]
-    @test @inferred(P[0.1,Base.OneTo(2)]) == [1.0,0.1]
-    for N = 1:10
-        @test @inferred(P[0.1,Base.OneTo(N)]) ≈ @inferred(P[0.1,1:N]) ≈ [cos(n*acos(0.1)) for n = 0:N-1]
-        @test @inferred(P[0.1,N]) ≈ cos((N-1)*acos(0.1))
-    end
-    @test P[0.1,[2,5,10]] ≈ [0.1,cos(4acos(0.1)),cos(9acos(0.1))]
-
-    P = Ultraspherical(1)
-    @test @inferred(P[0.1,Base.OneTo(0)]) == Float64[]
-    @test @inferred(P[0.1,Base.OneTo(1)]) == [1.0]
-    @test @inferred(P[0.1,Base.OneTo(2)]) == [1.0,0.2]
-    for N = 1:10
-        @test @inferred(P[0.1,Base.OneTo(N)]) ≈ @inferred(P[0.1,1:N]) ≈ [sin((n+1)*acos(0.1))/sin(acos(0.1)) for n = 0:N-1]
-        @test @inferred(P[0.1,N]) ≈ sin(N*acos(0.1))/sin(acos(0.1))
-    end
-    @test P[0.1,[2,5,10]] ≈ [0.2,sin(5acos(0.1))/sin(acos(0.1)),sin(10acos(0.1))/sin(acos(0.1))]
-
-    P = Ultraspherical(2)
-    @test @inferred(P[0.1,Base.OneTo(0)]) == Float64[]
-    @test @inferred(P[0.1,Base.OneTo(1)]) == [1.0]
-    @test @inferred(P[0.1,Base.OneTo(2)]) == [1.0,0.4]
-    @test @inferred(P[0.1,Base.OneTo(3)]) == [1.0,0.4,-1.88]
-end
-
-@testset "Collocation" begin
-    P = Chebyshev()
-    D = Derivative(axes(P,1))
-    n = 300
-    x = cos.((0:n-2) .* π ./ (n-2))
-    cfs = [P[-1,1:n]'; (D*P)[x,1:n] - P[x,1:n]] \ [exp(-1); zeros(n-1)]
-    u = P[:,1:n]*cfs
-    @test u[0.1] ≈ exp(0.1)
-
-    P = Chebyshev()
-    D = Derivative(axes(P,1))
-    D2 = D*(D*P) # could be D^2*P in the future
-    n = 300
-    x = cos.((1:n-2) .* π ./ (n-1)) # interior Chebyshev points 
-    C = [P[-1,1:n]';
-         D2[x,1:n] + P[x,1:n];
-         P[1,1:n]']
-    cfs = C \ [1; zeros(n-2); 2] # Chebyshev coefficients
-    u = P[:,1:n]*cfs  # interpret in basis
-    @test u[0.1] ≈ (3cos(0.1)sec(1) + csc(1)sin(0.1))/2
-end
-
-@testset "Auto-diff" begin
-    U = Ultraspherical(1)
-    C = Ultraspherical(2)
-
-    f = x -> Chebyshev{eltype(x)}()[x,5]
-    @test ForwardDiff.derivative(f,0.1) ≈ 4*U[0.1,4]
-    f = x -> Chebyshev{eltype(x)}()[x,5][1]
-    @test ForwardDiff.gradient(f,[0.1]) ≈ [4*U[0.1,4]]
-    @test ForwardDiff.hessian(f,[0.1]) ≈ [8*C[0.1,3]]
-
-    f = x -> Chebyshev{eltype(x)}()[x,1:5]
-    @test ForwardDiff.derivative(f,0.1) ≈ [0;(1:4).*U[0.1,1:4]]
-end
+include("test_orthogonalpolynomials.jl")