Skip to content

Instantly share code, notes, and snippets.

@stla
Last active Oct 24, 2020
Embed
What would you like to do?
Hypergeometric function of a matrix argument in Julia (Koev & Edelman's algorithm)
# this one works only for alpha = 1 !
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico, lambda)
#function _N(dico::Dict{Integer,Integer}, lambda::Vector{Integer})
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if(size(lambda,1) == 0)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(alpha, a, b, kappa, i)
if ((isempty(kappa)) || (kappa[1] == 0))
return 1
end
c = kappa[i] - 1 - (i-1)/alpha
d = alpha*kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f)/(l+h))
return prod1/prod2 * prod3 * prod4
end
function betaratio(kappa::Vector{Integer}, mu::Vector{Integer}, k::Int64, alpha::T) where T<:Real
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(map(x -> x/(x+alpha-1), u))
prod2 = prod(map(x -> (x+alpha)/x, v))
prod3 = prod(map(x -> (x+alpha)/x, w))
return alpha * prod1 * prod2 * prod3
end
function Hypergeo(m, alpha, a, b, x)
function summation(i, z, j)
function jack(k, beta, c, t, mu, Nmu)
for i in (max(k,1) : sum(mu .> 0))
if(size(mu,1) == i || mu[i] > mu[i+1])
d = Nmu
gamma = beta * betaratio(kappa, mu, i, alpha)
mu[i] = mu[i] - 1
println("$mu")
Nmu = _N(dico, mu) - 1
println("$mu")
if mu[i] > 0
jack(i, gamma, c+1, t, mu, Nmu)
else
if Nkappa > 1
f = 0
if(any(mu .> 0))
f = J[Nmu, t-1]
else
f = 1
end
J[Nkappa,t] = J[Nkappa,t] + gamma*x[t]^(c+1) * f
end
end
mu[i] = mu[i] + 1
Nmu = d
end
end
if k == 0
if Nkappa > 1
J[Nkappa,t] = J[Nkappa,t] + J[Nkappa,t-1]
end
else
J[Nkappa,t] = J[Nkappa,t] + beta * x[t]^c * J[Nmu,t-1]
end
end # end jack
r = 0
if i == 1
r = j
else
r = min(kappa[i-1],j)
end
Nkappa = 0
for kappai in 1:r
if(size(kappa,1) >= i)
kappa[i] = kappai
else
push!(kappa, kappai)
end
println("$kappa")
# kappa = vcat(kappa, [kappai])
Nkappa = _N(dico, kappa) - 1
println("Nkappa: $Nkappa")
z = z * _T(alpha, a, b, kappa, i)
if Nkappa > 1 && (size(kappa,1)==1 || kappa[2] == 0)
J[Nkappa,1] = x[1] * (1+alpha*(kappa[1]-1)) * J[Nkappa-1,1]
end
for t in 2:n
jack(0, 1.0, 0, t, kappa, Nkappa)
end
s <- s + z * J[Nkappa,n]
if j > kappai && i < n
summation(i+1, z, j-kappai)
end
end
kappa[i] = 0
end # end summation
#
#
#
n = size(x, 1)
(dico, Pmn) = DictParts(m, n)
s = 1
J = zeros(Pmn, n)
J[1,:] = cumsum(x)
kappa = Int64[]
summation(1, 1.0, m)
return s
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
x = Array{Int64}(undef,3,l)
for j in 1:l
x[1:3,j] = [fin+j, manque-j, j]
end
NewLast = hcat(NewLast, x)
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if isempty(lambda)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(
alpha::R,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
kappa::Vector{<:Integer}
) where {R<:Real}
isComplex = eltype(a) <: Complex || eltype(b) <: Complex
if isempty(kappa) || kappa[1] == 0
return isComplex ? Complex(R(1)) : R(1)
end
i = lastindex(kappa)
c = kappa[i] - 1 - (i-1)/alpha
d = alpha * kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> count(>=(j), kappa), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f) ./ (l+h)) # "/" was an error !
out = prod1/prod2 * prod3 * prod4
return isinf(out) || isnan(out) ? (isComplex ? Complex(R(0)) : R(0)) : out
end
function betaratio(
kappa::Vector{I},
mu::Vector{I},
k::I,
alpha::T
) where {I<:Integer, T<:Real}
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(u ./ (u .+ alpha .- 1))
prod2 = prod((v .+ alpha) ./ v)
prod3 = prod((w .+ alpha) ./ w)
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{T}) where T<:Integer
out = T[]
if !isempty(lambda)
for i in 1:lambda[1]
push!(out, count(>=(i), lambda))
end
end
return out
end
# ------------------------------------------------------------------------------
function hypergeomI(
m::I,
alpha::R,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
n::I,
x::Union{R,Complex{R}}
)::Union{R,Complex{R}} where {I<:Integer, R<:Real}
function summation(
i::Int64,
z::T,
j::I,
kappa::Vector{I}
) where {R<:Real, T<:Union{R,Complex{R}}}
function go(
kappai::I,
zz::T,
s::T
) where {R<:Real, T<:Union{R,Complex{R}}}
if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
return s
else
kappap = vcat(kappa, [kappai])
t = _T(alpha, a, b, filter(x -> x>0, kappap))
zp = zz * x * (n - i + alpha * (kappai-1)) * t
sp = j > kappai && i <= n ?
s + summation(i+1, zp, j-kappai, kappap) : s
spp = sp + zp
go(kappai + 1, zp, spp)
end
end # end go -----------------------------------------------------------
go(I(1), z, T(0))
end # end summation --------------------------------------------------------
isComplex =
eltype(a) <: Complex || eltype(b) <: Complex || eltype(x) <: Complex
the_one = isComplex ? Complex(R(1)) : R(1)
return the_one + summation(0, the_one, m, I[])
end
# ------------------------------------------------------------------------------
function hypergeom(
m::Integer,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
x::Vector{<:Union{R,Complex{R}}},
alpha::Union{Nothing,R}=nothing
)::Union{R,Complex{R}} where R<:Real
if isnothing(alpha)
alpha = R(2)
end
n = length(x)
if all(x .== x[1])
return hypergeomI(m, alpha, a, b, n, x[1])
end
function jack(
k::I,
beta::Union{R,Complex{R}},
c::I,
t::Integer,
mu::Vector{I},
jarray::Array{<:Union{R,Complex{R}},2},
kappa::Vector{I},
nkappa::I
) where {I<:Integer}
for i in ((max(k,1) : count(>(0), mu)))
u = mu[i]
if length(mu) == i || u > mu[i+1]
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
filter!(>(0), mup)
if length(mup) >= i && u > 1 # "!isempty(mup) &&" useless
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
if !isempty(mup) # && mup[1] > 0 #any(mup .> 0)
jarray[nkappa,t] +=
gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
else
jarray[nkappa,t] += gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
jarray[nkappa,t] += jarray[nkappa, t-1]
end
else
jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
end
return jarray
end # end jack -------------------------------------------------------------
function summation(
i::I,
z::T,
j::I,
kappa::Vector{I},
jarray::Array{T,2}
) where {I<:Integer, R<:Real, T<:Union{R,Complex{R}}}
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j ||
i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
jarray[nkappa,1] =
x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
end
for t in 2:n
jarray =
jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go -----------------------------------------------------------
go(I(1), T(z), T(0))
end # end summation --------------------------------------------------------
the_one = eltype(a) <: Complex || eltype(b) <: Complex ||
eltype(x) <: Complex ? Complex(R(1)) : R(1)
(dico, Pmn) = DictParts(m, n)
J = zeros(typeof(the_one), Pmn, n)
J[1,:] = cumsum(x)
return the_one + summation(0, the_one, m, Int64[], J)
end
# univariate -------------------------------------------------------------------
# both work
function hypergeom(
m::Integer,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
x::Union{R,Complex{R}}
)::Union{R,Complex{R}} where {R<:Real}
return hypergeomI(m, R(2), a, b, 1, x)
end
function hhypergeom(
m::Integer,
a::Vector{<:Union{R,T}},
b::Vector{<:Union{R,T}},
x::Union{R,T}
) where {R<:Real,T<:Complex{R}}
return hypergeomI(m, R(2), a, b, 1, x)
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
@inbounds manque = record[2]
@inbounds l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
x = Array{Int64}(undef,3,l)
for j in 1:l
@inbounds x[1:3,j] = [fin+j, manque-j, j]
end
NewLast = hcat(NewLast, x)
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if isempty(lambda)
return 1
end
@inbounds dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(
alpha::R,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
kappa::Vector{<:Integer}
) where {R<:Real}
isComplex = eltype(a) <: Complex || eltype(b) <: Complex
if isempty(kappa) || kappa[1] == 0
return isComplex ? Complex(R(1)) : R(1)
end
i = lastindex(kappa)
@inbounds c = kappa[i] - 1 - (i-1)/alpha
@inbounds d = alpha * kappa[i] - i
@inbounds s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> count(>=(j), kappa), s)
g = e .+ 1
ss = 1:(i-1)
@inbounds f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f) ./ (l+h)) # "/" was an error !
out = prod1/prod2 * prod3 * prod4
return isinf(out) || isnan(out) ? (isComplex ? Complex(R(0)) : R(0)) : out
end
function betaratio(
kappa::Vector{I},
mu::Vector{I},
k::I,
alpha::T
) where {I<:Integer, T<:Real}
@inbounds muk = mu[k]
t = k - alpha * muk
@inbounds u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
@inbounds v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
@inbounds w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(u ./ (u .+ alpha .- 1))
prod2 = prod((v .+ alpha) ./ v)
prod3 = prod((w .+ alpha) ./ w)
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{T}) where T<:Integer
out = T[]
if !isempty(lambda)
for i in 1:lambda[1]
push!(out, count(>=(i), lambda))
end
end
return out
end
# ------------------------------------------------------------------------------
function hypergeomI(
m::I,
alpha::R,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
n::I,
x::Union{R,Complex{R}}
)::Union{R,Complex{R}} where {I<:Integer, R<:Real}
function summation(
i::Int64,
z::T,
j::I,
kappa::Vector{I}
) where {R<:Real, T<:Union{R,Complex{R}}}
function go(
kappai::I,
zz::T,
s::T
) where {R<:Real, T<:Union{R,Complex{R}}}
if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
return s
else
kappap = vcat(kappa, [kappai])
t = _T(alpha, a, b, filter(x -> x>0, kappap))
zp = zz * x * (n - i + alpha * (kappai-1)) * t
sp = j > kappai && i <= n ?
s + summation(i+1, zp, j-kappai, kappap) : s
spp = sp + zp
go(kappai + 1, zp, spp)
end
end # end go -----------------------------------------------------------
go(I(1), z, T(0))
end # end summation --------------------------------------------------------
isComplex =
eltype(a) <: Complex || eltype(b) <: Complex || eltype(x) <: Complex
the_one = isComplex ? Complex(R(1)) : R(1)
return the_one + summation(0, the_one, m, I[])
end
# ------------------------------------------------------------------------------
function hypergeom(
m::Integer,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
x::Vector{<:Union{R,Complex{R}}},
alpha::Union{Nothing,R}=nothing
)::Union{R,Complex{R}} where R<:Real
if isnothing(alpha)
alpha = R(2)
end
n = length(x)
if all(x .== x[1])
return hypergeomI(m, alpha, a, b, n, x[1])
end
function jack(
k::I,
beta::Union{R,Complex{R}},
c::I,
t::Integer,
mu::Vector{I},
jarray::Array{<:Union{R,Complex{R}},2},
kappa::Vector{I},
nkappa::I
) where {I<:Integer}
for i in ((max(k,1) : count(>(0), mu)))
@inbounds u = mu[i]
if length(mu) == i || u > mu[i+1]
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
filter!(>(0), mup)
if length(mup) >= i && u > 1 # "!isempty(mup) &&" useless
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
if !isempty(mup) # && mup[1] > 0 #any(mup .> 0)
@inbounds jarray[nkappa,t] +=
gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
else
@inbounds jarray[nkappa,t] += gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
@inbounds jarray[nkappa,t] += jarray[nkappa, t-1]
end
else
@inbounds jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
end
return jarray
end # end jack -------------------------------------------------------------
function summation(
i::I,
z::T,
j::I,
kappa::Vector{I},
jarray::Array{T,2}
) where {I<:Integer, R<:Real, T<:Union{R,Complex{R}}}
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j ||
i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
@inbounds jarray[nkappa,1] =
x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
end
for t in 2:n
jarray =
jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go -----------------------------------------------------------
go(I(1), T(z), T(0))
end # end summation --------------------------------------------------------
the_one = eltype(a) <: Complex || eltype(b) <: Complex ||
eltype(x) <: Complex ? Complex(R(1)) : R(1)
(dico, Pmn) = DictParts(m, n)
J = zeros(typeof(the_one), Pmn, n)
J[1,:] = cumsum(x)
return the_one + summation(0, the_one, m, Int64[], J)
end
# univariate -------------------------------------------------------------------
# both work
function hypergeom(
m::Integer,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
x::Union{R,Complex{R}}
)::Union{R,Complex{R}} where {R<:Real}
return hypergeomI(m, R(2), a, b, 1, x)
end
function hhypergeom(
m::Integer,
a::Vector{<:Union{R,T}},
b::Vector{<:Union{R,T}},
x::Union{R,T}
) where {R<:Real,T<:Complex{R}}
return hypergeomI(m, R(2), a, b, 1, x)
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico, lambda)
#function _N(dico::Dict{Integer,Integer}, lambda::Vector{Integer})
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if(size(lambda,1) == 0)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(alpha, a, b, kappa, i)
if ((isempty(kappa)) || (kappa[1] == 0))
return 1
end
c = kappa[i] - 1 - (i-1)/alpha
d = alpha*kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f)/(l+h))
return prod1/prod2 * prod3 * prod4
end
#function betaratio(kappa::Vector{Integer}, mu::Vector{Integer}, k::Int64, alpha::T) where T<:Real
function betaratio(kappa, mu, k, alpha)
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(map(x -> x/(x+alpha-1), u))
prod2 = prod(map(x -> (x+alpha)/x, v))
prod3 = prod(map(x -> (x+alpha)/x, w))
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Array{Int64,1})
out = Int64[]
for i in 1:lambda[1]
push!(out, size(filter(x -> x>=i, lambda), 1))
end
return out
end
function summation(a, b, x, dico, n, alpha, i, z, j, kappa, jarray)
function go(kappai, zp, s)
if i == n || i == 0 && kappai > j || i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap, lastindex(kappap)) # !!! _T à modifier
lkappap = size(kappap, 1)
if nkappa > 1 && (lkappap == 1 || kappap[2] == 0)
entry = jarray[nkappa-1,1]
jarray[nkappa,1] = x[1] * (1 + alpha*(kappap[1]-1)) * entry
end
for t in 2:n
jarray = jack(alpha, x, dico, 0, 1.0, 0, t, kappap, jarray, kappap, nkappa)
end
entryp = jarray[nkappa,n]
sp = s + zpp * entryp
if j > kappai && i <= n
spp = summation(a, b, x, dico, n, alpha, i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go
go(1, z, 0.0)
end
function jack(alpha, x, dico, k, beta, c, t, mu, jarray, kappa, nkappa)
for i in ((max(k,1) : sum(mu .> 0)))
u = mu[i]
if(size(mu,1) == i || u > mu[i+1])
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
mup = mup[mup .> 0]
nmu= _N(dico, mup) - 1
if !isempty(mup) && size(mup,1) >= i && u > 1
jarray = jack(alpha, x, dico, i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
entryp = jarray[nkappa, t]
if any(mup .> 0)
entry = jarray[nmu, t-1]
jarray[nkappa,t] = entryp + gamma * entry * x[t]^(c+1)
else
jarray[nkappa,t] = entryp + gamma * x[t]^(c+1)
end
end
end
end
end
entry1 = jarray[nkappa,t]
if k == 0
if nkappa > 1
entry2 = jarray[nkappa, t-1]
jarray[nkappa,t] = entry1 + entry2
end
else
entry2 = jarray[_N(dico, mu)-1, t-1]
jarray[nkappa,t] = entry1 + beta * x[t]^c * entry2
end
return jarray
end
function hypergeom(m, alpha, a, b, x)
n = size(x, 1)
(dico, Pmn) = DictParts(m, n)
J = zeros(Pmn, n)
J[1,:] = cumsum(x)
kappa = Int64[]
return 1.0 + summation(a, b, x, dico, n, alpha, 0, 1.0, m, kappa, J)
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if(size(lambda,1) == 0)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(alpha::T, a::Vector{T}, b::Vector{T}, kappa::Vector{I}) where {I<:Integer, T<:Real}
if ((isempty(kappa)) || (kappa[1] == 0))
return 1
end
i = lastindex(kappa)
c = kappa[i] - 1 - (i-1)/alpha
d = alpha*kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f)/(l+h))
return prod1/prod2 * prod3 * prod4
end
function betaratio(kappa::Vector{I}, mu::Vector{I}, k::I, alpha::T) where {I<:Integer, T<:Real}
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(map(x -> x/(x+alpha-1), u))
prod2 = prod(map(x -> (x+alpha)/x, v))
prod3 = prod(map(x -> (x+alpha)/x, w))
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{Integer})
out = Int64[]
for i in 1:lambda[1]
push!(out, sum(lambda .>= i))# size(filter(x -> x>=i, lambda), 1))
end
return out
end
function hypergeom(m::Integer, alpha::T, a::Vector{T}, b::Vector{T}, x::Vector{T}) where T<:Real
function jack(k::I, beta::T, c::I, t::Integer, mu::Vector{I}, jarray::Array{T,2}, kappa::Vector{I}, nkappa::I) where I<:Integer
for i in ((max(k,1) : sum(mu .> 0)))
u = mu[i]
if(size(mu,1) == i || u > mu[i+1])
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
mup = mup[mup .> 0]
#nmu= _N(dico, mup) - 1
if !isempty(mup) && size(mup,1) >= i && u > 1
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
entryp = jarray[nkappa, t]
if any(mup .> 0)
entry = jarray[_N(dico, mup) - 1, t-1]
jarray[nkappa,t] = entryp + gamma * entry * x[t]^(c+1)
else
jarray[nkappa,t] = entryp + gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
jarray[nkappa,t] = jarray[nkappa,t] + jarray[nkappa, t-1]
end
else
entry = jarray[_N(dico, mu)-1, t-1]
jarray[nkappa,t] = jarray[nkappa,t] + beta * x[t]^c * entry
end
return jarray
end # end jack
function summation(i::I, z::T, j::I, kappa::Vector{I}, jarray::Array{T,2}) where I<:Integer
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j || i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
#lkappap = size(kappap, 1)
if nkappa > 1 && (size(kappap, 1) == 1 || kappap[2] == 0)
entry = jarray[nkappa-1,1]
jarray[nkappa,1] = x[1] * (1 + alpha*(kappap[1]-1)) * entry
end
for t in 2:n
jarray = jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go
go(I(1), T(z), T(0))
end # end summation
n = size(x, 1)
(dico, Pmn) = DictParts(m, n)
J = zeros(T, Pmn, n)
J[1,:] = cumsum(x)
return T(1) + summation(0, T(1), m, Int64[], J)
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if isempty(lambda)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(
alpha::T,
a::Vector{T},
b::Vector{T},
kappa::Vector{I}
) where {I<:Integer, T<:Real}
if isempty(kappa) || kappa[1] == 0
return T(1)
end
i = lastindex(kappa)
c = kappa[i] - 1 - (i-1)/alpha
d = alpha*kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f)/(l+h))
return prod1/prod2 * prod3 * prod4
end
function betaratio(
kappa::Vector{I},
mu::Vector{I},
k::I,
alpha::T
) where {I<:Integer, T<:Real}
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
# prod1 = prod(map(x -> x/(x+alpha-1), u))
# prod2 = prod(map(x -> (x+alpha)/x, v))
# prod3 = prod(map(x -> (x+alpha)/x, w))
prod1 = prod(u ./ (u .+ alpha .- 1))
prod2 = prod((v .+ alpha) ./ v)
prod3 = prod((w .+ alpha) ./ w)
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{T}) where T<:Integer
out = T[]
if !isempty(lambda)
for i in 1:lambda[1]
push!(out, sum(lambda .>= i)) # size(filter(x -> x>=i, lambda), 1))
end
end
return out
end
function hypergeom(
m::Integer,
alpha::T,
a::Vector{T},
b::Vector{T},
x::Vector{T}
) where T<:Real
function jack(
k::I,
beta::T,
c::I,
t::Integer,
mu::Vector{I},
jarray::Array{T,2},
kappa::Vector{I},
nkappa::I
) where I<:Integer
for i in ((max(k,1) : sum(mu .> 0)))
u = mu[i]
if size(mu,1) == i || u > mu[i+1]
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
mup = mup[mup .> 0]
#nmu= _N(dico, mup) - 1
if !isempty(mup) && size(mup,1) >= i && u > 1
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
if any(mup .> 0)
jarray[nkappa,t] +=
gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
else
jarray[nkappa,t] += gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
jarray[nkappa,t] += jarray[nkappa, t-1]
end
else
jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
end
return jarray
end # end jack
function summation(
i::I,
z::T,
j::I,
kappa::Vector{I},
jarray::Array{T,2}
) where I<:Integer
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j ||
i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
if nkappa > 1 && (size(kappap, 1) == 1 || kappap[2] == 0)
jarray[nkappa,1] =
x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
end
for t in 2:n
jarray =
jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go
go(I(1), T(z), T(0))
end # end summation
n = size(x, 1)
(dico, Pmn) = DictParts(m, n)
J = zeros(T, Pmn, n)
J[1,:] = cumsum(x)
return T(1) + summation(0, T(1), m, Int64[], J)
end
# this one allows complex numbers
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if isempty(lambda)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(
alpha::R,
a::Vector{T},
b::Vector{T},
kappa::Vector{I}
) where {R<:Real, I<:Integer, T<:Number}
if isempty(kappa) || kappa[1] == 0
return T(1)
end
i = lastindex(kappa)
c = kappa[i] - 1 - (i-1)/alpha
d = alpha * kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f)/(l+h))
return prod1/prod2 * prod3 * prod4
end
function betaratio(
kappa::Vector{I},
mu::Vector{I},
k::I,
alpha::T
) where {I<:Integer, T<:Real}
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(u ./ (u .+ alpha .- 1))
prod2 = prod((v .+ alpha) ./ v)
prod3 = prod((w .+ alpha) ./ w)
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{T}) where T<:Integer
out = T[]
if !isempty(lambda)
for i in 1:lambda[1]
push!(out, sum(lambda .>= i))
end
end
return out
end
function hypergeom(
m::Integer,
alpha::R,
a::Vector{T},
b::Vector{T},
x::Vector{T}
) where {R<:Real, T<:Number}
function jack(
k::I,
beta::T,
c::I,
t::Integer,
mu::Vector{I},
jarray::Array{T,2},
kappa::Vector{I},
nkappa::I
) where I<:Integer
for i in ((max(k,1) : sum(mu .> 0)))
u = mu[i]
if length(mu) == i || u > mu[i+1]
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
mup = mup[mup .> 0]
if !isempty(mup) && length(mup) >= i && u > 1
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
if any(mup .> 0)
jarray[nkappa,t] +=
gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
else
jarray[nkappa,t] += gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
jarray[nkappa,t] += jarray[nkappa, t-1]
end
else
jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
end
return jarray
end # end jack -------------------------------------------------------------
function summation(
i::I,
z::T,
j::I,
kappa::Vector{I},
jarray::Array{T,2}
) where I<:Integer
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j ||
i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
jarray[nkappa,1] =
x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
end
for t in 2:n
jarray =
jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go -----------------------------------------------------------
go(I(1), T(z), T(0))
end # end summation --------------------------------------------------------
n = length(x)
(dico, Pmn) = DictParts(m, n)
J = zeros(T, Pmn, n)
J[1,:] = cumsum(x)
return T(1) + summation(0, T(1), m, Int64[], J)
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if isempty(lambda)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(
alpha::R,
a::Vector{T},
b::Vector{T},
kappa::Vector{I}
) where {R<:Real, I<:Integer, T<:Number}
if isempty(kappa) || kappa[1] == 0
return T(1)
end
i = lastindex(kappa)
c = kappa[i] - 1 - (i-1)/alpha
d = alpha * kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f)/(l+h))
return prod1/prod2 * prod3 * prod4
end
function betaratio(
kappa::Vector{I},
mu::Vector{I},
k::I,
alpha::T
) where {I<:Integer, T<:Real}
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(u ./ (u .+ alpha .- 1))
prod2 = prod((v .+ alpha) ./ v)
prod3 = prod((w .+ alpha) ./ w)
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{T}) where T<:Integer
out = T[]
if !isempty(lambda)
for i in 1:lambda[1]
push!(out, sum(lambda .>= i))
end
end
return out
end
function hypergeomI(
m::I,
alpha::R,
a::Vector{T},
b::Vector{T},
n::I,
x::T
) where {I<:Integer, R<:Real, T<:Number}
function summation(i::Int64, z::T, j::I, kappa::Vector{I})
function go(kappai::I, zz::T, s::T)
if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
return s
else
kappap = vcat(kappa, [kappai])
t = _T(alpha, a, b, filter(x -> x>0, kappap))
zp = zz * x * (n - i + alpha * (kappai-1)) * t
sp = j > kappai && i <= n ?
s + summation(i+1, zp, j-kappai, kappap) : s
spp = sp + zp
go(kappai + 1, zp, spp)
end
end # end go -----------------------------------------------------------
go(I(1), z, T(0))
end # end summation --------------------------------------------------------
return 1 + summation(0, T(1), m, I[])
end
function hypergeom(
m::Integer,
alpha::R,
a::Vector{T},
b::Vector{T},
x::Vector{T}
) where {R<:Real, T<:Number}
n = length(x)
if all(x .== x[1])
return hypergeomI(m, alpha, a, b, n, x[1])
end
function jack(
k::I,
beta::T,
c::I,
t::Integer,
mu::Vector{I},
jarray::Array{T,2},
kappa::Vector{I},
nkappa::I
) where I<:Integer
for i in ((max(k,1) : sum(mu .> 0)))
u = mu[i]
if length(mu) == i || u > mu[i+1]
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
mup = mup[mup .> 0]
if !isempty(mup) && length(mup) >= i && u > 1
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
if any(mup .> 0)
jarray[nkappa,t] +=
gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
else
jarray[nkappa,t] += gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
jarray[nkappa,t] += jarray[nkappa, t-1]
end
else
jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
end
return jarray
end # end jack -------------------------------------------------------------
function summation(
i::I,
z::T,
j::I,
kappa::Vector{I},
jarray::Array{T,2}
) where I<:Integer
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j ||
i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
jarray[nkappa,1] =
x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
end
for t in 2:n
jarray =
jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go -----------------------------------------------------------
go(I(1), T(z), T(0))
end # end summation --------------------------------------------------------
(dico, Pmn) = DictParts(m, n)
J = zeros(T, Pmn, n)
J[1,:] = cumsum(x)
return T(1) + summation(0, T(1), m, Int64[], J)
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if isempty(lambda)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(
alpha::R,
a::Vector{T},
b::Vector{T},
kappa::Vector{I}
) where {R<:Real, I<:Integer, T<:Number}
if isempty(kappa) || kappa[1] == 0
return T(1)
end
i = lastindex(kappa)
c = kappa[i] - 1 - (i-1)/alpha
d = alpha * kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f)/(l+h))
return prod1/prod2 * prod3 * prod4
end
# _T(1//1, [1//1], [2//1],[1,1,1]) = -1//0
function betaratio(
kappa::Vector{I},
mu::Vector{I},
k::I,
alpha::T
) where {I<:Integer, T<:Real}
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(u ./ (u .+ alpha .- 1))
prod2 = prod((v .+ alpha) ./ v)
prod3 = prod((w .+ alpha) ./ w)
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{T}) where T<:Integer
out = T[]
if !isempty(lambda)
for i in 1:lambda[1]
push!(out, sum(lambda .>= i))
end
end
return out
end
# ------------------------------------------------------------------------------
function hypergeomI(
m::I,
alpha::R,
a::Vector{T},
b::Vector{T},
n::I,
x::T
) where {I<:Integer, R<:Real, T<:Number}
function summation(i::Int64, z::T, j::I, kappa::Vector{I})
function go(kappai::I, zz::T, s::T)
if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
return s
else
kappap = vcat(kappa, [kappai])
#println("kappap: $kappap")
t = _T(alpha, a, b, filter(x -> x>0, kappap))
#println("t: $t") # -1//0 when alpha=1//1 a = [1///1] b = [2//1] kappap = [1,1,1]
www = zz * x * (n - i + alpha * (kappai-1))
#println("www: $www") # 0//1
zp = zz * x * (n - i + alpha * (kappai-1)) * t
#println("zp: $zp")
sp = j > kappai && i <= n ?
s + summation(i+1, zp, j-kappai, kappap) : s
spp = sp + zp
go(kappai + 1, zp, spp)
end
end # end go -----------------------------------------------------------
go(I(1), z, T(0))
end # end summation --------------------------------------------------------
return 1 + summation(0, T(1), m, I[])
end
# hypergeomI(3, 1//1, [1//1], [2//1], 1, 1//2) error division per zero
# ------------------------------------------------------------------------------
function hypergeom(
m::Integer,
a::Vector{T},
b::Vector{T},
x::Vector{T},
alpha::Union{Nothing,R}=nothing
) where {R<:Real, T<:Union{R,Complex{R}}}
if isnothing(alpha)
alpha = R(2)
end
n = length(x)
if all(x .== x[1])
return hypergeomI(m, alpha, a, b, n, x[1])
end
function jack(
k::I,
beta::T,
c::I,
t::Integer,
mu::Vector{I},
jarray::Array{T,2},
kappa::Vector{I},
nkappa::I
) where I<:Integer
for i in ((max(k,1) : sum(mu .> 0)))
u = mu[i]
if length(mu) == i || u > mu[i+1]
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
mup = mup[mup .> 0]
if !isempty(mup) && length(mup) >= i && u > 1
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
if any(mup .> 0)
jarray[nkappa,t] +=
gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
else
jarray[nkappa,t] += gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
jarray[nkappa,t] += jarray[nkappa, t-1]
end
else
jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
end
return jarray
end # end jack -------------------------------------------------------------
function summation(
i::I,
z::T,
j::I,
kappa::Vector{I},
jarray::Array{T,2}
) where I<:Integer
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j ||
i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
jarray[nkappa,1] =
x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
end
for t in 2:n
jarray =
jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go -----------------------------------------------------------
go(I(1), T(z), T(0))
end # end summation --------------------------------------------------------
(dico, Pmn) = DictParts(m, n)
J = zeros(T, Pmn, n)
J[1,:] = cumsum(x)
return T(1) + summation(0, T(1), m, Int64[], J)
end
# univariate -------------------------------------------------------------------
# function hypergeom(
# m::Integer,
# a::Vector{T},
# b::Vector{T},
# x::T,
# alpha::Union{Nothing,R}=nothing # useless - that does not depend on alpha !
# ) where {R<:Real, T<:Union{R,Complex{R}}}
# if isnothing(alpha)
# alpha = R(1)
# end
# return hypergeomI(m, alpha, a, b, 1, x)
# end
function hypergeom(
m::Integer,
a::Vector{T},
b::Vector{T},
x::T
) where {R<:Real, T<:Union{R,Complex{R}}}
return hypergeomI(m, R(2), a, b, 1, x)
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if isempty(lambda)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(
alpha::R,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
kappa::Vector{I}
) where {R<:Real, I<:Integer, T<:Union{Real,Complex{Real}}}
if isempty(kappa) || kappa[1] == 0
return T(1)
end
i = lastindex(kappa)
c = kappa[i] - 1 - (i-1)/alpha
d = alpha * kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f)/(l+h))
return prod1/prod2 * prod3 * prod4
end
# _T(1//1, [1//1], [2//1],[1,1,1]) = -1//0
function betaratio(
kappa::Vector{I},
mu::Vector{I},
k::I,
alpha::T
) where {I<:Integer, T<:Real}
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(u ./ (u .+ alpha .- 1))
prod2 = prod((v .+ alpha) ./ v)
prod3 = prod((w .+ alpha) ./ w)
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{T}) where T<:Integer
out = T[]
if !isempty(lambda)
for i in 1:lambda[1]
push!(out, sum(lambda .>= i))
end
end
return out
end
# ------------------------------------------------------------------------------
function hypergeomI(
m::I,
alpha::R,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
n::I,
x::Union{R,Complex{R}}
)::Union{R,Complex{R}} where {I<:Integer, R<:Real}
function summation(i::Int64, z::T, j::I, kappa::Vector{I}) where {R<:Real,T<:Union{R,Complex{R}}}
function go(kappai::I, zz::T, s::T) where {R<:Real,T<:Union{R,Complex{R}}}
if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
return s
else
kappap = vcat(kappa, [kappai])
#println("kappap: $kappap")
t = _T(alpha, a, b, filter(x -> x>0, kappap))
#println("t: $t") # -1//0 when alpha=1//1 a = [1///1] b = [2//1] kappap = [1,1,1]
www = zz * x * (n - i + alpha * (kappai-1))
#println("www: $www") # 0//1
zp = zz * x * (n - i + alpha * (kappai-1)) * t
#println("zp: $zp")
sp = j > kappai && i <= n ?
s + summation(i+1, zp, j-kappai, kappap) : s
spp = sp + zp
go(kappai + 1, zp, spp)
end
end # end go -----------------------------------------------------------
go(I(1), z, T(0))
end # end summation --------------------------------------------------------
return 1.0 + summation(0, 1.0, m, I[])
end
# hypergeomI(3, 1//1, [1//1], [2//1], 1, 1//2) error division per zero
# ------------------------------------------------------------------------------
function hypergeom(
m::Integer,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
x::Vector{T},
alpha::Union{Nothing,R}=nothing
)::Union{R,Complex{R}} where {R<:Real, T<:Union{R,Complex{R}}}
if isnothing(alpha)
alpha = R(2)
end
n = length(x)
if all(x .== x[1])
return hypergeomI(m, alpha, a, b, n, x[1])
end
function jack(
k::I,
beta::Union{R,Complex{R}},
c::I,
t::Integer,
mu::Vector{I},
jarray::Array{<:Union{R,Complex{R}},2},
kappa::Vector{I},
nkappa::I
) where {I<:Integer}
for i in ((max(k,1) : sum(mu .> 0)))
u = mu[i]
if length(mu) == i || u > mu[i+1]
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
mup = mup[mup .> 0]
if !isempty(mup) && length(mup) >= i && u > 1
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
if any(mup .> 0)
jarray[nkappa,t] +=
gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
else
jarray[nkappa,t] += gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
jarray[nkappa,t] += jarray[nkappa, t-1]
end
else
jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
end
return jarray
end # end jack -------------------------------------------------------------
function summation(
i::I,
z::T,
j::I,
kappa::Vector{I},
jarray::Array{T,2}
) where {I<:Integer, R<:Real, T<:Union{R,Complex{R}}}
#z = convert(typeof(x),z)
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j ||
i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
jarray[nkappa,1] =
x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
end
for t in 2:n
jarray =
jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go -----------------------------------------------------------
go(I(1), T(z), T(0))
end # end summation --------------------------------------------------------
the_one = eltype(a) <: Complex || eltype(b) <: Complex || eltype(x) <: Complex ? Complex(R(1)) : R(1)
the_type = typeof(the_one)
(dico, Pmn) = DictParts(m, n)
#println("theone: $the_one")
J = zeros(the_type, Pmn, n)
J[1,:] = cumsum(x)
return convert(the_type,1) + summation(0, the_one, m, Int64[], J)
end
# univariate -------------------------------------------------------------------
# both work
function hypergeom(
m::Integer,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
x::Union{R,Complex{R}}
)::Union{R,Complex{R}} where {R<:Real}
return hypergeomI(m, R(2), a, b, 1, x)
end
function hhypergeom(
m::Integer,
a::Vector{<:Union{R,T}},
b::Vector{<:Union{R,T}},
x::Union{R,T}
) where {R<:Real,T<:Complex{R}}
return hypergeomI(m, R(2), a, b, 1, x)
end
function DictParts(m::Integer, n::Integer)
D = Dict{Int64,Int64}()
Last = hcat([0,m,m])
fin = 0
for i in 1:n
NewLast = Array{Int64}(undef,3,0)
for record in eachcol(Last)
manque = record[2]
l = min(manque, record[3])
if l > 0
D[record[1]+1] = fin + 1
for j in 1:l
NewLast = hcat(NewLast, [fin+j, manque-j, j])
end
fin += l
end
end
Last = NewLast
end
return (D, fin)
end
function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
if isempty(lambda)
return 1
end
dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat
function _T(
alpha::R,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
kappa::Vector{<:Integer}
) where {R<:Real}
isComplex = eltype(a) <: Complex || eltype(b) <: Complex
if isempty(kappa) || kappa[1] == 0
return isComplex ? Complex(R(1)) : R(1)
end
i = lastindex(kappa)
c = kappa[i] - 1 - (i-1)/alpha
d = alpha * kappa[i] - i
s = collect(1:(kappa[i]-1))
e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
g = e .+ 1
ss = 1:(i-1)
f = alpha .* kappa[ss] - collect(ss) .- d
h = f .+ alpha
l = h .* f
prod1 = prod(a .+ c)
prod2 = prod(b .+ c)
prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
prod4 = prod((l-f) ./ (l+h)) # "/" was an error !
out = prod1/prod2 * prod3 * prod4
return isinf(out) || isnan(out) ? (isComplex ? Complex(R(0)) : R(0)) : out
end
function betaratio(
kappa::Vector{I},
mu::Vector{I},
k::I,
alpha::T
) where {I<:Integer, T<:Real}
muk = mu[k]
t = k - alpha * muk
u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
v = map(i -> t - i + alpha * mu[i], 1:(k-1))
muPrime = dualPartition(mu)
w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
prod1 = prod(u ./ (u .+ alpha .- 1))
prod2 = prod((v .+ alpha) ./ v)
prod3 = prod((w .+ alpha) ./ w)
return alpha * prod1 * prod2 * prod3
end
function dualPartition(lambda::Vector{T}) where T<:Integer
out = T[]
if !isempty(lambda)
for i in 1:lambda[1]
push!(out, sum(lambda .>= i))
end
end
return out
end
# ------------------------------------------------------------------------------
function hypergeomI(
m::I,
alpha::R,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
n::I,
x::Union{R,Complex{R}}
)::Union{R,Complex{R}} where {I<:Integer, R<:Real}
function summation(
i::Int64,
z::T,
j::I,
kappa::Vector{I}
) where {R<:Real, T<:Union{R,Complex{R}}}
function go(
kappai::I,
zz::T,
s::T
) where {R<:Real, T<:Union{R,Complex{R}}}
if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
return s
else
kappap = vcat(kappa, [kappai])
t = _T(alpha, a, b, filter(x -> x>0, kappap))
zp = zz * x * (n - i + alpha * (kappai-1)) * t
sp = j > kappai && i <= n ?
s + summation(i+1, zp, j-kappai, kappap) : s
spp = sp + zp
go(kappai + 1, zp, spp)
end
end # end go -----------------------------------------------------------
go(I(1), z, T(0))
end # end summation --------------------------------------------------------
isComplex =
eltype(a) <: Complex || eltype(b) <: Complex || eltype(x) <: Complex
the_one = isComplex ? Complex(R(1)) : R(1)
return the_one + summation(0, the_one, m, I[])
end
# ------------------------------------------------------------------------------
function hypergeom(
m::Integer,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
x::Vector{<:Union{R,Complex{R}}},
alpha::Union{Nothing,R}=nothing
)::Union{R,Complex{R}} where R<:Real
if isnothing(alpha)
alpha = R(2)
end
n = length(x)
if all(x .== x[1])
return hypergeomI(m, alpha, a, b, n, x[1])
end
function jack(
k::I,
beta::Union{R,Complex{R}},
c::I,
t::Integer,
mu::Vector{I},
jarray::Array{<:Union{R,Complex{R}},2},
kappa::Vector{I},
nkappa::I
) where {I<:Integer}
for i in ((max(k,1) : sum(mu .> 0)))
u = mu[i]
if length(mu) == i || u > mu[i+1]
gamma = beta * betaratio(kappa, mu, i, alpha)
mup = copy(mu)
mup[i] = u - 1
mup = mup[mup .> 0]
if length(mup) >= i && u > 1 # "!isempty(mup) &&" useless
jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
else
if nkappa > 1
if !isempty(mup) && mup[1] > 0 #any(mup .> 0)
jarray[nkappa,t] +=
gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
else
jarray[nkappa,t] += gamma * x[t]^(c+1)
end
end
end
end
end
if k == 0
if nkappa > 1
jarray[nkappa,t] += jarray[nkappa, t-1]
end
else
jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
end
return jarray
end # end jack -------------------------------------------------------------
function summation(
i::I,
z::T,
j::I,
kappa::Vector{I},
jarray::Array{T,2}
) where {I<:Integer, R<:Real, T<:Union{R,Complex{R}}}
function go(kappai::I, zp::T, s::T)
if i == n || i == 0 && kappai > j ||
i > 0 && kappai > min(last(kappa), j)
return s
end
kappap = vcat(kappa, [kappai])
nkappa = _N(dico, kappap) - 1
zpp = zp * _T(alpha, a, b, kappap)
if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
jarray[nkappa,1] =
x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
end
for t in 2:n
jarray =
jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
end
sp = s + zpp * jarray[nkappa,n]
if j > kappai && i <= n
spp = summation(i+1, zpp, j-kappai, kappap, jarray)
go(kappai+1, zpp, sp + spp)
else
go(kappai+1, zpp, sp)
end
end # end go -----------------------------------------------------------
go(I(1), T(z), T(0))
end # end summation --------------------------------------------------------
the_one = eltype(a) <: Complex || eltype(b) <: Complex ||
eltype(x) <: Complex ? Complex(R(1)) : R(1)
(dico, Pmn) = DictParts(m, n)
J = zeros(typeof(the_one), Pmn, n)
J[1,:] = cumsum(x)
return the_one + summation(0, the_one, m, Int64[], J)
end
# univariate -------------------------------------------------------------------
# both work
function hypergeom(
m::Integer,
a::Vector{<:Union{R,Complex{R}}},
b::Vector{<:Union{R,Complex{R}}},
x::Union{R,Complex{R}}
)::Union{R,Complex{R}} where {R<:Real}
return hypergeomI(m, R(2), a, b, 1, x)
end
function hhypergeom(
m::Integer,
a::Vector{<:Union{R,T}},
b::Vector{<:Union{R,T}},
x::Union{R,T}
) where {R<:Real,T<:Complex{R}}
return hypergeomI(m, R(2), a, b, 1, x)
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment