#! /usr/bin/Rscript
source("fetching_functions.r")
source("metapop_functions.r")

args =  commandArgs(TRUE)
argSpecies = args[1]
print(argSpecies)
argStep  = args[2]
print(argStep)
argInput = args[3]
print(argInput)

argSelected=as.integer(strsplit(args[4],",")[[1]])
print(argSelected)


argTraits = args[c(-1,-2,-3,-4)]
print(argTraits)
traitName = substring(argTraits, 1, nchar(argTraits)-6)
argTraitsAdd = paste(traitName, "-add", sep='')
print(argTraitsAdd)
argTraitsDz = paste(traitName, "-dz", sep='')

argTraitsAddByLoc = paste(substring(argTraits, 1, nchar(argTraits)-6), "-add-bylocus", sep='')
print(argTraitsAddByLoc)
#Set correct directory
setwd(argInput)

### GET PARAMETER ###
dimMap = getMapDim(argInput)#get value automatically from the MTP/user_input/conf.lua
mapWidth = dimMap[1]
mapHeight = dimMap[2]

print(mapWidth)
print(mapHeight)



##################################################
# Construct a general list of files of interset  #
#  = .csv, 3 fields separated by "_"             #
##################################################
#Get all files
filesList = list.files()

fileAndsub = findCsv(filesList, FALSE)
s3 = fileAndsub[[1]]
filesList = fileAndsub[[2]]
sublists = fileAndsub[[3]]

#print(paste("fL", filesList))

###############################################################################
# Checks

#Checking consistency of arguments
checkConsistencyPheno(filesList, argSpecies, argTraits, argStep)

###############################################################################
###############################################################################
#########                         Analysis                            #########
###############################################################################
###############################################################################


traits_pheno_Qst = sapply(argTraits, function(x) calc_qst_generic(argSpecies, x, argStep, FALSE))
traits_pheno_Means = lapply(argTraits, function(x) pop_mean_generic(argSpecies, x, argStep))

traits_dz_Means = pop_mean_generic(argSpecies, argTraitsDz, argStep)

for(nTrait in 1:length(traits_pheno_Means))
{
    print(paste("Pop pheno: ", argTraits[[nTrait]]))
    print(traits_pheno_Means[[nTrait]])
    print("Mean pheno :")
    print(traits_pheno_Qst[[nTrait]])
}

paste("[F] :",traitName,"-> DZ :", mean(traits_dz_Means))

traits_Qst = sapply(argTraitsAdd, function(x) calc_qst_generic(argSpecies, x, argStep, TRUE))
traits_Means = lapply(argTraitsAdd, function(x) pop_mean_generic(argSpecies, x, argStep))

for(nTrait in 1:length(traits_Means))
{
    print(paste("Populations means for trait: ", argTraits[[nTrait]]))
    print(traits_Means[[nTrait]])
    print("Qst value:")
    print(traits_Qst[[nTrait]])
}

paste("[F] :",traitName,"-> Mean Qst :", mean(traits_Qst))

#########################################
#       PARSING OF ADDITIVE VALUE       #
#########################################

valadd2 = function(species, trait, step)
{
    targetFile = filesList[[argSpecies]][[trait]][[argStep]]
    if (file.info(targetFile)$size > 0)
    {
        vp = split_vals_generic(targetFile)
        values = vp[[1]]
        pops = vp[[2]]

        Vmean = sapply(values, mean)
        Vvar = sapply(values, var)
        pos = 1
        print(paste("T:", trait))
        if (grepl("add", trait))
        {
                for (i in pops)
                {
                    print(paste("[F] :", traitName,"1-ADDMEAN", traitName, "_", i, " : ", Vmean[pos]))
                    print(paste("[F] :", traitName,"2-ADDVAR", traitName, "_", i, " : ", Vvar[pos]))
                    pos <- pos + 1
                }
        }
        else
        {
                for (i in pops)
                {
                    print(paste("[F] :", traitName,"1-PHENOMEAN", traitName, "_", i, " : ", Vmean[pos]))
                    print(paste("[F] :", traitName,"2-PHENOVAR", traitName, "_", i, " : ", Vvar[pos]))
                    pos <- pos + 1
                }
        }

        # FILL EMPTY POSITIONS
        VmeanFilled = rep(NaN, mapWidth*mapHeight)

        pos = 1
        while (pos <= length(pops))
        {
            VmeanFilled[pops[pos]+1] = Vmean[pos]
            pos <- pos + 1
        }
        return(VmeanFilled)
    }
}


madd = t(matrix(valadd2(argSpecies, argTraitsAdd, argStep), ncol = mapHeight, nrow = mapWidth))
# WORKS ONLY FOR ONE TRAIT
print(madd)


#########################################
#    PARSING OF ADDITIVE VALUE BY LOC   #
#########################################

valaddByloc = function(species, trait, step)
{
    targetFile = filesList[[argSpecies]][[trait]][[argStep]]
    print(paste("TF:", targetFile))
    if(!length(targetFile)){
	print("NO ADDBYLOC FILE, no theta_b nor theta_w computed")
	return(0)
    }
    if (file.info(targetFile)$size > 0)
    {
        return(valByloc_generic(targetFile, argSelected, traitName, seq(length(argSelected))))
    }
}

#traits_theta = sapply(argTraitsAddByLoc, function(x) valaddByloc(argSpecies, x, argStep))
traits_theta = valaddByloc(argSpecies, argTraitsAddByLoc, argStep)
print(traits_theta)

# get additive val vector
meanAdd = apply(madd, 1, mean, na.rm = TRUE)
#mean N
print(paste("[F] :",traitName,"-> Addval Mean North :", meanAdd[1]))
#mean S
print(paste("[F] :",traitName,"-> Addval Mean South :", tail(meanAdd,1)))

varAdd = apply(madd, 1, var, na.rm = TRUE)
#std N
print(paste("[F] :",traitName,"-> Addval Var North :", varAdd[1]))
#std S
print(paste("[F] :",traitName,"-> Addval Var South :", tail(varAdd,1)))

#REAL STD COMPUTATION
valvar = function(species, trait, step)
{
    targetFile = filesList[[argSpecies]][[trait]][[argStep]]
    if (file.info(targetFile)$size > 0)
    {
        vp = split_vals_generic(targetFile)
        values = vp[[1]]
        pops = vp[[2]]
        Va = sapply(values, var)

        # FILL EMPTY CELLS
        pos = 1
        VaFilled = rep(NaN, mapWidth*mapHeight)
        while (pos <= length(pops))
        {
            VaFilled[pops[pos]+1] = Va[pos]
            pos <- pos + 1
        }
        return(VaFilled)
    }
}
mvar = t(matrix(valvar(argSpecies, argTraitsAdd, argStep), ncol = mapHeight))
print(mvar)
print(paste("[F] :",traitName,"-> Mean(Addval Var) :", mean(mvar, na.rm = TRUE)))

#std
print(paste("[F] :",traitName,"-> Addval Var Glob :", var(as.vector(madd), na.rm = TRUE)))
print(paste("[F] :",traitName,"-> Addval Mean Glob :", mean(as.vector(madd), na.rm = TRUE)))
