Evaluating SemiclassicalJacobi is not type stable

[DanielVandH]

julia> P = SemiclassicalJacobi(2.0, 1/2, 1/2, 1/2)
SemiclassicalJacobi with weight x^0.5 * (1-x)^0.5 * (2.0-x)^0.5 on 0..1

julia> @inferred P[0.2, 1]
ERROR: return type Float64 does not match inferred return type Any
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] top-level scope
   @ REPL[14]:

Probably because

julia> @inferred recurrencecoefficients(P)
ERROR: return type Tuple{BroadcastVector{Float64, typeof(inv), Tuple{ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}, BroadcastVector{Float64, typeof(-), Tuple{BroadcastVector{Float64, typeof(/), Tuple{ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}}}, BroadcastVector{Float64, typeof(/), Tuple{ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}} does not match inferred return type Tuple{Any, Any, Any}
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] top-level scope
   @ REPL[15]:1

since

julia> @inferred jacobimatrix(P)
ERROR: return type LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}} does not match inferred return type AbstractMatrix{Float64}
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] top-level scope
   @ REPL[24]:1

I don't think it's possible to make jacobimatrix type stable in the current setup so probably the easiest approach is to slap a ::T in a view areas of wherever this all gets evaluated

[DanielVandH]

The @descend call with Cthulhu shows

julia> @descend P[0.2, 1]
getindex(P::OrthogonalPolynomial, x::Number, n::Number) @ ClassicalOrthogonalPolynomials C:\Users\djv23\.julia\packages\ClassicalOrthogonalPolynomials\68t49\src\clenshaw.jl:54
54 function getindex(P::SemiclassicalJacobi{Float64}::OrthogonalPolynomial, x::Float64::Number, n::Int64::Number)::Any
55     @boundscheck checkbounds(P::SemiclassicalJacobi{Float64}, x::Float64, n::Int64)
56     Base.unsafe_getindex(P::SemiclassicalJacobi{Float64}, x::Float64, n::Int64)::Any
57 end
Select a call to descend into or ↩ to ascend. [q]uit. [b]ookmark.
Toggles: [w]arn, [h]ide type-stable statements, [t]ype annotations, [s]yntax highlight for Source/LLVM/Native, [j]ump to source always, [v]scode: inlay types, [V]scode: diagnostics.
Show: [S]ource code, [A]ST, [T]yped code, [L]LVM IR, [N]ative code
Actions: [E]dit source code, [R]evise and redisplay
   %5 = checkbounds(::SemiclassicalJacobi{Float64},::Float64,::Int64)::Any
 • Base.unsafe_getindex(P::SemiclassicalJacobi{Float64}, x::Float64, n::Int64)
   ↩
unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::Number) @ ClassicalOrthogonalPolynomials C:\Users\djv23\.julia\packages\ClassicalOrthogonalPolynomials\68t49\src\clenshaw.jl:68
68 Base.unsafe_getindex(P::SemiclassicalJacobi{Float64}::OrthogonalPolynomial, x::Float64::Number, n::Int64::Number)::Any = initiateforwardrecurrence((n::Int64-1)::Int64, recurrencecoefficients(P::SemiclassicalJacobi{Float64})::Tuple{Any, Any, Any}..., x::Float64, _p0(P::SemiclassicalJacobi{Float64})::Core.Const(1.0))::Tuple{Any, Any}[end]
Select a call to descend into or ↩ to ascend. [q]uit. [b]ookmark.
Toggles: [w]arn, [h]ide type-stable statements, [t]ype annotations, [s]yntax highlight for Source/LLVM/Native, [j]ump to source always, [v]scode: inlay types, [V]scode: diagnostics.
Show: [S]ource code, [A]ST, [T]yped code, [L]LVM IR, [N]ative code
Actions: [E]dit source code, [R]evise and redisplay
   n::Int64-1
   recurrencecoefficients(P::SemiclassicalJacobi{Float64})
   _p0(P::SemiclassicalJacobi{Float64})
 • initiateforwardrecurrence((n::Int64-1)::Int64, recurrencecoefficients(P::SemiclassicalJacobi{Float64})::Tuple{Any, Any, Any}..., x::Float64, _p0(P::SemiclassicalJacobi{Float64})::Core.Const(1.0))
   initiateforwardrecurrence((n::Int64-1)::Int64, recurrencecoefficients(P::SemiclassicalJacobi{Float64})::Tuple{Any, Any, Any}..., x::Float64, _p0(P::SemiclassicalJacobi{Float64})::Core.Const(1.0))::Tuple{Any, Any}[end]
   %9 = < constprop > getindex(::Tuple{Any, Any},::Core.Const(2))::Any
   ↩
