#! /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)


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

#Set correct directory
#setwd("../raw_outputs/values/")
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]]

### SPECIES (after findCsv processing) ###

#Split names list according to types
#speList = lapply(sublists, function(x) s3[x])
#names(speList) = speNames


################### TYPES #################
# Same, but for types inside each species #
###########################################
#filesList = findTraitCsv(speList)
#print(paste("fL", filesList))

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

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


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

getNlocByLt = function(fileRef)
{
  # get name of loci for each column based on first line of the addbylocus file
  header = readLines(fileRef, n=1)
  pos = unlist(gregexpr(pattern =':',header))
  SPLT = strsplit(unlist(strsplit(substr(header, pos+2, nchar(header)), ";")), " ")
  lengthByLt = lapply(SPLT, length)
  totL = 0
  for (i in seq(length(lengthByLt)))
  {
    totL = totL + lengthByLt[[i]]
  }
  return(list(totL, lengthByLt))
}

nameltype = c("a", "b", "c")

calc_qst_lt = function(species, trait, step)
{
    ltypeQst = c()
    currentNumLt = 1
    currentPosSlice = 1
    for (ltype in nameltype)
    {
        ltypeTrait <- paste(trait, "-", ltype, sep="")
        targetFiles = filesList[[argSpecies]][[ltypeTrait]]
        if (! (is.null(targetFiles)))
        {
            vp = split_vals_generic(targetFiles[[argStep]])
            values = vp[[1]]
            pops = vp[[2]]

            Vmean = sapply(values, mean)
            Vvar = sapply(values, var)
            pos = 1
            for (i in pops)
            {
                print(paste("[F] :", traitName, ltype, "3-ADDMEAN", traitName, ltype, "_", i, " : ", Vmean[pos]))
                print(paste("[F] :", traitName, ltype, "4-ADDVAR", traitName, ltype, "_", i, " : ", Vvar[pos]))
                pos <- pos + 1
            }

            Vb = var(sapply(values, mean))
            Vw = mean(sapply(values, var, na.rm=TRUE))
            print(paste("[F] :", traitName,"->",paste("Vb",ltype, sep="_"),":", Vb))
            print(paste("[F] :", traitName,"->",paste("Vw",ltype, sep="_"),":", Vw))
            if (is.na(Vw))
            {
                print("[Warning] Vw is NaN (maybe because one pop only has one individual)")
                ltypeQst = c(ltypeQst, "")
            }
            else
            {
                Qst = Vb / (Vb + 2 * Vw)
                ltypeQst = c(ltypeQst, Qst)

                # PRINT THETAW AND THETAB FOR GIVEN LTYPE
                addbylocFile = filesList[[argSpecies]][[paste(traitName, "-add-bylocus", sep='')]][[argStep]]

                lociinfo = getNlocByLt(addbylocFile)
                nloc = lociinfo[[1]]
                selectedLoci = lociinfo[[2]]
                print(paste("SLOCii", ltype, nloc, selectedLoci))
                print(paste("CUR SLICE", currentPosSlice))
                currentPosSliceEnd = currentPosSlice + selectedLoci[[currentNumLt]] - 1
                lt_vbyloc = valByloc_generic(addbylocFile, seq(nloc), paste(traitName, ltype, sep=" "), seq(currentPosSlice, currentPosSliceEnd))
                currentPosSlice = currentPosSliceEnd + 1
            }
        } else {
            print(paste("ltype",ltype,"absent."))
            ltypeQst = c(ltypeQst, "")
        }
        currentNumLt = currentNumLt + 1
    }
    names(ltypeQst) <- nameltype
    return(ltypeQst)
}

traits_Qst = sapply(argTraitsAdd, function(x) calc_qst_lt(argSpecies, x, argStep))

if (!is.na(traits_Qst))
{
    print(traits_Qst)

    for(t in argTraitsAdd)
    # in fact only one trait given at a time, no really need to go through a loop here
    {
        print(t)
        for (p in seq(3))
        {
            if (traits_Qst[p] != "")
		{
		    qstlt = paste("Qst",nameltype[p], sep="_")
		    print(paste("[F] :",substring(argTraits, 1, nchar(t)-4),"->",qstlt,":", traits_Qst[p]))
		}
	    }
    }
}
