SemiclassicalJacobiFamily for a range of SemiclassicalJacobi (#45)

* Fix typo in SemiclassicalJacobi{T} constructor

* v0.2.1

* Add SemiclassicalJacobiFamily

* SemiclassicalJacobi{T} family

* D * x^a * (1-x)^b

* ab half-weighted

* two-weighted derivatives done!

* Update test_derivative.jl

* Update runtests.jl

* divmul to avoid recomputation of SemiclassicalJacobi

* Add divmul

* more divmul

* fix em code

* Update derivatives.jl

* fix bugs

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "SemiclassicalOrthogonalPolynomials"
 uuid = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

examples/equilibriummeasure.jl

@@ -4,9 +4,9 @@ import SemiclassicalOrthogonalPolynomials: Weighted
 import ClassicalOrthogonalPolynomials: associated, affine
 
 Base.floatmin(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T,V,N}(floatmin(V))
+Base.big(d::Dual{T,V,N}) where {T,V,N} = Dual{T}(big(d.value), ForwardDiff.Partials(map(big,d.partials.values)))
 
-
-####
+#### 
 #
 # We first do a single interval.
 # Equilibrium measures for a symmetric potential
@@ -61,7 +61,6 @@ for _ = 1:10
 end
 
 plot(equilibrium(b))
-ClassicalOrthogonalPolynomials.plotgrid(equilibrium(b))
 
 
 #####
@@ -99,8 +98,8 @@ function equilibriumcoefficients(P,a,b)
 end
 function equilibriumendpointvalues(ab::SVector{2})
     a,b = ab
-    P = TwoBandJacobi(a/b, -1/2, -1/2, 1/2)
-    P[[a/b,1],:] * equilibriumcoefficients(P,a,b)
+    P = TwoBandJacobi(a/b, -one(a)/2, -one(a)/2, one(a)/2)
+    Vector(P[[a/b,1],:] * equilibriumcoefficients(P,a,b))
 end
 
 function equilibrium(ab)
@@ -109,7 +108,7 @@ function equilibrium(ab)
     Weighted(P) * equilibriumcoefficients(P,a,b)
 end
 
-(a,b) = ab = SVector(2.,3.)
+ab = SVector(2.,3.)
 ab -= jacobian(equilibriumendpointvalues,ab) \ equilibriumendpointvalues(ab)
 
 

src/SemiclassicalOrthogonalPolynomials.jl

@@ -46,7 +46,7 @@ function getindex(P::SemiclassicalJacobiWeight, x::Real)
 end
 
 function sum(P::SemiclassicalJacobiWeight{T}) where T
-    (t,a,b,c) = BigFloat.((P.t,P.a,P.b,P.c))
+    (t,a,b,c) = map(big, map(float, (P.t,P.a,P.b,P.c)))
     return convert(T, t^c * beta(1+a,1+b) * _₂F₁general2(1+a,-c,2+a+b,1/t))
 end
 
@@ -163,7 +163,7 @@ end
 
 SemiclassicalJacobi(t, a, b, c, P::SemiclassicalJacobi) = SemiclassicalJacobi(t, a, b, c, semiclassical_jacobimatrix(P, a, b, c))
 SemiclassicalJacobi(t, a, b, c) = SemiclassicalJacobi(t, a, b, c, semiclassical_jacobimatrix(t, a, b, c))
-SemiclassicalJacobi{T}(t, a, b, c) where T = SemiclassicalJacobi(convert(T,t), convert(T,a), convert(T,b), convert(T,b))
+SemiclassicalJacobi{T}(t, a, b, c) where T = SemiclassicalJacobi(convert(T,t), convert(T,a), convert(T,b), convert(T,c))
 
 WeightedSemiclassicalJacobi(t, a, b, c, P...) = SemiclassicalJacobiWeight(t, a, b, c) .* SemiclassicalJacobi(t, a, b, c, P...)
 
@@ -399,35 +399,66 @@ convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:a,T,<:Semiclassic
 convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:b,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, zero(T),Q.P.b,zero(T)) .* Q.P
 convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:c,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, zero(T),zero(T),Q.P.c) .* Q.P
 
+convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:ab,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, Q.P.a,Q.P.b,zero(T)) .* Q.P
+convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:bc,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, zero(T),Q.P.b,Q.P.c) .* Q.P
+convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:ac,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, Q.P.a,zero(T),Q.P.c) .* Q.P
+
+
 include("derivatives.jl")
 
 
 ###
 # Hierarchy
+#
+# here we build the operators lazily
 ###
 
-function Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, ar::AbstractUnitRange, b::Number, c::Number)
-    Ps = [SemiclassicalJacobi(t, first(ar), b, c)]
-    for a in ar[2:end]
-        push!(Ps, SemiclassicalJacobi(t, a, b, c, Ps[end]))
-    end
-    Ps
+mutable struct SemiclassicalJacobiFamily{T, A, B, C} <: AbstractCachedVector{SemiclassicalJacobi{T}}
+    data::Vector{SemiclassicalJacobi{T}}
+    t::T
+    a::A
+    b::B
+    c::C
+    datasize::Tuple{Int}
 end
 
-function Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, br::AbstractUnitRange, c::Number)
-    Ps = [SemiclassicalJacobi(t, a, first(br), c)]
-    for b in br[2:end]
-        push!(Ps, SemiclassicalJacobi(t, a, b, c, Ps[end]))
-    end
-    Ps
+size(P::SemiclassicalJacobiFamily) = (max(length(P.a), length(P.b), length(P.c)),)
+
+_checkrangesizes() = ()
+_checkrangesizes(a::Number, b...) = _checkrangesizes(b...)
+_checkrangesizes(a, b...) = (length(a), _checkrangesizes(b...)...)
+
+_isequal() = true
+_isequal(a) = true
+_isequal(a,b,c...) = a == b && _isequal(b,c...)
+
+checkrangesizes(a...) = _isequal(_checkrangesizes(a...)...) || throw(DimensionMismatch())
+
+function SemiclassicalJacobiFamily{T}(data::Vector, t, a, b, c) where T
+    checkrangesizes(a, b, c)
+    SemiclassicalJacobiFamily{T,typeof(a),typeof(b),typeof(c)}(data, t, a, b, c, (length(data),))
 end
 
-function Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, b::Number, cr::AbstractUnitRange)
-    Ps = [SemiclassicalJacobi(t, a, b, first(cr))]
-    for c in cr[2:end]
-        push!(Ps, SemiclassicalJacobi(t, a, b, c, Ps[end]))
+SemiclassicalJacobiFamily(t, a, b, c) = SemiclassicalJacobiFamily{float(promote_type(typeof(t),eltype(a),eltype(b),eltype(c)))}(t, a, b, c)
+SemiclassicalJacobiFamily{T}(t, a, b, c) where T = SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c))], t, a, b, c)
+
+Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, b::Number, c::Number) = SemiclassicalJacobi(t, a, b, c)
+Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Number, b::Number, c::Number) where T = SemiclassicalJacobi{T}(t, a, b, c)
+Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) = 
+    SemiclassicalJacobiFamily(t, a, b, c)
+Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) where T = 
+    SemiclassicalJacobiFamily{T}(t, a, b, c)
+
+
+_broadcast_getindex(a,k) = a[k]
+_broadcast_getindex(a::Number,k) = a
+
+function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily, inds::AbstractUnitRange)
+    t,a,b,c = P.t,P.a,P.b,P.c
+    for k in inds
+        P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), P.data[k-1])
     end
-    Ps
+    P
 end
 
 include("twoband.jl")

src/derivatives.jl

@@ -28,43 +28,66 @@ function LazyArrays.cache_filldata!(K::MulAddAccumulate, inds)
     end
 end
 
-@simplify function *(D::Derivative, P::SemiclassicalJacobi)
-    Q = SemiclassicalJacobi(P.t, P.a+1,P.b+1,P.c+1,P)
+"""
+    divmul(A, B, C)
+
+is equivalent to A \\ (B*C)
+"""
+function divmul(Q::SemiclassicalJacobi, D::Derivative, P::SemiclassicalJacobi)
     A,B,_ = recurrencecoefficients(P)
     α,β,_ = recurrencecoefficients(Q)
 
     d = AccumulateAbstractVector(*, A ./ Vcat(1,α))
     v1 = AccumulateAbstractVector(+, B ./ A)
     v2 = MulAddAccumulate(Vcat(0,0,α[2:∞]) ./ α, Vcat(0,β ./ α) ./ α);
     v3 = AccumulateAbstractVector(*, Vcat(A[1]A[2], A[3:∞] ./ α))
-    Q * (_BandedMatrix(Vcat(((1:∞) .* d)', (((1:∞) .* (v1 .+ B[2:end]./A[2:end]) .- (2:∞) .* (α .* v2 .+ β ./ α)) .* v3)'), ∞, 2,-1))'
+    _BandedMatrix(Vcat(((1:∞) .* d)', (((1:∞) .* (v1 .+ B[2:end]./A[2:end]) .- (2:∞) .* (α .* v2 .+ β ./ α)) .* v3)'), ∞, 2,-1)'
+end
+
+@simplify function *(D::Derivative, P::SemiclassicalJacobi)
+    Q = SemiclassicalJacobi(P.t, P.a+1,P.b+1,P.c+1,P)
+    Q * divmul(Q, D, P)
+end
+
+
+
+
+function divmul(wP::Weighted{<:Any,<:SemiclassicalJacobi}, D::Derivative, wQ::Weighted{<:Any,<:SemiclassicalJacobi})
+    Q,P = wQ.P,wP.P
+    ((-sum(orthogonalityweight(Q))/sum(orthogonalityweight(P))) * (Q \ (D * P))')
 end
 
 @simplify function *(D::Derivative, wQ::Weighted{<:Any,<:SemiclassicalJacobi})
-    Q = wQ.P
-    P = SemiclassicalJacobi(Q.t, Q.a-1,Q.b-1,Q.c-1)
-    Weighted(P) * ((-sum(orthogonalityweight(Q))/sum(orthogonalityweight(P))) * (Q \ (D * P))')
+    wP = Weighted(SemiclassicalJacobi(wQ.P.t, wQ.P.a-1,wQ.P.b-1,wQ.P.c-1))
+    wP * divmul(wP, D, wQ)
 end
 
-@simplify function *(D::Derivative, HP::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi})
+
+##
+# One-Weighted
+##
+
+function divmul(HQ::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi}, D::Derivative, HP::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi})
+    Q = HQ.P
     P = HP.P
-    Q = SemiclassicalJacobi(P.t, P.a-1, P.b+1, P.c+1)
+    t = P.t
     a = Q.a
     A,B,C = recurrencecoefficients(P)
     α,β,γ = recurrencecoefficients(Q)
     d = AccumulateAbstractVector(*, A ./ α)
     v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
     v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
 
