--------------------------------------------------------------------------------
------------------------ General Purpose Functions -----------------------------
--------------------------------------------------------------------------------




-- resize a matrix, the input matrix is centered in the new one
-- If output bigger than input, new cell padded with value indicated in buffer
local function resize_matrix(matrix, outputDim1, outputDim2, buffer)

    local myOutput  = make_matrix(buffer, outputDim1, outputDim2)

    -- Get dimensions
    local inputDim1 = #matrix
    local inputDim2 = #matrix[1]
    
    --Check Dim : dimensions must remain odd or even after transformation
    if (inputDim1 + outputDim1) % 2 == 1 or (inputDim2 + outputDim2) % 2 == 1
    then
        print("resize_matrix, odd/event mixture in adimension")
        os.exit()  
    end
    
    --Differential of dimentions
    local deltaDim1 = outputDim1 - inputDim1
    local deltaDim2 = outputDim2 - inputDim2
    
    -- Iterate on output & copy input if in the right range
    for out1 = 1, outputDim1, 1
    do
        local in1 = out1 - (deltaDim1/2)
        for out2 = 1, outputDim2, 1
        do
            local in2 = out2 - (deltaDim2/2)
            if (in1 > 0 and in1 <= inputDim1) and (in2 > 0 and in2 <= inputDim2)
            then 
                myOutput[out1][out2] = matrix[in1][in2]
            end
        end
    end
    
    return myOutput            
end


