stdupas / EnvironmentalDemogeneticsABC

Statistical models of coalescent on a graph
0 stars 3 forks source link

Mutations code pieces #8

Open Becheler opened 9 years ago

Becheler commented 9 years ago

function to build the resultant :

resultant <- function(nbrMutations, mutationModel, args){
args <- c(nbrMutations, args)
return(do.call(what = mutationModel, args = args))
}

and, I think the apply on the table would be something like :

resultantOnCoalTab <- function(coalTable, step){
coalTable[Resultant] <- apply(X=coaltable[,"MutationNbr"], FUN=resultant) * step
return(coalTable)
}
stdupas commented 9 years ago
stepWiseMutationModel <- function(mutations)
{
  # Simulates a change in the genetic value depending on number of mutations
  # assuming stepwise mutation model
  # Args: 
  #   mutations= number of mutations
  #
  # Value
  #   The resultant in microsatellite repetition number using binomial rules
  # 
  # Example:
  # stepWiseMutationModel(5)
  #
  2*rbinom(length(mutations),mutations,.5)-mutations
}
Becheler commented 9 years ago

Yop la syntaxe suivante ne change rien, mais j'ai tendance à éviter tout ce qui est implicite, et donc source d'erreurs ou de manque de lisibilité (expliciter return, expliciter les noms d'arguments pour toutes les fonctions etc...)

stepWiseMutationModel <- function(mutations)
{
  # Compute a resultant according to step wise mutation model
  #
  # Args: 
  #   mutations: number of mutations
  #
  # Returns:
  #   The resultant in microsatellite repetition number using binomial rules

  res <- 2*rbinom(n = length(mutations), size = mutations, prob = 0.5) - mutations
  return(res)

  # Example:
  # stepWiseMutationModel(mutations=5)
  # stepWiseMutationModel(mutations=c(2,3))
}

Je ne sais pas si c'est à dessein, mais dans la version antérieur du bout de code, il y avait un /2 qui a sauté ici :

res <- 2*rbinom(n = length(mutations), size = mutations, prob = 0.5) - mutations
  res <- 2*rbinom(n = length(mutations), size = mutations, prob = 0.5) - mutations/2

Mais dans tous les cas on dirait que ça fonctionne, dans le sens où tu récupères un vecteur de résultantes de même longueur que le vecteur d'entrée. Du coup il est possible de faire :

 coalTable[[,"Resultant"]] <- stepWiseMutationModel(coalTable[["Mutations"]]

ce qui revient à faire un :

mutationModel <- "stepWiseMutationModel"
args <- list()
 coalTable[[,"Resultant"]] <- do.call(what = mutationModel, 
                                                      args = c(coalTable[["Mutations"]], args))

spécifier une liste vide peut sembler inutile pour ce modèle précis, mais c'est pour ne pas avoir à fournir une fonction explicite de mutation (MutationModel correspond en fait au modèle spécifié par l'utilisateur qu'on ira chopper dans ParamList, avec ses arguments-prior stockés dans args, qu'on concatène avec la partie manquante de arguments (à savoir le nombre de mutations).

et je propose de cacher toute cette méchante machinerie dans une fonction toute propre qu'on peut appeller dans le main :

resultantFunction <- function(coalTable, mutationModel, args){
 coalTable[[,"Resultant"]] <- do.call(what = mutationModel, 
                                                     args = c(coalTable[["Mutations"]], args))
}

Et ça devrait fonctionner, à deux ou trois ajustement syntaxiques près...

stdupas commented 9 years ago

OK super je push ce truc alors 1) Pour le /2 il s'appliquait mutations, qu'on multipliait ensuite par 2 dans la parenthèse, donc 2/2=1. Je coirs que c'est bon. coaltable$resultant = 2_(rbinom(dim(coaltable)[1],coaltable[,"mutations"],.5)-coaltable[,"mutations"]/2) en fait c'est pareil que res <- 2_rbinom(n = length(mutations), size = mutations, prob = 0.5) - mutations 2) pour le return(res), je pensais qu'il fallait mieux pas avoir à faire une attribution de variable inutile. Ca gagne du temps. Mais c'est du détail. Pour l'instant faisons les choses claires. 3) Pour le resultant function: pourquoi le mettre dans le main? 4) Raccroche le téléphone

Becheler commented 9 years ago

4) Lol décidément je ne sais pas me servir de ce téléphone, je l'ai raccroché 46 fois... C'est bizarre ce truc. Sacré téléphone, je ne m'y fais pas. Premier essai pour t'appeler, il voulait m'envoyer chez la police ... 1) Oki, j'avais pas checké l'ensemble de l'expression, j'ai juste vu que ça avait sauté :) Donc c'était finalement à dessein :) 2) Pour la déclaration de variable, j'avais choisi de la créer pour pas que l'expression soit difficilement lisible si on fait un return( 2rbinom(n = length(mutations), size = mutations, prob = 0.5) - mutations)

3) Appeller la fonction dans le main pour qu'elle agisse sur l'objet coalescent. Pas déclarer la fonction dans le main (j'étais pas clair, désolé, j'avais corrigé mais après coup).