Skip to content

Commit 0d60c96

Browse files
dlfivefiftydlfivefifty
committed
move family code to seperate file
1 parent 15a0db4 commit 0d60c96

File tree

2 files changed

+147
-149
lines changed

2 files changed

+147
-149
lines changed

‎src/SemiclassicalOrthogonalPolynomials.jl

+1-149
Original file line numberDiff line numberDiff line change
@@ -607,155 +607,7 @@ convert(::Type{WeightedBasis}, Q::HalfWeighted{:bc,T,<:SemiclassicalJacobi}) whe
607607
convert(::Type{WeightedBasis}, Q::HalfWeighted{:ac,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, Q.P.a,zero(T),Q.P.c) .* Q.P
608608

609609
include("derivatives.jl")
610-
611-
612-
###
613-
# Hierarchy
614-
#
615-
# here we build the operators lazily
616-
###
617-
618-
mutable struct SemiclassicalJacobiFamily{T, A, B, C} <: AbstractCachedVector{SemiclassicalJacobi{T}}
619-
data::Vector{SemiclassicalJacobi{T}}
620-
t::T
621-
a::A
622-
b::B
623-
c::C
624-
datasize::Tuple{Int}
625-
end
626-
627-
isnormalized(J::SemiclassicalJacobi) = J.b ≠ -1 # there is no normalisation for b == -1
628-
size(P::SemiclassicalJacobiFamily) = (max(length(P.a), length(P.b), length(P.c)),)
629-
630-
_checkrangesizes() = ()
631-
_checkrangesizes(a::Number, b...) = _checkrangesizes(b...)
632-
_checkrangesizes(a, b...) = (length(a), _checkrangesizes(b...)...)
633-
634-
_isequal() = true
635-
_isequal(a) = true
636-
_isequal(a,b,c...) = a == b && _isequal(b,c...)
637-
638-
checkrangesizes(a...) = _isequal(_checkrangesizes(a...)...) || throw(DimensionMismatch())
639-
640-
function SemiclassicalJacobiFamily{T}(data::Vector, t, a, b, c) where T
641-
checkrangesizes(a, b, c)
642-
SemiclassicalJacobiFamily{T,typeof(a),typeof(b),typeof(c)}(data, t, a, b, c, (length(data),))
643-
end
644-
645-
function _getsecondifpossible(v)
646-
length(v) > 1 && return v[2]
647-
return v[1]
648-
end
649-
650-
SemiclassicalJacobiFamily(t, a, b, c) = SemiclassicalJacobiFamily{float(promote_type(typeof(t),eltype(a),eltype(b),eltype(c)))}(t, a, b, c)
651-
function SemiclassicalJacobiFamily{T}(t, a, b, c) where T
652-
# We need to start with a hierarchy containing two entries
653-
return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, _getsecondifpossible(a), _getsecondifpossible(b), _getsecondifpossible(c))], t, a, b, c)
654-
end
655-
656-
Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, b::Number, c::Number) = SemiclassicalJacobi(t, a, b, c)
657-
Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Number, b::Number, c::Number) where T = SemiclassicalJacobi{T}(t, a, b, c)
658-
Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) =
659-
SemiclassicalJacobiFamily(t, a, b, c)
660-
Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) where T =
661-
SemiclassicalJacobiFamily{T}(t, a, b, c)
662-
663-
664-
_broadcast_getindex(a,k) = a[k]
665-
_broadcast_getindex(a::Number,k) = a
666-
667-
function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily, inds::AbstractUnitRange)
668-
t,a,b,c = P.t,P.a,P.b,P.c
669-
isrange = P.b isa AbstractUnitRange
670-
for k in inds
671-
# If P.data[k-2] is not normalised (aka b = -1), cholesky fails. With the current design, this is only a problem if P.b
672-
# is a range since we can translate between polynomials that both have b = -1.
673-
Pprev = (isrange && P.b[k-2] == -1) ? P.data[k-1] : P.data[k-2] # isrange && P.b[k-2] == -1 could also be !isnormalized(P.data[k-2])
674-
P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), Pprev)
675-
end
676-
P
677-
end
678-
679-
###
680-
# here we construct hierarchies of c weight sums by means of contiguous recurrence relations
681-
###
682-
683-
""""
684-
A SemiclassicalJacobiCWeightFamily
685-
686-
is a vector containing a sequence of weights of the form `x^a * (1-x)^b * (t-x)^c` where `a` and `b` are scalars and `c` is a range of values with integer spacing; where `x in 0..1`. It is automatically generated when calling `SemiclassicalJacobiWeight.(t,a,b,cmin:cmax)`.
687-
"""
688-
struct SemiclassicalJacobiCWeightFamily{T, C} <: AbstractVector{SemiclassicalJacobiWeight{T}}
689-
data::Vector{SemiclassicalJacobiWeight{T}}
690-
t::T
691-
a::T
692-
b::T
693-
c::C
694-
datasize::Tuple{Int}
695-
end
696-
697-
getindex(W::SemiclassicalJacobiCWeightFamily, inds) = getindex(W.data, inds)
698-
699-
size(W::SemiclassicalJacobiCWeightFamily) = (length(W.c),)
700-
701-
function SemiclassicalJacobiCWeightFamily{T}(data::Vector, t, a, b, c) where T
702-
checkrangesizes(a, b, c)
703-
SemiclassicalJacobiCWeightFamily{T,typeof(c)}(data, t, a, b, c, (length(data),))
704-
end
705-
706-
SemiclassicalJacobiCWeightFamily(t, a, b, c) = SemiclassicalJacobiCWeightFamily{float(promote_type(typeof(t),eltype(a),eltype(b),eltype(c)))}(t, a, b, c)
707-
function SemiclassicalJacobiCWeightFamily{T}(t::Number, a::Number, b::Number, c::Union{AbstractUnitRange,Number}) where T
708-
return SemiclassicalJacobiCWeightFamily{T}(SemiclassicalJacobiWeight.(t,a:a,b:b,c), t, a, b, c)
709-
end
710-
711-
Base.broadcasted(::Type{SemiclassicalJacobiWeight}, t::Number, a::Number, b::Number, c::Union{AbstractUnitRange,Number}) =
712-
SemiclassicalJacobiCWeightFamily(t, a, b, c)
713-
714-
_unweightedsemiclassicalsum(a,b,c,t) = pFq((a+1,-c),(a+b+2, ), 1/t)
715-
716-
function Base.broadcasted(::typeof(sum), W::SemiclassicalJacobiCWeightFamily{T}) where T
717-
a = W.a; b = W.b; c = W.c; t = W.t;
718-
cmin = minimum(c); cmax = maximum(c);
719-
@assert isinteger(cmax) && isinteger(cmin)
720-
# This is needed at high parameter values.
721-
# Manually setting setprecision(2048) allows accurate computation even for very high c.
722-
t,a,b = convert(BigFloat,t),convert(BigFloat,a),convert(BigFloat,b)
723-
F = zeros(BigFloat,cmax+1)
724-
F[1] = _unweightedsemiclassicalsum(a,b,0,t) # c=0
725-
cmax == 0 && return abs.(convert.(T,t.^c.*exp(loggamma(a+1)+loggamma(b+1)-loggamma(a+b+2)).*getindex(F,1:1)))
726-
F[2] = _unweightedsemiclassicalsum(a,b,1,t) # c=1
727-
@inbounds for n in 1:cmax-1
728-
F[n+2] = ((n-1)/t+1/t-n)/(n+a+b+2)*F[n]+(a+b+4+2*n-2-(n+a+1)/t)/(n+a+b+2)*F[n+1]
729-
end
730-
return abs.(convert.(T,t.^c.*exp(loggamma(a+1)+loggamma(b+1)-loggamma(a+b+2)).*getindex(F,W.c.+1)))
731-
end
732-
733-
""""
734-
sumquotient(wP, wQ) computes sum(wP)/sum(wQ) by taking into account cancellations, allowing more stable computations for high weight parameters.
735-
"""
736-
function sumquotient(wP::SemiclassicalJacobiWeight{T},wQ::SemiclassicalJacobiWeight{T}) where T
737-
@assert wP.t ≈ wQ.t
738-
@assert isinteger(wP.c) && isinteger(wQ.c)
739-
a = wP.a; b = wP.b; c = Int(wP.c); t = wP.t;
740-
# This is needed at high parameter values.
741-
t,a,b = convert(BigFloat,t),convert(BigFloat,a),convert(BigFloat,b)
742-
F = zeros(BigFloat,max(2,c+1))
743-
F[1] = _unweightedsemiclassicalsum(a,b,0,t) # c=0
744-
F[2] = _unweightedsemiclassicalsum(a,b,1,t) # c=1
745-
@inbounds for n in 1:c-1
746-
F[n+2] = ((n-1)/t+1/t-n)/(n+a+b+2)*F[n]+(a+b+4+2*n-2-(n+a+1)/t)/(n+a+b+2)*F[n+1]
747-
end
748-
a = wQ.a; b = wQ.b; c = Int(wQ.c);
749-
t,a,b = convert(BigFloat,t),convert(BigFloat,a),convert(BigFloat,b)
750-
G = zeros(BigFloat,max(2,c+1))
751-
G[1] = _unweightedsemiclassicalsum(a,b,0,t) # c=0
752-
G[2] = _unweightedsemiclassicalsum(a,b,1,t) # c=1
753-
@inbounds for n in 1:c-1
754-
G[n+2] = ((n-1)/t+1/t-n)/(n+a+b+2)*G[n]+(a+b+4+2*n-2-(n+a+1)/t)/(n+a+b+2)*G[n+1]
755-
end
756-
return abs.(convert.(T,t.^(Int(wP.c)-c).*exp(loggamma(wP.a+1)+loggamma(wP.b+1)-loggamma(wP.a+wP.b+2)-loggamma(a+1)-loggamma(b+1)+loggamma(a+b+2))*F[Int(wP.c)+1]/G[c+1]))
757-
end
758-
610+
include("family.jl")
759611
include("neg1c.jl")
760612
include("deprecated.jl")
761613

