# \[Help needed\] Translate 10 lines from Maple to Julia

**URL:** <https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650>\
**Category:** New to Julia\
**Created:** [May 6, 2021, 1:41pm UTC](https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650 "2021-05-06T13:41:49Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 6, 2021, 1:41pm UTC](https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650/1 "2021-05-06T13:41:49Z")

</div>

Hey,

I have a somewhat weird request: I do not speak Maple at all, but i need to translate the 10 lines of code from [this document](https://iris.unito.it/retrieve/handle/2318/1561388/399488/FASTpolykaysMAPLE.pdf), to be able to use the function `ks()` that they created on page 2. Notes on underlying other functions are on page 3 and 4, the rest of the document gives generalisations that I do not need.

This function computes a polynomial in `S_1,...,S_m` that corresponds to the expression of the k-statistic (unbiaised estimator of the kth cumulant):

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/1/9164b6e6773c9129abcbb25b3306a6c18dfb444a.png)

Details on the corresponding theory, at least for the first few cases, is available there: [k-Statistic -- from Wolfram MathWorld](https://mathworld.wolfram.com/k-Statistic.html)

By reading the Maple code, do you think it’s going to be complicated ?

Edit: I found the partition and stirling2 functions in the Combinatorics parckage. This is what i have for the moment:

```julia
import Combinatorics

# The partition Maple function is replaed by : 
# Combinatorics.integer_partitions(k)

# stirling2 is replaced by 
# Combinatorics.stirlings2(k,n)

# Let's try maketab: 
function makeTab_ss(N)
    

end

```

I still need to undertand wtf is doing the makeTab\_ss function 😉

---

<div class="post-metadata">

**Author:** ![KZiemian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kziemian/32/9020_2.png) [@KZiemian](https://discourse.julialang.org/u/KZiemian)\
**Post date:** [May 6, 2021, 1:57pm UTC](https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650/2 "2021-05-06T13:57:06Z")

</div>

> [@lrnv](#):
>
> By reading the Maple code, do you think it’s going to be complicated?

If you need just naive Julia implementation (not-as-fast-as-it-could-be) only problem is to decipher Maple code. I quite rusty in things like Maple, but, if my time allows, I can try to help you with that. You can send me message if you want so more help.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 6, 2021, 1:58pm UTC](https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650/3 "2021-05-06T13:58:36Z")

</div>

> [@KZiemian](#):
>
> If you need just naive Julia implementation (not-as-fast-as-it-could-be) only problem is to decipher Maple code. I quite rusty in things like Maple, but, if my time allows, I can try to help you with that. You can send me message if you want so more help.

I need something that works and i do not need it to be fast for the moment, as it would be a one time computation that I then store somewhere.

If you can help me understand the Maple code it would be great !

---

<div class="post-metadata">

**Author:** ![KZiemian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kziemian/32/9020_2.png) [@KZiemian](https://discourse.julialang.org/u/KZiemian)\
**Post date:** [May 6, 2021, 2:09pm UTC](https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650/4 "2021-05-06T14:09:01Z")

</div>

Ok, how much time we have? I’m now quite busy, but I probably can sit down in few hours. Here is now 16:08, around 19 I should have some time.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 6, 2021, 2:10pm UTC](https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650/5 "2021-05-06T14:10:00Z")

</div>

Well, i’ll be happy to show you what i have by 19 ! Thanks a lot.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 6, 2021, 3:51pm UTC](https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650/6 "2021-05-06T15:51:22Z")

</div>

Ok i do have the following code:

```julia
# The goal of this file is to compute the k_statistics polynomials ad given by 
# https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650
# and 
# https://iris.unito.it/retrieve/handle/2318/1561388/399488/FASTpolykaysMAPLE.pdf
# also see 
# Di Nardo E., G. Guarino, D. Senato (2007), A new method for fast computing unbiased estimators of cumulants. In press Statistics and Computing. http://www.springer.com/statistics/computational/journal/11222 (download from http://www.unibas.it/utenti/dinardo/lavori.html)

import Combinatorics
import Polynomials
import StatsBase
import Random
import TypedPolynomials

function makeTab_ss(N::T) where T<:Integer
    part = Combinatorics.integer_partitions(N) 
    dict = Dict{Int64, Int64}()
    Nf = factorial(N)
    result = Array{Tuple{Rational{Int64}, Vector{Int64}}}(undef,0)
    for y in part
        denom = T(1)
        dict = StatsBase.countmap(y)
        for (val,numoccur) in dict
            denom *= factorial(val)^numoccur * factorial(numoccur)
        end
        push!(result,(Nf//denom,y))
    end
    return result
end

function Stirling_triangle(N::T) where T<:Integer
    arr = zeros(T,N,N)
    for i in 1:N
        for j in 1:i
            arr[i,j] = Combinatorics.stirlings2(i,j)
        end
    end
    return arr
end

function makeK_s(N::T) where T<:Integer
    result = Array{Polynomials.Polynomial{Float64}}(undef,0)
    Stirling2 = Stirling_triangle(N)
    for i in 1:N
        rez = zeros(i+1)
        for j in 1:i
            rez[j+1] = Combinatorics.stirlings2(i,j) * (-1)^(j-1) * factorial(j-1)
        end
        push!(result,Polynomials.Polynomial(rez))
    end
    return result
end

function expand_eval(v,u)
    result = Array{eltype(u)}(undef,0)
    for (number,array) in v
        rez = number * prod(u[i] for i in array)::eltype(u)
        push!(result,rez)
    end
    return result
end

function fd(j,k)
    Polynomials.fromroots([i+k for i in 0:(j-k)])
end
# checks from the note in the source: 
@assert fd(3,0) == Polynomials.Polynomial([0,-6,11,-6,1])
@assert fd(3,1) == Polynomials.Polynomial([-6,11,-6,1])
@assert fd(3,2) == Polynomials.Polynomial([6,-5,1])

function build_u2(N)
    result = []
    for i in 1:N
        rez = (-1)^(i-1)*factorial(i-1)*fd(N-1,i)
        push!(result,rez)
    end
    return result
end

function compute_k_representation(n::T,k_idx) where T

    v = makeTab_ss(k_idx)
    u = makeK_s(k_idx)
    v2 = expand_eval(v,u)
    poly = zero(eltype(u))
    u2 = build_u2(k_idx)
    common_denominator = prod(n - x for x in 0:(k_idx-1))
    factors = []
    for j in 1:size(v2,1)
        poly = sum((v2[j].coeffs[i+1]*u2[i]) for i in 1:(length(v2[j].coeffs)-1))
        numerator = poly(n)
        push!(factors,numerator/common_denominator)
    end
    which_S = [v[i][2] for i in 1:length(v)]
    return (factors,which_S)
end

function compute_k(N,S)
    m = length(S)
    result = []
    for k in 1:m
        repr = compute_k_representation(N,k)
        rez = sum(repr[1][i] * prod([S[j] for j in repr[2][i]]) for i in 1:size(repr[1],1))
        push!(result,rez)
    end
    return result
end

# Simple experiment: 
# Vary the number of observations N, and the mean and variance of the normal random vector. 
# The first cumulant should be the mean, 
# the second should be the variance 
# and all others should be 0 since it's a normal. 
# the approximation should get better and better as you increase the number of pointS. 
mean = 0
var = 1
m=7
N_obs = 100000
Random.seed!(123)
data = mean .+ sqrt(var) .* randn(N_obs)
S_obs = [sum(data .^ i) for i in 1:m]
cumulants_direct = compute_k(big.(N_obs),big.(S_obs)) # directly without the poynomials.

# As a simple check, we can verify that the polynomials holds : 
TypedPolynomials.@polyvar N
TypedPolynomials.@polyvar S[1:7] ##########" put m here in a fixed way, as the infrastructure does not allow for variables.
result = compute_k(N,S);
@assert result[1] == S[1]/N
@assert result[2] == (N * S[2] - S[1]^2)/(N*(N-1))
@assert result[3] == (2S[1]^3 - 3N*S[1]*S[2] +N^2 * S[3]) / (N*(N-1)*(N-2))
@assert result[4] == (-6*S[1]^4 + 12 * N * S[1]^2 * S[2] - 3 * N * (N-1) * S[2]^2 - 4 * N * (N+1) * S[1] * S[3] + N^2*(N+1)*S[4])/(N*(N-1)*(N-2)*(N-3))
cumulants_poly = [result[i](N=>big.(N_obs),S=>big.(S_obs)) for i in 1:m]

display([cumulants_direct cumulants_poly])

```

The code is clearly suboptimal, however with the introduction of TypedPolynomials, i was able to check that it is correct.

Now the question is :

- Is it possible to compute once and for all the polynomials, for any N and S, and store them to disk somewhere in a package that i can load easily ?
- Is it possible to do the same thing (pre-compute the polynomials) for a given N, on the fly ?

The use case is that i might compute the value of these polynomials a lot of time, for varying S and for fixed N.

---

<div class="post-metadata">

**Author:** ![cescalantec](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cescalantec/32/18795_2.png) [@cescalantec](https://discourse.julialang.org/u/cescalantec)\
**Post date:** [October 12, 2022, 6:31pm UTC](https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650/7 "2022-10-12T18:31:48Z")

</div>

Puede servir darle una mirada a

> **[Julia - Maple Help](https://www.maplesoft.com/support/help/maple/view.aspx?path=CodeGeneration%2FJulia)**
>
> CodeGeneration Julia translate Maple code to Julia code Calling Sequence Parameters Description Examples Compatibility Calling Sequence Julia( x , cgopts ) Parameters x - expression, list, rtable, procedure, or module cgopts - (optional) one or...