[adegenet-commits] r624 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 11 10:10:43 CEST 2010
Author: jombart
Date: 2010-05-11 10:10:43 +0200 (Tue, 11 May 2010)
New Revision: 624
Modified:
pkg/R/haploGen.R
Log:
cleaning the code, renaming arguments consistently...
Modified: pkg/R/haploGen.R
===================================================================
--- pkg/R/haploGen.R 2010-05-11 07:47:51 UTC (rev 623)
+++ pkg/R/haploGen.R 2010-05-11 08:10:43 UTC (rev 624)
@@ -9,9 +9,9 @@
## mean.repro, sd.repro: average number of transmissions and its standard deviation (normal dist)
##
haploGen <- function(seq.length=1000, mu=0.0001,
- Tmax=50, mean.gen.time=5, sd.gen.time=1,
+ t.max=50, mean.gen.time=5, sd.gen.time=1,
mean.repro=2, sd.repro=1, max.nb.haplo=1e3,
- geo.sim=TRUE, grid.size=5, lambda.xy=0.5, matConnect=NULL,
+ geo.sim=TRUE, grid.size=5, lambda.xy=0.5, mat.connect=NULL,
ini.n=1, ini.xy=NULL){
## CHECKS ##
@@ -83,7 +83,7 @@
}
## where does a transmission occur (destination)?
- if(is.null(matConnect)){ # use universal lambda param
+ if(is.null(mat.connect)){ # use universal lambda param
xy.dupli <- function(cur.xy, nbLoc){
mvt <- rpois(2*nbLoc, lambda.xy) * sample(c(-1,1), size=2*nbLoc, replace=TRUE)
res <- t(matrix(mvt, nrow=2) + as.vector(cur.xy))
@@ -92,14 +92,14 @@
return(res)
}
} else { # use location-dependent proba of dispersal between locations
- if(any(matConnect < 0)) stop("Negative values in matConnect (probabilities expected!)")
- matConnect <- prop.table(matConnect,1)
+ if(any(mat.connect < 0)) stop("Negative values in mat.connect (probabilities expected!)")
+ mat.connect <- prop.table(mat.connect,1)
xy.dupli <- function(cur.xy, nbLoc){
- ## lambda.xy <- matConnect[cur.xy[1] , cur.xy[2]]
+ ## lambda.xy <- mat.connect[cur.xy[1] , cur.xy[2]]
## mvt <- rpois(2*nbLoc, lambda.xy) * sample(c(-1,1), size=2*nbLoc, replace=TRUE)
## res <- t(matrix(mvt, nrow=2) + as.vector(cur.xy))
idxAncesLoc <- myGrid[cur.xy[1], cur.xy[2]]
- newLoc <- sample(1:grid.size^2, size=nbLoc, prob=matConnect[idxAncesLoc,], replace=TRUE) # get new locations
+ newLoc <- sample(1:grid.size^2, size=nbLoc, prob=mat.connect[idxAncesLoc,], replace=TRUE) # get new locations
res <- cbind(row(myGrid)[newLoc] , col(myGrid)[newLoc]) # get coords of new locations
return(res)
}
@@ -153,7 +153,7 @@
nbDes <- nb.desc()
if(nbDes==0) return(NULL) # stop if no descendant
newDates <- sapply(1:nbDes, function(i) date.dupli(date)) # find dates for descendants
- newDates <- newDates[newDates <= Tmax] # don't store future sequences
+ newDates <- newDates[newDates <= t.max] # don't store future sequences
nbDes <- length(newDates)
if(nbDes==0) return(NULL) # stop if no suitable date
newSeq <- lapply(1:nbDes, function(i) seq.dupli(seq)) # generate new sequences
@@ -174,7 +174,7 @@
nbDes <- nb.desc()
if(nbDes==0) return(NULL) # stop if no descendant
newDates <- sapply(1:nbDes, function(i) date.dupli(date)) # find dates for descendants
- newDates <- newDates[newDates <= Tmax] # don't store future sequences
+ newDates <- newDates[newDates <= t.max] # don't store future sequences
nbDes <- length(newDates)
if(nbDes==0) return(NULL) # stop if no suitable date
newSeq <- lapply(1:nbDes, function(i) seq.dupli(seq)) # generate new sequences
@@ -227,9 +227,9 @@
## PERFORM SIMULATIONS - SPATIAL CASE ##
if(geo.sim){
## some checks
- if(!is.null(matConnect)) {
- if(nrow(matConnect) != ncol(matConnect)) stop("matConnect is not a square matrix")
- if(nrow(matConnect) != grid.size^2) stop("dimension of matConnect does not match grid size")
+ if(!is.null(mat.connect)) {
+ if(nrow(mat.connect) != ncol(mat.connect)) stop("mat.connect is not a square matrix")
+ if(nrow(mat.connect) != grid.size^2) stop("dimension of mat.connect does not match grid size")
}
## initialization
@@ -260,7 +260,7 @@
resize.result.xy()
## VERBOSE OUTPUT FOR DEBUGGING ##
## cat("\nNb strains:",length(res$ances),"/",max.nb.haplo)
- ## cat("\nLatest date:", max(res$dates),"/",Tmax)
+ ## cat("\nLatest date:", max(res$dates),"/",t.max)
## cat("\nRemaining strains to duplicate", sum(toExpand))
## cat("\n",append=TRUE,file="haploGenTime.out")
## iter.time <- as.numeric(difftime(Sys.time(),time.previous,unit="sec"))
@@ -437,7 +437,7 @@
################
## plotHaploGen
################
-plotHaploGen <- function(x, annot=FALSE, dateRange=NULL, col=NULL, bg="grey", add=FALSE, ...){
+plotHaploGen <- function(x, annot=FALSE, date.range=NULL, col=NULL, bg="grey", add=FALSE, ...){
## SOME CHECKS ##
if(class(x)!="haploGen") stop("x is not a haploGen object")
@@ -463,7 +463,7 @@
## CALL TO PLOTSEQTRACK ##
- plotSeqTrack(x=res, xy=xy, annot=annot, dateRange=dateRange,
+ plotSeqTrack(x=res, xy=xy, annot=annot, date.range=date.range,
col=col, bg=bg, add=add, ...)
return(invisible(res))
@@ -478,7 +478,8 @@
###################
## sample.haploGen
###################
-sample.haploGen <- function(x, n, rDate=.rTimeSeq, arg.rDate=NULL){
+sample.haploGen <- function(x, n){
+##sample.haploGen <- function(x, n, rDate=.rTimeSeq, arg.rDate=NULL){
## EXTRACT THE SAMPLE ##
res <- x[sample(1:nrow(x$seq), n, replace=FALSE)]
@@ -497,19 +498,17 @@
L <- 1000
}
- truedates <- res$dates
- daterange <- diff(range(res$dates,na.rm=TRUE))
+ ## truedates <- res$dates
+ ## daterange <- diff(range(res$dates,na.rm=TRUE))
- if(identical(rDate,.rTimeSeq)){
- sampdates <- .rTimeSeq(n=length(truedates), mu=mu0, L=L, maxNbDays=daterange/2)
- } else{
- arg.rDate$n <- n
- sampdates <- do.call(rDate, arg.rDate)
- }
- sampdates <- truedates + abs(sampdates)
+ ## if(identical(rDate,.rTimeSeq)){
+ ## sampdates <- .rTimeSeq(n=length(truedates), mu=mu0, L=L, maxNbDays=daterange/2)
+ ## } else{
+ ## arg.rDate$n <- n
+ ## sampdates <- do.call(rDate, arg.rDate)
+ ## }
+ ## sampdates <- truedates + abs(sampdates)
- res$dates <- sampdates
-
return(res)
} # end sample.haploGen
@@ -551,6 +550,17 @@
+
+
+
+
+
+
+
+
+
+
+
## #####################
## ## seqTrackG.haploGen
## #####################
More information about the adegenet-commits
mailing list