##! /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)]
argTraits = args[c(-1,-2,-3)]
print(argTraits)

traitName = substring(argTraits, 1, nchar(argTraits)-6)

argTraitsAdd = paste(traitName, "-add", sep='')
print(argTraitsAdd)

argTraitsEps = paste(traitName, "-epsi", sep='')
print(argTraitsEps)

#Set correct directory
setwd(argInput)
    
##################################################    
# 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_Qst = sapply(argTraitsAdd, function(x) calc_qst_generic(argSpecies, x, argStep, FALSE))
traits_Means = lapply(argTraitsAdd, function(x) pop_mean_generic(argSpecies, x, argStep))


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


valadd3 = function(species, tr, step, wantmean, wantP)
{
    traitA= paste(tr, "-add", sep='')
    rawvaluesA = read.table(filesList[[argSpecies]][[traitA]][[argStep]], sep = " ", header = F)
    valuesA = rawvaluesA[order(rawvaluesA[,2]),]

    #Get pop numbers index
    popIndex = valuesA[,1]
    #Get pops ids
    pops = unique(popIndex)
    #Make a mask to select inds for each pops
    splitter = (lapply(pops, '==', popIndex))
    #Isolate values
    valuesA = valuesA[,3]
    #Separate pops
    values = lapply(splitter, function(x) valuesA[x])
    if(wantP){
	    traitE = paste(tr, "-epsi", sep='')
	    rawvaluesE = read.table(filesList[[argSpecies]][[traitE]][[argStep]], sep = " ", header = F)
	    valuesE = rawvaluesE[order(rawvaluesE[,2]),]
	    valuesE = valuesE[,3]
            values = lapply(splitter, function(x) valuesA[x] + valuesE[x])
    }   
    if(wantmean)
    {
        Vb = sapply(values, mean)
        return(Vb)
    }
    else
    {
        return(values)
    }
}


madd = sapply(traitName, function(x) valadd3(argSpecies, x, argStep, FALSE, FALSE))
print(madd)
print(madd[1,1])

print(dim(madd))

# I.a R intra for 2 carac
intraR = sapply(seq(dim(madd)[1]), function(x) rat_corr_generic(madd[x,1][[1]], madd[x,2][[1]]) )
print(intraR)
print(sum(intraR))


# I.b R inter for 2 carac
maddMean = sapply(traitName, function(x) valadd3(argSpecies, x, argStep, TRUE, FALSE))
#interR = cov(maddMean[,1], maddMean[,2]) / (var(maddMean[,1])*var(maddMean[,2]))
interR = cov(maddMean[,1], maddMean[,2]) / (sd(maddMean[,1])*sd(maddMean[,2]))
print(interR)

# II.[G] for all possible combinations
nt = length(argTraitsAdd)
print(nt)
Ginter = matrix(nrow=nt, ncol=nt)
Gintra = matrix(nrow=nt, ncol=nt)


for(i in seq(nt))
{
    for(j in seq(nt))
    {
        Ginter[i,j] = cov(maddMean[,i], maddMean[,j]) / (sd(maddMean[,i])*sd(maddMean[,j]))
        matrGintra = sapply(seq(dim(madd)[1]), function(x) rat_corr_generic(madd[x,i][[1]], madd[x,j][[1]]) )
        print(paste(i, j))
        print(matrGintra)
        Gintra[i,j] = 1/(dim(madd)[1])*sum(matrGintra)
    }
}

for(i in seq(nt))
{
    for(j in seq(i,nt))
    {
        print(paste("[F] : Ginter[",i,",",j,"] :", Ginter[i,j]))
        print(paste("[F] : Gintra[",i,",",j,"] :", Gintra[i,j]))
    }
}

# II.b [P] for all possible combinations
#formula : P1 = A1 + epsilon
# Should the P computation be redone or use the Tn-pheno file?

# use files argTraitsEps

maddP = sapply(traitName, function(x) valadd3(argSpecies, x, argStep, FALSE, TRUE))
print(maddP)
print(maddP[1,1])

print(dim(maddP))

intraRP = sapply(seq(dim(maddP)[1]), function(x) rat_corr_generic(maddP[x,1][[1]], maddP[x,2][[1]]) )
print(intraRP)
print(sum(intraRP))

maddMeanP = sapply(traitName, function(x) valadd3(argSpecies, x, argStep, TRUE, TRUE))
#interR = cov(maddMean[,1], maddMean[,2]) / (var(maddMean[,1])*var(maddMean[,2]))
interRP = cov(maddMeanP[,1], maddMeanP[,2]) / (sd(maddMeanP[,1])*sd(maddMeanP[,2]))
print(interRP)

nt = length(argTraitsAdd)
print(nt)
Pinter = matrix(nrow=nt, ncol=nt)
Pintra = matrix(nrow=nt, ncol=nt)

for(i in seq(nt))
{
    for(j in seq(nt))
    {
        Pinter[i,j] = cov(maddMeanP[,i], maddMeanP[,j]) / (sd(maddMeanP[,i])*sd(maddMeanP[,j]))
        Pintra[i,j] = 1/dim(madd)[1]*sum(sapply(seq(dim(maddP)[1]), function(x) rat_corr_generic(maddP[x,i][[1]], maddP[x,j][[1]]) ))
    }
}

for(i in seq(nt))
{
    for(j in seq(i,nt))
    {
        print(paste("[F] : Pinter[",i,",",j,"] :", Pinter[i,j]))
        print(paste("[F] : Pintra[",i,",",j,"] :", Pintra[i,j]))
    }
}
