-- Statistics and math module

-- Sum of a vector 
local function sum_vector(vector)
    local sum = 0
    for i = 1, #vector do
        sum = sum + vector[i]
    end
    return sum
end

-- Sum of a matrix 
local function sum_matrix(matrix)
    local sum = 0
    for i = 1, #matrix do
        sum = sum + sum_vector(matrix[i])
    end
    return sum
end

-- Mean of a vector 
local function mean_vector(vector)
    return sum_vector(vector)/#vector
end

-- Mean of a matrix 
local function mean_matrix(matrix)
    return sum_matrix(matrix)/(#matrix*#matrix[1])
end

-- Applies a function to each value of a vector, returns a new vector
-- func must take data as the fist argument
-- Additional arguments of func can be passed at the end (...)
local function apply_vector(vector, func, ...)
    local new = {}
    for i = 1, #vector do
        new[i] = func(vector[i], ...)
    end
    return new
end

-- Applies a function to each value of a matrix, returns a new matrix
local function apply_matrix(matrix, func, ...)
    local new = {}
    for i = 1, #matrix do
        new[i] = apply_vector(matrix[i], func, ...)
    end
    return new
end

-- TODO: add in place variants of apply

-- Variance of a vector
local function variance_vector(vector)
    local mean = mean_vector(vector)
    local func = function(x, y) return (y - x)^2 end
    return mean_vector(apply_vector(vector, func, mean))
end

-- Variance of a matrix
local function variance_matrix(matrix)
    local mean = mean_matrix(matrix)
    local func = function(x, y) return (y - x)^2 end
    return mean_matrix(apply_matrix(matrix, func, mean))
end

-- Normalize a vector
local function normalize(vector)
    local sum = sum_vector(vector)
    local func = function(x, y) return x/y end
    return apply_vector(vector, func, sum)
end

----------------------------- Random functions ---------------------------------
--TODO : move all this in CPP

function Random() 
    --/* uniforme entre 0 compris et 1 non compris */
    -- Générateur de nombre pseudo-aléatoires. Utilise les seeds S1 et S2 et les modifie à chaque appel.
    --
    -- Reference : L'ECUYER, 1988 in Goudet, J., 1993, The genetics of geographically
    -- structured populations, School of Biological Sciences University College of North
    -- Wales.    
    local i = 0
    local r = 0
    local z = 0
    r = S1 / 53668
    r = math.floor(r)
    S1 = ( 40014 * (S1 - (r * 53668) ) ) - (r * 12211)
    if (S1 < 0) then S1 = S1 + 2147483563 end
    r = S2 / 52774
    r = math.floor(r)
    S2 = ( 40692 * (S2 - (r * 52774) ) ) - (r * 3791)
    if (S2 < 0) then S2 = S2 + 2147483399 end
    z = S1 - S2
    if (z < 1) then z = z + 2147483562 end
    return (z * 4.656613e-10)
end





--Box-Muller transformation
flag = 0
r2 = 0
function gasdev()   
        local fac = 0
        local rsq = 0
        local v1 = 0
        local v2 = 0
        if(flag == 0) then 
                repeat
                    v1 = ( 2.0 * Random() ) - 1.0
                    v2 = ( 2.0 * Random() ) - 1.0
                    rsq = math.pow(v1,2) + math.pow(v2,2)
                until ( not ((rsq >= 1.0) or (rsq == 0.0)) )
                fac = math.sqrt( (-2.0 * math.log(rsq)) / rsq )
                r2 = v1 * fac
                flag = 1
                return (v2 * fac)
        else
                flag = 0
                return r2
        end
end

-- Tirage de k éléments dans [a, b] sans remise
function get_k(k, a, b)
    local vect  = {}
    local limit = (b - a + 1) -- +1 to include sup' limit
    local selected = limit
    assert(k >= limit, "get_k: number of elements k("..k..") < limit ("..limit..")")
    -- Creates values to pick from
    for i = 1, limit do
        vect[i] = a + i - 1 -- (-1) to start at a and not a + 1
    end
    -- Pick k values and put them at the end of the table
    for i = 1, k do
        indice = math.floor(Random() * selected + 1)
        vect[indice], vect[selected] =  vect[selected], vect[indice]
        selected = selected - 1
        
        range = b - i - a
    end
    -- Copy the selected values in a table and return it
    local values = {}
    for i = selected, limit do
        table.insert(values, vect[i])
    end
    return values
end

--------------------------------------------------------------------------------
local stats = {}
stats.sum_vector = sum_vector
stats.sum_matrix = sum_matrix
stats.mean_vector = mean_vector
stats.mean_matrix = mean_matrix
stats.apply_vector = apply_vector
stats.apply_matrix = apply_matrix 
stats.variance_vector = variance_vector
stats.variance_matrix = variance_matrix
stats.normalize = normalize
return stats
