{{ message }}

Instantly share code, notes, and snippets.

# stla/hypergeomPQ.jl

Last active Oct 24, 2020
Hypergeometric function of a matrix argument in Julia (Koev & Edelman's algorithm)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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 file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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