initiateforwardrecurrence(N, A, B, C, x, μ) @ RecurrenceRelationshipArrays C:\Users\djv23\.julia\packages\RecurrenceRelationshipArrays\SCsob\src\recurrence.jl:26
26 function initiateforwardrecurrence(N::Int64, A::Any, B::Any, C::Any, x::Float64, μ::Float64)::Tuple{Any, Any}
27     T::Any = promote_type(eltype(A::Any)::Any, eltype(B::Any)::Any, eltype(C::Any)::Any, typeof(x::Float64)::Type{Float64})::Any
28     p0::Any = convert(T::Any, μ::Float64)::Any
29     (N::Int64 == 0)::Bool && return zero(T::Any)::Any, p0::Any::Tuple{Any, Any}
30     p1::Any = convert(T::Any, (muladd(A::Any[1]::Any,x::Float64,B::Any[1]::Any)::Any*p0::Any)::Any)::Any
31     @inbounds for n::Int64 = (2:N::Int64)::Int64::Union{Nothing, Tuple{Int64, Int64}}
32         p1::Any,p0::Any = forwardrecurrence_next(n::Int64, A::Any, B::Any, C::Any, x::Float64, p0, p1)::Any,p1::Any
33     end
34     p0::Any,p1::Any::Tuple{Any, Any}
35 end
Select a call to descend into or ↩ to ascend. [q]uit. [b]ookmark.
Toggles: [w]arn, [h]ide type-stable statements, [t]ype annotations, [s]yntax highlight for Source/LLVM/Native, [j]ump to source always, [v]scode: inlay types, [V]scode: diagnostics.
Show: [S]ource code, [A]ST, [T]yped code, [L]LVM IR, [N]ative code
Actions: [E]dit source code, [R]evise and redisplay
 • runtime eltype(A::Any)
   runtime eltype(B::Any)
   runtime eltype(C::Any)
   promote_type(eltype(A::Any)::Any, eltype(B::Any)::Any, eltype(C::Any)::Any, typeof(x::Float64)::Type{Float64})
   runtime convert(T::Any, μ::Float64)
   N::Int64 == 0
   runtime zero(T::Any)
   runtime A::Any[1]
   runtime B::Any[1]
v  runtime muladd(A::Any[1]::Any,x::Float64,B::Any[1]::Any)

Maybe RecurrenceRelationshipArrays needs some explicit type annotations to really help with cases like this. I took a look but there's a lot of machinery that I'm not familiar with at this point

[TSGut]

That's a Lanczos case, how about the integer ones? Asking because I would be happy to look into the decomposition methods but less inclined to spend time improving the Lanczos implementation.

[DanielVandH]

It's independent of the parameters (t, a, b, c) since jacobimatrix is type unstable for all parameters since it's not concretely typed which might be the real reason

[DanielVandH]

Ah, yes it is the real reason. If I make it concretely typed than everything is inferred (obviously, since abstract types of course won't be (: ). So I guess there would be two possible solutions

• Give X a concrete type in struct SemiclassicalJacobi. This would give a compilation penalty though which would be why it's not concretely typed to begin with
• Add more type annotations (or function barriers and just let dynamic dispatch do its thing, putting type annotations only in the outermost function) to RecurrenceRelationshipArrays. This will still give some internal type instabilities, but the final result won't be and the allocations could be significantly reduced.

I guess the latter is preferable since it'll probably help some other packages.

I think this used to be type stable at least a couple months ago when I last checked.. so maybe something else is going on e.g. a function barrier broke somewhere