-    HalfWeighted{:a}(Q) * _BandedMatrix(
+    _BandedMatrix(
         Vcat(
         ((a:∞) .* v2 .- ((a+1):∞) .* Vcat(1,v1[2:end] .* d))',
         (((a+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
 end
 
-@simplify function *(D::Derivative, HP::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi})
+function divmul(HQ::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi}, D::Derivative, HP::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi})
+    Q = HQ.P
     P = HP.P
-    Q = SemiclassicalJacobi(P.t, P.a+1, P.b-1, P.c+1)
+    t = P.t
     b = Q.b
     A,B,C = recurrencecoefficients(P)
     α,β,γ = recurrencecoefficients(Q)
@@ -73,26 +96,124 @@ end
     v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
     v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
 
-    HalfWeighted{:b}(Q) * _BandedMatrix(
+    _BandedMatrix(
         Vcat(
         (-(b:∞) .* v2 .+ ((b+1):∞) .* Vcat(1,v1[2:end] .* d) .+ Vcat(0,(1:∞) .* d2))',
         (-((b+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
 end
 
-@simplify function *(D::Derivative, HP::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi})
+function divmul(HQ::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi}, D::Derivative, HP::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi})
+    Q = HQ.P
     P = HP.P
     t = P.t
-    Q = SemiclassicalJacobi(t, P.a+1, P.b+1, P.c-1)
     c = Q.c
     A,B,C = recurrencecoefficients(P)
     α,β,γ = recurrencecoefficients(Q)
     d = AccumulateAbstractVector(*, A ./ α)
     d2 = AccumulateAbstractVector(*, A ./ Vcat(1,α))
     v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
     v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
-
-    HalfWeighted{:c}(Q) * _BandedMatrix(
+    _BandedMatrix(
         Vcat(
         (-(c:∞) .* v2 .+ ((c+1):∞) .* Vcat(1,v1[2:end] .* d) .+ Vcat(0,(t:t:∞) .* d2))',
         (-((c+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
-end
+end
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    t = P.t
+    HQ = HalfWeighted{:a}(SemiclassicalJacobi(t, P.a-1, P.b+1, P.c+1))
+    HQ * divmul(HQ, D, HP)
+end
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    t = P.t
+    HQ = HalfWeighted{:b}(SemiclassicalJacobi(t, P.a+1, P.b-1, P.c+1))
+    HQ * divmul(HQ, D, HP)
+end
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    t = P.t
+    HQ = HalfWeighted{:c}(SemiclassicalJacobi(t, P.a+1, P.b+1, P.c-1))
+    HQ * divmul(HQ, D, HP)
+end
+
+##
+# Double-Weighted
+##
+
+function divmul(HQ::HalfWeighted{:ab}, D::Derivative, HP::HalfWeighted{:ab})
+    Q = HQ.P
+    P = HP.P
+    A,B,_ = recurrencecoefficients(P)
+    α,β,_ = recurrencecoefficients(Q)
+    a,b = Q.a,Q.b
+
+    d = AccumulateAbstractVector(*, Vcat(1,A) ./ α)
+    e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
+    f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
+    g = cumsum(β ./ α)
+    _BandedMatrix(Vcat((((a+1):∞) .* e .- ((b+a+1):∞).*f .+ ((a+b+2):∞) .* e .* g )',
+                           (-((a+b+2):∞)  .* d)'),ℵ₀,1,0)
+end
+
+
+function divmul(HQ::HalfWeighted{:bc}, D::Derivative, HP::HalfWeighted{:bc})
+    Q = HQ.P
+    P = HP.P
+    t = P.t
+    A,B,_ = recurrencecoefficients(P)
+    α,β,_ = recurrencecoefficients(Q)
+    b,c = Q.b,Q.c
+
+    d = AccumulateAbstractVector(*, Vcat(1,A) ./ α)
+    e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
+    f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
+    g = cumsum(β ./ α)
+    _BandedMatrix(Vcat((-((t+1)* (0:∞) .+ (t*(b+1) + c+1)) .* e .+ ((c+b+1):∞).*f .- ((b+c+2):∞) .* e .* g )',
+                        (((b+c+2):∞)  .* d)'),ℵ₀,1,0)
+end
+
+function divmul(HQ::HalfWeighted{:ac}, D::Derivative, HP::HalfWeighted{:ac})
+    Q = HQ.P
+    P = HP.P
+    t = P.t
+    A,B,_ = recurrencecoefficients(P)
+    α,β,_ = recurrencecoefficients(Q)
+    a,c = Q.a,Q.c
+
+    d = AccumulateAbstractVector(*, Vcat(1,A) ./ α)
+    e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
+    f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
+    g = cumsum(β ./ α)
+    _BandedMatrix(Vcat((t* ((a+1):∞) .* e .- ((c+a+1):∞).*f .+ ((a+c+2):∞) .* e .* g )',
+            (-((a+c+2):∞)  .* d)'),ℵ₀,1,0)
+end
+
+
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:ab,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    t = P.t
+    HQ = HalfWeighted{:ab}(SemiclassicalJacobi(t, P.a-1,P.b-1,P.c+1))
+    HQ * divmul(HQ, D, HP)
+end
+
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:ac,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    t = P.t
+    HQ = HalfWeighted{:ac}(SemiclassicalJacobi(t, P.a-1,P.b+1,P.c-1))
+    HQ * divmul(HQ, D, HP)
+end
+
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:bc,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    t = P.t
+    HQ = HalfWeighted{:bc}(SemiclassicalJacobi(t, P.a+1,P.b-1,P.c-1))
+    HQ * divmul(HQ, D, HP)
+end
+