-- Subsets a matrix from a bigger one and returns result as vector (add safeties?), 
local function matrixSubsetVect(startDim1, startDim2, lenDim1, lenDim2, matrix)
    local tempVect = {}
    local count = 1
    
    if((startDim1 + (lenDim1 - 1) > #matrix) or (startDim2 + (lenDim2 - 1) > #matrix[1]))
    then
        os.exit()
    end
    
    for i = startDim1, startDim1 + (lenDim1 - 1), 1 
    do
        for j = startDim2, startDim2 + (lenDim2 - 1), 1 
        do
            tempVect[count] = matrix[i][j]
            count = count + 1
        end    
    end
    return tempVect
end


-- Extracts a submatrix from a matrix, as vector
-- The center of the matrix is aligned with indicated coord of submatrix
local function subsetFromCenter(matrix, subDim1, subDim2, centerAlign1, centerAlign2)
    -- Dimensions of input matrix
    local dimMat1 = #matrix
    local dimMat2 = #matrix[1]
    
    -- Center of input matrix
    local centerMat1 = math.ceil(dimMat1/2)
    local centerMat2 = math.ceil(dimMat2/2)
    
    -- Beginning positions of the subsample
    local startDim1 = centerMat1 - (centerAlign1 - 1)
    local startDim2 = centerMat2 - (centerAlign2 - 1)
    
    --Create subset matrix
    return matrixSubsetVect(startDim1, startDim2, subDim1, subDim2, matrix)
end  



-- Rotates a square matrix on spot
local function rotateMatrix(matrix)
    if #matrix ~= #matrix[1]
    then
        print("rotateMatrix: matrix not squared")    
        os.exit()
    end
    
    for i = 1, #matrix, 1
    do
        for j = i, #matrix, 1
        do
            matrix[i][j], matrix[j][i] = matrix[j][i], matrix[i][j] 
        end
    end

end


--------------------------------------------------------------------------------
---------------------  Predefined flux pattern functions -----------------------
--------------------------------------------------------------------------------
-- Island model, migration rate gives the % of elements that spread outside of the cell
-- Also takes the dimensions of the landscape as defined in the config file
local function make_island(ncols, nlines, migrationRate)
    --print("makeisland", migrationRate, ncols, nlines)
    -- Nlines and ncols of the flux matrix
    local fPatCols = (ncols*2)-1
    local fPatLines = (nlines*2)-1
    -- Size of the pattern grid (twice the landscape minus one)
    local size = fPatCols * fPatLines
    -- Value of each non center cell
    --local extValue = migrationRate/(size-1)
    local extValue = migrationRate/(ncols*nlines-1)
    --Create matrix filled with extValue
    local fluxPattern = make_matrix(extValue, fPatLines, fPatCols) 
    --Change center value (= 1 - migrationRate)
    fluxPattern[math.ceil(fPatLines/2)][math.ceil(fPatCols/2)] = 1-migrationRate
    return fluxPattern
end


--Stepping stone, pattern is of fixed size (3x3)
local function make_step_stone(migrationRate)
    -- Value of direct neighbours
    local neighBourValue  = migrationRate/4
    -- Create matrix filled with 0
    local fluxPattern = make_matrix(0.0, 3, 3)
    -- Fill neighbours
    fluxPattern[1][2] = neighBourValue
    fluxPattern[2][1] = neighBourValue
    fluxPattern[2][3] = neighBourValue
    fluxPattern[3][2] = neighBourValue
    -- Fill center 
    fluxPattern[2][2] = 1 - migrationRate
    return fluxPattern
end



-- Generate flux matrix from a function mane:
--      Redirects to the right function type by name
-- Gather parameters and check for type / existence before calling them
-- Loads for both sexes
local function make_flux_pattern(funcName, ms, mp, mapWidth, mapLength, propagule)
    local fluxPattern = {}

    -- Island model takes dimensions and migration rates of both gametes types
    if funcName == "island" then
        -- Function specific parameters
        local nCols = assert_number(mapWidth)
        local nLines = assert_number(mapLength)
        local m = 0
        -- Gather propagule-specific parameter (mp & sp)
        if propagule == "seed" then m = assert_number(ms)
        elseif  propagule == "pollen" then m = assert_number(mp)
        else print("ERROR, load_flux_function, invalid propagule "..propagule); os.exit()
        end
        -- Flux pattern generation
        fluxPattern = make_island(nCols, nLines, m)

    -- Stepping stone only takes the migration rates
    elseif funcName == "steppingStone" then
        -- Gather propagule-specific parameter (mp & ms)
        if propagule == "seed" then m = assert_number(ms)
        elseif  propagule == "pollen" then m = assert_number(mp)
        else print("ERROR, load_flux_function, invalid propagule "..propagule); os.exit()
        end
        -- Flux pattern generation
        fluxPattern = make_step_stone(m)
        
    ----------------------------------------
    -- Add other functions here if needed...
    ----------------------------------------
    
    -- Error if not a function nameknown namevector
    else 
        print("Unknown flux function:" .. funcName) 
        os.exit()
    end
    return fluxPattern
end





--------------------------------------------------------------------------------
-------------------  Converting flux Pattern to flux Matrix --------------------
--------------------------------------------------------------------------------

-- Function to create the flux matrix, returns new table
-- The first dim represents donnors, the second dim the recivers
local function makeFluxMat(pattern, nlongi, nlati)
    -- Init of an empty flux matrix (right dimensions)
    local dim = nlati*nlongi
    local fluxMatrix = make_matrix(0, dim, dim)
    
    -- pattern must be adjusted to the size of (2*dim)-1
    -- Passing with 0.0
    local rPattern = resize_matrix(pattern, 2 * nlati - 1, 2 * nlongi - 1, 0.0)

    -- Iterating on donnors (all pop/cell)
    for donnor = 1, #fluxMatrix, 1 
    do
        -- Calc coordinates of the donnor on the map (we could iterate on longi and lati instead)
        local lati = math.ceil(donnor / nlongi)
        local longi = (donnor-1) % nlongi + 1  --(donnor-1) allows to start at zero, then 1 is added

        -- Extracts relevant part of the rPattern matrix as a vector
        local tempVect = subsetFromCenter(rPattern, nlati, nlongi, lati, longi)
        -- Each line is nomalized, but we also normalize columns in the end (more important)
        tempVect = stats.normalize(tempVect)
        -- Add the vector to the flux matrix
        fluxMatrix[donnor] = tempVect
    end
   
    -- Rotate flux matrix
    rotateMatrix(fluxMatrix)
    
    -- Normalize recivers (now on lines) TODO:check this
    for i = 1, #fluxMatrix, 1
    do
       fluxMatrix[i] = stats.normalize(fluxMatrix[i])
    end
    
    -- A new matrix is returned
    return fluxMatrix
end

--------------------------------------------------------------------------------
-------------------------------  Module interface ------------------------------
--------------------------------------------------------------------------------

local module = {}
module.makeFluxMat = makeFluxMat
module.make_flux_pattern = make_flux_pattern

return module
