#!/usr/bin/env Rscript
#"Tirage des frequences initiales d'allèles à l'equilibre du modèle de migration en iles pour n loci nucléaires standard 
#[Goudet et Buchi, 2005]
#s nb de populations
#k nb total d'allèles
#kobs nb d'allèles initialement présents
#Nloc nb de loci




#Typical launch command example
#./dirichlet_standard.r filename="S1_genome.csv", Npop=55, Nm=10.2, Seed=1349 > test

#Nm related somehow to mp




###### Function for a better parsing of arguments, Can be safely ignored #######
###### From: http://cwcode.wordpress.com/2013/04/16/the-joys-of-rscript/ #######
getArgs <- function(verbose=FALSE, defaults=NULL) {
    myargs <- gsub("^--", "", commandArgs(TRUE))
    setopts <- !grepl("=", myargs)
    if(any(setopts))
        myargs[setopts] <- paste(myargs[setopts], "=notset", sep="")
    myargs.list <- strsplit(myargs, "=")
    myargs <- lapply(myargs.list, "[[", 2 )
    names(myargs) <- lapply(myargs.list, "[[" , 1)

    ## logicals
    if(any(setopts))
        myargs[setopts] <- TRUE

    ## defaults
    if(!is.null(defaults)) 
    {
        defs.needed <- setdiff(names(defaults), names(myargs))
        if(length(defs.needed)) 
        {
            myargs[ defs.needed ] <- defaults[ defs.needed ]
        }
    }

    ## verbage
    if(verbose) {
    cat("read",length(myargs),"named args:\n")
    print(myargs)
    }
    myargs
}
################################################################################





#Get args using the function above
myargs = getArgs()

# Reads the matrix of the specie
filename = myargs$filename
mat = read.table(filename, skip = 1 , sep = " ")

# Collects the number of locus and the number of alleles for each loci
Nlocus = dim(mat)[1]
Nallels = sapply(1:Nlocus, function(x) length(mat[x,!is.na(mat[x,])]))


# Contains dirichlet function
suppressMessages(library(MCMCpack))

# Some other param
set.seed(as.numeric(myargs$Seed)*100) #Random seed
s = as.numeric(myargs$Npop) #55
Nm = as.numeric(myargs$Nm )#10.2 #What is that
varname = myargs$varname

cat(varname, " = \n{", sep = "")
for(i in 1:Nlocus){  #For each locus
    #Lua compatible syntax
    cat("\n    {\n")
    #Vector of positions of selected locus
    allObs = which(mat[i,] == 1)
    #Number of selected locus
    kobs = length(allObs)
   
    #Creating enpty matrix
    pmat2 <- matrix(rep(0, s*Nallels[i]), nrow=s)
   
    if(kobs != 0)
    {
        #Generate distibution (matrix kobs*npop)
        P <- rdirichlet(1, rep(1,kobs))
        pmat <- rdirichlet(s, P*4*Nm)
        #Insert frequencies in the right place
        pmat2[,allObs] = pmat
    }

    #Print results
    for(j in 1:s){  #For each pop
        cat("        { ")
        cat(round(pmat2[j,]*100,2), sep = ", ")
        cat(" },\n")
    }
    cat("    },\n")
}
cat("}")
cat("\nreturn ", varname, sep = "")