‎src/family.jl

+146
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
###
2+
# Hierarchy
3+
#
4+
# here we build the operators lazily
5+
###
6+
7+
mutable struct SemiclassicalJacobiFamily{T, A, B, C} <: AbstractCachedVector{SemiclassicalJacobi{T}}
8+
data::Vector{SemiclassicalJacobi{T}}
9+
t::T
10+
a::A
11+
b::B
12+
c::C
13+
datasize::Tuple{Int}
14+
end
15+
16+
isnormalized(J::SemiclassicalJacobi) = J.b ≠ -1 # there is no normalisation for b == -1
17+
size(P::SemiclassicalJacobiFamily) = (max(length(P.a), length(P.b), length(P.c)),)
18+
19+
_checkrangesizes() = ()
20+
_checkrangesizes(a::Number, b...) = _checkrangesizes(b...)
21+
_checkrangesizes(a, b...) = (length(a), _checkrangesizes(b...)...)
22+
23+
_isequal() = true
24+
_isequal(a) = true
25+
_isequal(a,b,c...) = a == b && _isequal(b,c...)
26+
27+
checkrangesizes(a...) = _isequal(_checkrangesizes(a...)...) || throw(DimensionMismatch())
28+
29+
function SemiclassicalJacobiFamily{T}(data::Vector, t, a, b, c) where T
30+
checkrangesizes(a, b, c)
31+
SemiclassicalJacobiFamily{T,typeof(a),typeof(b),typeof(c)}(data, t, a, b, c, (length(data),))
32+
end
33+
34+
function _getsecondifpossible(v)
35+
length(v) > 1 && return v[2]
36+
return v[1]
37+
end
38+
39+
SemiclassicalJacobiFamily(t, a, b, c) = SemiclassicalJacobiFamily{float(promote_type(typeof(t),eltype(a),eltype(b),eltype(c)))}(t, a, b, c)
40+
function SemiclassicalJacobiFamily{T}(t, a, b, c) where T
41+
# We need to start with a hierarchy containing two entries
42+
return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, _getsecondifpossible(a), _getsecondifpossible(b), _getsecondifpossible(c))], t, a, b, c)
43+
end
44+
45+
Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, b::Number, c::Number) = SemiclassicalJacobi(t, a, b, c)
46+
Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Number, b::Number, c::Number) where T = SemiclassicalJacobi{T}(t, a, b, c)
47+
Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) =
48+
SemiclassicalJacobiFamily(t, a, b, c)
49+
Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) where T =
50+
SemiclassicalJacobiFamily{T}(t, a, b, c)
51+
52+
53+
_broadcast_getindex(a,k) = a[k]
54+
_broadcast_getindex(a::Number,k) = a
55+
56+
function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily, inds::AbstractUnitRange)
57+
t,a,b,c = P.t,P.a,P.b,P.c
58+
isrange = P.b isa AbstractUnitRange
59+
for k in inds
60+
# If P.data[k-2] is not normalised (aka b = -1), cholesky fails. With the current design, this is only a problem if P.b
61+
# is a range since we can translate between polynomials that both have b = -1.
62+
Pprev = (isrange && P.b[k-2] == -1) ? P.data[k-1] : P.data[k-2] # isrange && P.b[k-2] == -1 could also be !isnormalized(P.data[k-2])
63+
P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), Pprev)
64+
end
65+
P
66+
end
67+
68+
###
69+
# here we construct hierarchies of c weight sums by means of contiguous recurrence relations
70+
###
71+
72+
""""
73+
A SemiclassicalJacobiCWeightFamily
74+
75+
is a vector containing a sequence of weights of the form `x^a * (1-x)^b * (t-x)^c` where `a` and `b` are scalars and `c` is a range of values with integer spacing; where `x in 0..1`. It is automatically generated when calling `SemiclassicalJacobiWeight.(t,a,b,cmin:cmax)`.
76+
"""
77+
struct SemiclassicalJacobiCWeightFamily{T, C} <: AbstractVector{SemiclassicalJacobiWeight{T}}
78+
data::Vector{SemiclassicalJacobiWeight{T}}
79+
t::T
80+
a::T
81+
b::T
82+
c::C
83+
datasize::Tuple{Int}
84+
end
85+
86+
getindex(W::SemiclassicalJacobiCWeightFamily, inds) = getindex(W.data, inds)
87+
88+
size(W::SemiclassicalJacobiCWeightFamily) = (length(W.c),)
89+
90+
function SemiclassicalJacobiCWeightFamily{T}(data::Vector, t, a, b, c) where T
91+
checkrangesizes(a, b, c)
92+
SemiclassicalJacobiCWeightFamily{T,typeof(c)}(data, t, a, b, c, (length(data),))
93+
end
94+
95+
SemiclassicalJacobiCWeightFamily(t, a, b, c) = SemiclassicalJacobiCWeightFamily{float(promote_type(typeof(t),eltype(a),eltype(b),eltype(c)))}(t, a, b, c)
96+
function SemiclassicalJacobiCWeightFamily{T}(t::Number, a::Number, b::Number, c::Union{AbstractUnitRange,Number}) where T
97+
return SemiclassicalJacobiCWeightFamily{T}(SemiclassicalJacobiWeight.(t,a:a,b:b,c), t, a, b, c)
98+
end
99+
100+
Base.broadcasted(::Type{SemiclassicalJacobiWeight}, t::Number, a::Number, b::Number, c::Union{AbstractUnitRange,Number}) =
101+
SemiclassicalJacobiCWeightFamily(t, a, b, c)
102+
103+
_unweightedsemiclassicalsum(a,b,c,t) = pFq((a+1,-c),(a+b+2, ), 1/t)
104+
105+
function Base.broadcasted(::typeof(sum), W::SemiclassicalJacobiCWeightFamily{T}) where T
106+
a = W.a; b = W.b; c = W.c; t = W.t;
107+
cmin = minimum(c); cmax = maximum(c);
108+
@assert isinteger(cmax) && isinteger(cmin)
109+
# This is needed at high parameter values.
110+
# Manually setting setprecision(2048) allows accurate computation even for very high c.
111+
t,a,b = convert(BigFloat,t),convert(BigFloat,a),convert(BigFloat,b)
112+
F = zeros(BigFloat,cmax+1)
113+
F[1] = _unweightedsemiclassicalsum(a,b,0,t) # c=0
114+
cmax == 0 && return abs.(convert.(T,t.^c.*exp(loggamma(a+1)+loggamma(b+1)-loggamma(a+b+2)).*getindex(F,1:1)))
115+
F[2] = _unweightedsemiclassicalsum(a,b,1,t) # c=1
116+
@inbounds for n in 1:cmax-1
117+
F[n+2] = ((n-1)/t+1/t-n)/(n+a+b+2)*F[n]+(a+b+4+2*n-2-(n+a+1)/t)/(n+a+b+2)*F[n+1]
118+
end
119+
return abs.(convert.(T,t.^c.*exp(loggamma(a+1)+loggamma(b+1)-loggamma(a+b+2)).*getindex(F,W.c.+1)))
120+
end
121+
122+
""""
123+
sumquotient(wP, wQ) computes sum(wP)/sum(wQ) by taking into account cancellations, allowing more stable computations for high weight parameters.
124+
"""
125+
function sumquotient(wP::SemiclassicalJacobiWeight{T},wQ::SemiclassicalJacobiWeight{T}) where T
126+
@assert wP.t ≈ wQ.t
127+
@assert isinteger(wP.c) && isinteger(wQ.c)
128+
a = wP.a; b = wP.b; c = Int(wP.c); t = wP.t;
129+
# This is needed at high parameter values.
130+
t,a,b = convert(BigFloat,t),convert(BigFloat,a),convert(BigFloat,b)
131+
F = zeros(BigFloat,max(2,c+1))
132+
F[1] = _unweightedsemiclassicalsum(a,b,0,t) # c=0
133+
F[2] = _unweightedsemiclassicalsum(a,b,1,t) # c=1
134+
@inbounds for n in 1:c-1
135+
F[n+2] = ((n-1)/t+1/t-n)/(n+a+b+2)*F[n]+(a+b+4+2*n-2-(n+a+1)/t)/(n+a+b+2)*F[n+1]
136+
end
137+
a = wQ.a; b = wQ.b; c = Int(wQ.c);
138+
t,a,b = convert(BigFloat,t),convert(BigFloat,a),convert(BigFloat,b)
139+
G = zeros(BigFloat,max(2,c+1))
140+
G[1] = _unweightedsemiclassicalsum(a,b,0,t) # c=0
141+
G[2] = _unweightedsemiclassicalsum(a,b,1,t) # c=1
142+
@inbounds for n in 1:c-1
143+
G[n+2] = ((n-1)/t+1/t-n)/(n+a+b+2)*G[n]+(a+b+4+2*n-2-(n+a+1)/t)/(n+a+b+2)*G[n+1]
144+
end
145+
return abs.(convert.(T,t.^(Int(wP.c)-c).*exp(loggamma(wP.a+1)+loggamma(wP.b+1)-loggamma(wP.a+wP.b+2)-loggamma(a+1)-loggamma(b+1)+loggamma(a+b+2))*F[Int(wP.c)+1]/G[c+1]))
146+
end

0 commit comments

Comments
 (0)