[adegenet-commits] r597 - in pkg: . R man misc/bug-report.1.2-3.02
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 15 19:38:03 CEST 2010
Author: jombart
Date: 2010-04-15 19:38:03 +0200 (Thu, 15 Apr 2010)
New Revision: 597
Added:
pkg/misc/bug-report.1.2-3.02/About-propShared.html
Removed:
pkg/misc/bug-report.1.2-3.02/About propShared.html
Modified:
pkg/DESCRIPTION
pkg/R/dapc.R
pkg/R/find.clust.R
pkg/R/haploGen.R
pkg/R/haploPop.R
pkg/R/seqTrack.R
pkg/man/dapc.Rd
pkg/man/find.clusters.Rd
pkg/man/seqTrack.Rd
Log:
On the way to making a new release. Commented lots of stuff under development and not used for now, with undecided status.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-04-15 17:05:01 UTC (rev 596)
+++ pkg/DESCRIPTION 2010-04-15 17:38:03 UTC (rev 597)
@@ -1,10 +1,10 @@
Package: adegenet
-Version: 1.3-0
-Date: 2009/04/01
+Version: 1.2-4
+Date: 2009/04/15
Title: adegenet: a R package for the multivariate analysis of genetic markers.
Author: Thibaut Jombart <t.jombart at imperial.ac.uk>, with contributions from Peter Solymos and Francois Balloux <f.balloux at imperial.ac.uk>
Maintainer: Thibaut Jombart <t.jombart at imperial.ac.uk>
-Suggests: ade4, genetics, hierfstat, spdep, tripack, ape
+Suggests: ade4, genetics, hierfstat, spdep, tripack, ape, pegas, graph, RBGL
Depends: methods
Description: Classes and functions for genetic data analysis within the multivariate framework.
License: GPL (>=2)
Modified: pkg/R/dapc.R
===================================================================
--- pkg/R/dapc.R 2010-04-15 17:05:01 UTC (rev 596)
+++ pkg/R/dapc.R 2010-04-15 17:38:03 UTC (rev 597)
@@ -8,7 +8,7 @@
#################
dapc.data.frame <- function(x, grp, n.pca=NULL, n.da=NULL,
center=TRUE, scale=FALSE, var.contrib=FALSE,
- pca.select=c("nbEig","percVar"), perc.pca=NULL){
+ pca.select=c("nbEig","percVar"), perc.pca=NULL, ...){
## FIRST CHECKS
if(!require(ade4, quiet=TRUE)) stop("ade4 library is required.")
@@ -124,7 +124,7 @@
#############
dapc.genind <- function(x, pop=NULL, n.pca=NULL, n.da=NULL,
scale=FALSE, scale.method=c("sigma", "binom"), truenames=TRUE, all.contrib=FALSE,
- pca.select=c("nbEig","percVar"), perc.pca=NULL){
+ pca.select=c("nbEig","percVar"), perc.pca=NULL, ...){
## FIRST CHECKS
if(!require(ade4, quiet=TRUE)) stop("ade4 library is required.")
Modified: pkg/R/find.clust.R
===================================================================
--- pkg/R/find.clust.R 2010-04-15 17:05:01 UTC (rev 596)
+++ pkg/R/find.clust.R 2010-04-15 17:38:03 UTC (rev 597)
@@ -7,7 +7,7 @@
## find.clusters.data.frame
######################
find.clusters.data.frame <- function(x, clust=NULL, n.pca=NULL, n.clust=NULL, stat=c("BIC", "AIC", "WSS"), choose.n.clust=TRUE, criterion=c("min","diff", "conserv"),
- max.n.clust=round(nrow(x)/10), n.iter=1e3, n.start=10, center=TRUE, scale=TRUE){
+ max.n.clust=round(nrow(x)/10), n.iter=1e3, n.start=10, center=TRUE, scale=TRUE, ...){
## CHECKS ##
if(!require(ade4, quiet=TRUE)) stop("ade4 library is required.")
Modified: pkg/R/haploGen.R
===================================================================
--- pkg/R/haploGen.R 2010-04-15 17:05:01 UTC (rev 596)
+++ pkg/R/haploGen.R 2010-04-15 17:38:03 UTC (rev 597)
@@ -421,24 +421,24 @@
-#####################
-## seqTrackG.haploGen
-#####################
-seqTrackG.haploGen <- function(x, optim=c("min","max"), ...){
- myX <- dist.dna(x$seq, model="raw")
- x.names <- labels(x)
- x.dates <- as.POSIXct(x)
- seq.length <- ncol(x$seq)
- myX <- myX * seq.length
- prevCall <- as.list(x$call)
- if(is.null(prevCall$mu)){
- mu0 <- 0.0001
- } else {
- mu0 <- eval(prevCall$mu)
- }
- res <- seqTrackG(myX, x.names=x.names, x.dates=x.dates, best=optim,...)
- return(res)
-}
+## #####################
+## ## seqTrackG.haploGen
+## #####################
+## seqTrackG.haploGen <- function(x, optim=c("min","max"), ...){
+## myX <- dist.dna(x$seq, model="raw")
+## x.names <- labels(x)
+## x.dates <- as.POSIXct(x)
+## seq.length <- ncol(x$seq)
+## myX <- myX * seq.length
+## prevCall <- as.list(x$call)
+## if(is.null(prevCall$mu)){
+## mu0 <- 0.0001
+## } else {
+## mu0 <- eval(prevCall$mu)
+## }
+## res <- seqTrackG(myX, x.names=x.names, x.dates=x.dates, best=optim,...)
+## return(res)
+## }
@@ -448,28 +448,28 @@
##############################
## optimize.seqTrack.haploGen
##############################
-optimize.seqTrack.haploGen <- function(x, thres=0.2, optim=c("min","max"),
- typed.chr=NULL, mu0=NULL, chr.length=NULL,
- prox.mat=NULL, nstep=10, step.size=1e3,
- rDate=.rTimeSeq, arg.rDate=NULL, rMissDate=.rUnifTimeSeq, ...){
+## optimize.seqTrack.haploGen <- function(x, thres=0.2, optim=c("min","max"),
+## typed.chr=NULL, mu0=NULL, chr.length=NULL,
+## prox.mat=NULL, nstep=10, step.size=1e3,
+## rDate=.rTimeSeq, arg.rDate=NULL, rMissDate=.rUnifTimeSeq, ...){
- x.names <- labels(x)
- x.dates <- as.POSIXct(x)
- seq.length <- ncol(x$seq)
- myX <- dist.dna(x$seq, model="raw") * seq.length
- prevCall <- as.list(x$call)
- if(is.null(prevCall$mu)){
- mu0 <- 0.0001
- } else {
- mu0 <- eval(prevCall$mu)
- }
+## x.names <- labels(x)
+## x.dates <- as.POSIXct(x)
+## seq.length <- ncol(x$seq)
+## myX <- dist.dna(x$seq, model="raw") * seq.length
+## prevCall <- as.list(x$call)
+## if(is.null(prevCall$mu)){
+## mu0 <- 0.0001
+## } else {
+## mu0 <- eval(prevCall$mu)
+## }
- res <- optimize.seqTrack.default(x=myX, x.names=x.names, x.dates=x.dates,
- typed.chr=typed.chr, mu0=mu0, chr.length=chr.length,
- thres=thres, optim=optim, prox.mat=prox.mat,
- nstep=nstep, step.size=step.size,
- rDate=rDate, arg.rDate=arg.rDate, rMissDate=rMissDate, ...)
-} # end optimize.seqTrack.haploGen
+## res <- optimize.seqTrack.default(x=myX, x.names=x.names, x.dates=x.dates,
+## typed.chr=typed.chr, mu0=mu0, chr.length=chr.length,
+## thres=thres, optim=optim, prox.mat=prox.mat,
+## nstep=nstep, step.size=step.size,
+## rDate=rDate, arg.rDate=arg.rDate, rMissDate=rMissDate, ...)
+## } # end optimize.seqTrack.haploGen
@@ -572,7 +572,7 @@
daterange <- diff(range(res$dates,na.rm=TRUE))
if(identical(rDate,.rTimeSeq)){
- sampdates <- .rTimeSeq(mu0=mu0, L=L, n=length(truedates), maxNbDays=daterange/2)
+ sampdates <- .rTimeSeq(n=length(truedates), mu=mu0, L=L, maxNbDays=daterange/2)
} else{
arg.rDate$n <- n
sampdates <- do.call(rDate, arg.rDate)
Modified: pkg/R/haploPop.R
===================================================================
--- pkg/R/haploPop.R 2010-04-15 17:05:01 UTC (rev 596)
+++ pkg/R/haploPop.R 2010-04-15 17:38:03 UTC (rev 597)
@@ -309,7 +309,7 @@
##################
## print.haploPop
##################
-print.haploPop <- function(object, ...){
+print.haploPop <- function(x, ...){
x <- object
myCall <- x$call
Modified: pkg/R/seqTrack.R
===================================================================
--- pkg/R/seqTrack.R 2010-04-15 17:05:01 UTC (rev 596)
+++ pkg/R/seqTrack.R 2010-04-15 17:38:03 UTC (rev 597)
@@ -5,14 +5,14 @@
UseMethod("seqTrack")
}
-seqTrackG <- function(...){
- UseMethod("seqTrackG")
-}
+## seqTrackG <- function(...){
+## UseMethod("seqTrackG")
+## }
-optimize.seqTrack <- function(...){
- UseMethod("optimize.seqTrack")
-}
+## optimize.seqTrack <- function(...){
+## UseMethod("optimize.seqTrack")
+## }
get.likelihood <- function(...){
@@ -74,7 +74,7 @@
rangeSize <- as.integer(difftime(dateMax,dateMin, units="days"))
nbDays <- round(runif(n, min=0, max=rangeSize))
res <- dateMin + nbDays*3600*24
- res <- as.POSIXct(round(res, units="days"))
+ res <- as.POSIXct(round.POSIXt(res, units="days"))
return(res)
}
@@ -326,132 +326,132 @@
-##############
-## seqTrackG - graph version of SeqTrack
-##############
-##
-## - x is a matrix giving weights x[i,j] such that:
-## 'i is ancestor j'; the algo looks for maximal weight branching
-##
-## - prox.mat is a directed proximity measure, so that prox.mat[i,j] is
-## the 'proximity when going from i to j'
-##
-seqTrackG.default <- function(x, x.names, x.dates, best=c("min","max"), force.temporal.order=TRUE,
- res.type=c("seqTrack", "graphNEL"), ...){
+## ##############
+## ## seqTrackG - graph version of SeqTrack
+## ##############
+## ##
+## ## - x is a matrix giving weights x[i,j] such that:
+## ## 'i is ancestor j'; the algo looks for maximal weight branching
+## ##
+## ## - prox.mat is a directed proximity measure, so that prox.mat[i,j] is
+## ## the 'proximity when going from i to j'
+## ##
+## seqTrackG.default <- function(x, x.names, x.dates, best=c("min","max"), force.temporal.order=TRUE,
+## res.type=c("seqTrack", "graphNEL"), ...){
- ## CHECKS ##
- if(!require("graph")) stop("the graph package is not installed")
- if(!require("RBGL")) stop("the RBGL package is not installed")
- if(!exists("edmondsOptimumBranching")) {
- stop("edmondsOptimumBranching does not exist; \nmake sure to use the latest Bioconductor (not CRAN) version of RBGL")
- cat("\nWould you like to try and install latest version of RBGL (needs internet connection)\n y/n: ")
- ans <- tolower(as.character(readLines(n = 1)))
- if(ans=="y"){
- source("http://bioconductor.org/biocLite.R")
- biocLite("RBGL")
- }
- }
+## ## CHECKS ##
+## if(!require("graph")) stop("the graph package is not installed")
+## if(!require("RBGL")) stop("the RBGL package is not installed")
+## if(!exists("edmondsOptimumBranching")) {
+## stop("edmondsOptimumBranching does not exist; \nmake sure to use the latest Bioconductor (not CRAN) version of RBGL")
+## cat("\nWould you like to try and install latest version of RBGL (needs internet connection)\n y/n: ")
+## ans <- tolower(as.character(readLines(n = 1)))
+## if(ans=="y"){
+## source("http://bioconductor.org/biocLite.R")
+## biocLite("RBGL")
+## }
+## }
- best <- match.arg(best)
- res.type <- match.arg(res.type)
+## best <- match.arg(best)
+## res.type <- match.arg(res.type)
- if(length(x.names) != length(x.dates)){
- stop("inconsistent length for x.dates")
- }
+## if(length(x.names) != length(x.dates)){
+## stop("inconsistent length for x.dates")
+## }
- if(is.character(x.dates)){
- msg <- paste("x.dates is a character vector; " ,
- "please convert it as dates using 'as.POSIXct'" ,
- "\n(making sure dates are given as 'YYYY/MM/DD' or 'YYYY-MM-DD').", sep="")
- stop(msg)
- }
+## if(is.character(x.dates)){
+## msg <- paste("x.dates is a character vector; " ,
+## "please convert it as dates using 'as.POSIXct'" ,
+## "\n(making sure dates are given as 'YYYY/MM/DD' or 'YYYY-MM-DD').", sep="")
+## stop(msg)
+## }
- x <- as.matrix(x)
+## x <- as.matrix(x)
- x.dates <- as.POSIXct(round.POSIXt(x.dates,units="days")) # round dates to the day
+## x.dates <- as.POSIXct(round.POSIXt(x.dates,units="days")) # round dates to the day
- if(length(x.names) != nrow(x)){
- stop("inconsistent dimension for x")
- }
+## if(length(x.names) != nrow(x)){
+## stop("inconsistent dimension for x")
+## }
- ## HANDLE BEST==MIN ##
- ## reverse x by a translation
- x.old <- x
- areZero <- x<1e-14
- x <- (max(x)+1) - x
- x[areZero] <- 0
- diag(x) <- 0
+## ## HANDLE BEST==MIN ##
+## ## reverse x by a translation
+## x.old <- x
+## areZero <- x<1e-14
+## x <- (max(x)+1) - x
+## x[areZero] <- 0
+## diag(x) <- 0
- ## BUILD THE GRAPH ##
- ## make prox matrix with temporal consistency
- if(force.temporal.order){
- x[outer(myDates, myDates, ">=")] <- 0
- }
+## ## BUILD THE GRAPH ##
+## ## make prox matrix with temporal consistency
+## if(force.temporal.order){
+## x[outer(myDates, myDates, ">=")] <- 0
+## }
- ## tweak to get around a bug in as(...,"graphNEL") - looses edge weights
- ## to replace with commented line once fixed in CRAN release
- ## myGraph <- as(x, "graphNEL")
- myGraph <- as(new("graphAM", x, edgemode="directed", values=list(weight=1)), "graphNEL")
+## ## tweak to get around a bug in as(...,"graphNEL") - looses edge weights
+## ## to replace with commented line once fixed in CRAN release
+## ## myGraph <- as(x, "graphNEL")
+## myGraph <- as(new("graphAM", x, edgemode="directed", values=list(weight=1)), "graphNEL")
- ## CALL EDMONDSOPTIMUMBRANCHING ##
- temp <- edmondsOptimumBranching(myGraph)
+## ## CALL EDMONDSOPTIMUMBRANCHING ##
+## temp <- edmondsOptimumBranching(myGraph)
- ## SHAPE OUTPUT ##
- if(res.type=="seqTrack"){
- ## reorder output from edmondsOptimumBranching
- N <- length(x.names)
- myLev <- x.names
- ances <- as.integer(factor(temp$edgeList["from",], levels=myLev))
- desc <- as.integer(factor(temp$edgeList["to",], levels=myLev))
- newOrd <- order(desc)
- desc <- 1:N
- ances <- ances[newOrd]
- weights <- as.numeric(temp$weights)[newOrd]
+## ## SHAPE OUTPUT ##
+## if(res.type=="seqTrack"){
+## ## reorder output from edmondsOptimumBranching
+## N <- length(x.names)
+## myLev <- x.names
+## ances <- as.integer(factor(temp$edgeList["from",], levels=myLev))
+## desc <- as.integer(factor(temp$edgeList["to",], levels=myLev))
+## newOrd <- order(desc)
+## desc <- 1:N
+## ances <- ances[newOrd]
+## weights <- as.numeric(temp$weights)[newOrd]
- ## create the data.frame
- res <- data.frame(id=1:N)
- hasNoAnces <- difftime(myDates, min(myDates), units="secs") < 1 # 1 sec resolution for dates
- res$weight <- res$ances <- 1:N
- res$date <- x.dates
+## ## create the data.frame
+## res <- data.frame(id=1:N)
+## hasNoAnces <- difftime(myDates, min(myDates), units="secs") < 1 # 1 sec resolution for dates
+## res$weight <- res$ances <- 1:N
+## res$date <- x.dates
- ## fill in the d.f. with correct values
- res$ances[!hasNoAnces] <- ances
- res$ances[hasNoAnces] <- NA
+## ## fill in the d.f. with correct values
+## res$ances[!hasNoAnces] <- ances
+## res$ances[hasNoAnces] <- NA
- res$weight[!hasNoAnces] <- weights
- res$weight[hasNoAnces] <- NA
+## res$weight[!hasNoAnces] <- weights
+## res$weight[hasNoAnces] <- NA
- res$ances.date <- x.dates[res$ances]
+## res$ances.date <- x.dates[res$ances]
- res[is.na(res)] <- NA # have clean NAs
- row.names(res) <- x.names
- class(res) <- c("seqTrack","data.frame")
+## res[is.na(res)] <- NA # have clean NAs
+## row.names(res) <- x.names
+## class(res) <- c("seqTrack","data.frame")
- ## handle best==min
- if(best=="min"){
- res$weight <- max(x.old)+1 - res$weight
- }
+## ## handle best==min
+## if(best=="min"){
+## res$weight <- max(x.old)+1 - res$weight
+## }
- }
+## }
- if(res.type=="graphNEL"){
- ## handle optim==min
- if(best=="min"){
- temp$weights <- max(x.old)+1 - temp$weights
- }
- res <- ftM2graphNEL(t(temp$edgeList), W=temp$weights, edgemode="directed")
- }
+## if(res.type=="graphNEL"){
+## ## handle optim==min
+## if(best=="min"){
+## temp$weights <- max(x.old)+1 - temp$weights
+## }
+## res <- ftM2graphNEL(t(temp$edgeList), W=temp$weights, edgemode="directed")
+## }
- ## RETURN RESULT ##
- return(res)
+## ## RETURN RESULT ##
+## return(res)
-} # end seqTrackG
+## } # end seqTrackG
@@ -682,404 +682,404 @@
## 2) VECTORIZE mu and seq.length, recycle if needed with a warning
## 3) uncomment, adapt, and test code for missing data
##
-optimize.seqTrack.default <- function(x, x.names, x.dates, typed.chr=NULL, mu=NULL, seq.length=NULL,
- thres=0.2, best=c("min","max"), prox.mat=NULL, nstep=10, step.size=1e3,
- rDate=.rTimeSeq, arg.rDate=NULL, rMissDate=.rUnifDate, ...){
+## optimize.seqTrack.default <- function(x, x.names, x.dates, typed.chr=NULL, mu=NULL, seq.length=NULL,
+## thres=0.2, best=c("min","max"), prox.mat=NULL, nstep=10, step.size=1e3,
+## rDate=.rTimeSeq, arg.rDate=NULL, rMissDate=.rUnifDate, ...){
- ## CHECKS ##
- best <- match.arg(best)
- if(best=="min"){
- which.best <- which.min
- } else {
- which.best <- which.max
- }
+## ## CHECKS ##
+## best <- match.arg(best)
+## if(best=="min"){
+## which.best <- which.min
+## } else {
+## which.best <- which.max
+## }
- if(length(x.names) != length(x.dates)){
- stop("inconsistent length for x.dates")
- }
+## if(length(x.names) != length(x.dates)){
+## stop("inconsistent length for x.dates")
+## }
- if(is.character(x.dates)){
- msg <- paste("x.dates is a character vector; " ,
- "please convert it as dates using 'as.POSIXct'" ,
- "\n(making sure dates are given as 'YYYY/MM/DD' or 'YYYY-MM-DD').", sep="")
- stop(msg)
- }
+## if(is.character(x.dates)){
+## msg <- paste("x.dates is a character vector; " ,
+## "please convert it as dates using 'as.POSIXct'" ,
+## "\n(making sure dates are given as 'YYYY/MM/DD' or 'YYYY-MM-DD').", sep="")
+## stop(msg)
+## }
- isMissDate <- is.na(x.dates)
+## isMissDate <- is.na(x.dates)
- if(!identical(rDate, .rTimeSeq)){
- if(is.null(arg.rDate)){
- warning("Specific time distribution specified without arguments.")
- arg.rDate <- list(n=step.size)
- } else {
- if(!is.list(arg.rDate)) stop("If provided, arg.rDate must be a list.")
- if(!is.null(arg.rDate$n)) {
- warning("arg.rDate$n is provided, but will be replaced by step.size.")
- }
- arg.rDate$n <- step.size
- }
- }
+## if(!identical(rDate, .rTimeSeq)){
+## if(is.null(arg.rDate)){
+## warning("Specific time distribution specified without arguments.")
+## arg.rDate <- list(n=step.size)
+## } else {
+## if(!is.list(arg.rDate)) stop("If provided, arg.rDate must be a list.")
+## if(!is.null(arg.rDate$n)) {
+## warning("arg.rDate$n is provided, but will be replaced by step.size.")
+## }
+## arg.rDate$n <- step.size
+## }
+## }
- N <- length(x.names)
- id <- 1:N
- ## if(length(mu) < N) { # recycle mu
- ## mu <- rep(mu, length=N)
- ## }
- ## if(length(seq.length) < N) {# recycle seq.length
- ## seq.length <- rep(seq.length, length=N)
- ## }
+## N <- length(x.names)
+## id <- 1:N
+## ## if(length(mu) < N) { # recycle mu
+## ## mu <- rep(mu, length=N)
+## ## }
+## ## if(length(seq.length) < N) {# recycle seq.length
+## ## seq.length <- rep(seq.length, length=N)
+## ## }
- ## handle typed.chr, mu, seq.length
- if(identical(rDate, .rTimeSeq)){
- if(is.null(typed.chr)|is.null(mu)|is.null(seq.length)){
- stop("typed.chr, mu, and seq.length must be provided if rDate is .rTimeSeq")
- }
+## ## handle typed.chr, mu, seq.length
+## if(identical(rDate, .rTimeSeq)){
+## if(is.null(typed.chr)|is.null(mu)|is.null(seq.length)){
+## stop("typed.chr, mu, and seq.length must be provided if rDate is .rTimeSeq")
+## }
- if(!is.list(typed.chr)) {
- stop("typed.chr must be a list")
- }
- if(length(typed.chr)!=N) {
- stop("typed.chr has an inconsistent length")
- }
+## if(!is.list(typed.chr)) {
+## stop("typed.chr must be a list")
+## }
+## if(length(typed.chr)!=N) {
+## stop("typed.chr has an inconsistent length")
+## }
- if(is.null(names(mu))) stop("mu has no names")
- if(is.null(names(seq.length))) stop("seq.length has no names")
- if(any(mu > 1)) stop("mu has values > 1")
- if(any(mu < 0)) stop("mu has negative values")
+## if(is.null(names(mu))) stop("mu has no names")
+## if(is.null(names(seq.length))) stop("seq.length has no names")
+## if(any(mu > 1)) stop("mu has values > 1")
+## if(any(mu < 0)) stop("mu has negative values")
- if(!identical(names(mu) , names(seq.length))) stop("Names of mu and seq.length differ.")
- if(any(!unique(unlist(typed.chr)) %in% names(mu))) {
- stop("Some chromosomes indicated in typed.chr are not in mu.")
- }
+## if(!identical(names(mu) , names(seq.length))) stop("Names of mu and seq.length differ.")
+## if(any(!unique(unlist(typed.chr)) %in% names(mu))) {
+## stop("Some chromosomes indicated in typed.chr are not in mu.")
+## }
- list.mu <- lapply(typed.chr, function(e) mu[e])
- list.seq.length <- lapply(typed.chr, function(e) seq.length[e])
- }
+## list.mu <- lapply(typed.chr, function(e) mu[e])
+## list.seq.length <- lapply(typed.chr, function(e) seq.length[e])
+## }
- x.dates <- as.POSIXct(round.POSIXt(x.dates,units="days")) # round dates to the day
+## x.dates <- as.POSIXct(round.POSIXt(x.dates,units="days")) # round dates to the day
- x <- as.matrix(x)
+## x <- as.matrix(x)
- if(!is.null(prox.mat) && !identical(dim(prox.mat),dim(x))) {
- stop("prox.mat is provided but its dimensions are inconsistent with that of x")
- }
+## if(!is.null(prox.mat) && !identical(dim(prox.mat),dim(x))) {
+## stop("prox.mat is provided but its dimensions are inconsistent with that of x")
+## }
- ## rename dimensions using id
- colnames(x) <- rownames(x) <- id
+## ## rename dimensions using id
+## colnames(x) <- rownames(x) <- id
- if(length(x.names) != nrow(x)){
- stop("inconsistent dimension for x")
- }
+## if(length(x.names) != nrow(x)){
+## stop("inconsistent dimension for x")
+## }
- ## SET THRESHOLD IF NEEDED ## ## NO LONGER USED
- ## if(is.null(thres)){
- ## thres <- sum(seqTrack(x.names=x.names, x.dates=x.dates, W=W,
- ## best=best, prox.mat=prox.mat, ...)$weight, na.rm=TRUE)
- ## }
+## ## SET THRESHOLD IF NEEDED ## ## NO LONGER USED
+## ## if(is.null(thres)){
+## ## thres <- sum(seqTrack(x.names=x.names, x.dates=x.dates, W=W,
+## ## best=best, prox.mat=prox.mat, ...)$weight, na.rm=TRUE)
+## ## }
- ## AUXILIARY FUNCTIONS ##
+## ## AUXILIARY FUNCTIONS ##
- ## to compare results -> returns a list of length two: logical, and the value of the res
- val.res <- function(res){
- return(sum(res$weight, na.rm=TRUE))
- }
+## ## to compare results -> returns a list of length two: logical, and the value of the res
+## val.res <- function(res){
+## return(sum(res$weight, na.rm=TRUE))
+## }
- ## DO THE OPTIMISATION ##
- RANGE.DATES <- as.integer(round(diff(range(x.dates, na.rm=TRUE)))) # time window of the sample, in days
- NB.DATES.TO.SIM <- sum(!isMissDate)
+## ## DO THE OPTIMISATION ##
+## RANGE.DATES <- as.integer(round(diff(range(x.dates, na.rm=TRUE)))) # time window of the sample, in days
+## NB.DATES.TO.SIM <- sum(!isMissDate)
- ## for loop is not to slow for < 1e6 rep
- ## and allows not to handle huge objects
- ## (which would grow exponentially)
+## ## for loop is not to slow for < 1e6 rep
+## ## and allows not to handle huge objects
+## ## (which would grow exponentially)
- ## res.best <- res.ini # initialization
+## ## res.best <- res.ini # initialization
- ## DEFINE OUTPUTS ##
- ances <- integer(0)
- date <- character(0)
- ances.date <- character(0)
- valRes <- numeric(0)
+## ## DEFINE OUTPUTS ##
+## ances <- integer(0)
+## date <- character(0)
+## ances.date <- character(0)
+## valRes <- numeric(0)
- ## DEFAULT CASE: NO MISSING DATES
- if(!any(isMissDate)){
- ## dates initialisation, taken from initial prior
- ## If dates distrib is .rTimeSeq
- if(identical(rDate, .rTimeSeq)){
- newDates <- sapply(1:N, function(i)
- rDate(n=step.size, mu=list.mu[[i]], L=list.seq.length[[i]],
- maxNbDays=RANGE.DATES))
- } else { ## Else, any other distrib with free arguements
- newDates <- sapply(1:N, function(i) do.call(rDate, arg.rDate))
- }
+## ## DEFAULT CASE: NO MISSING DATES
+## if(!any(isMissDate)){
+## ## dates initialisation, taken from initial prior
+## ## If dates distrib is .rTimeSeq
+## if(identical(rDate, .rTimeSeq)){
+## newDates <- sapply(1:N, function(i)
+## rDate(n=step.size, mu=list.mu[[i]], L=list.seq.length[[i]],
+## maxNbDays=RANGE.DATES))
+## } else { ## Else, any other distrib with free arguements
+## newDates <- sapply(1:N, function(i) do.call(rDate, arg.rDate))
+## }
- newDates <- t(newDates)*24*3600 + x.dates
+## newDates <- t(newDates)*24*3600 + x.dates
- ## >> one step of 'step.size' simulations, all with same prior << ##
- for(i in 1:nstep){
- ## >> each step contains 'step.size' iterations << ##
- for(j in 1:step.size){
- myDates <- as.POSIXct(newDates[,j])
+## ## >> one step of 'step.size' simulations, all with same prior << ##
+## for(i in 1:nstep){
+## ## >> each step contains 'step.size' iterations << ##
+## for(j in 1:step.size){
+## myDates <- as.POSIXct(newDates[,j])
- res.new <- seqTrack(x, x.names=x.names, x.dates=myDates,
- best=best, prox.mat=prox.mat, ...)
+## res.new <- seqTrack(x, x.names=x.names, x.dates=myDates,
+## best=best, prox.mat=prox.mat, ...)
- ##ances <- cbind(ances, res.new$ances) # not needed now
- date <- cbind(date, as.character(res.new$date))
- ##ances.date <- cbind(ances.date, as.character(res.new$ances.date)) # not needed now
- valRes <- c(valRes, val.res(res.new))
- ##}
- } # end for j
+## ##ances <- cbind(ances, res.new$ances) # not needed now
+## date <- cbind(date, as.character(res.new$date))
+## ##ances.date <- cbind(ances.date, as.character(res.new$ances.date)) # not needed now
+## valRes <- c(valRes, val.res(res.new))
+## ##}
+## } # end for j
- ## retain a given % (thres) of the dates ##
- toKeep <- valRes <= quantile(valRes, thres) ## NOT WORKING FOR optim==max !!!
- valRes <- valRes[toKeep]
+## ## retain a given % (thres) of the dates ##
+## toKeep <- valRes <= quantile(valRes, thres) ## NOT WORKING FOR optim==max !!!
+## valRes <- valRes[toKeep]
- date <- date[,toKeep,drop=FALSE] # retained posterior
+## date <- date[,toKeep,drop=FALSE] # retained posterior
- ## DEBUGING ##
- ## cat("\ntoKeep:\n")
- ## print(toKeep)
- ## cat("\nhead date (posterior):\n")
- ## print(head(date))
- ## END DEBUGING ##
+## ## DEBUGING ##
+## ## cat("\ntoKeep:\n")
+## ## print(toKeep)
+## ## cat("\nhead date (posterior):\n")
+## ## print(head(date))
+## ## END DEBUGING ##
- newDates <- apply(date, 1, function(vec)
- sample(vec, size=step.size, replace=TRUE)) # new prior
- newDates <- t(newDates)
+## newDates <- apply(date, 1, function(vec)
+## sample(vec, size=step.size, replace=TRUE)) # new prior
+## newDates <- t(newDates)
- ## stop if all dates are fixed
- if(all(apply(newDates, 1, function(r) length(unique(r))==1))){
- cat("\nConvergence reached at step",i,"\n")
- break # stop the algorithm
- }
+## ## stop if all dates are fixed
+## if(all(apply(newDates, 1, function(r) length(unique(r))==1))){
+## cat("\nConvergence reached at step",i,"\n")
+## break # stop the algorithm
+## }
- ## re-initialize posterior distributions
- if(i<nstep){
- ## ances <- integer(0) # not needed now
- date <- character(0)
- ## ances.date <- character(0) # not needed now
- valRes <- numeric(0)
- } # end if
+## ## re-initialize posterior distributions
+## if(i<nstep){
+## ## ances <- integer(0) # not needed now
+## date <- character(0)
+## ## ances.date <- character(0) # not needed now
+## valRes <- numeric(0)
+## } # end if
- } # end for i
+## } # end for i
- ## ## dates: new prior taken from obtained posterior
- ## if(length(valRes)==0) { # if no simul are retained
- ## warning(paste("No simulation was retained at the given threshold at step",i))
- ## } else {
- ## if(optim=="min"){ # define weights for further samplings
- ## w <- max(valRes,na.rm=TRUE) - valRes
- ## w <- w/sum(w)
- ## } else {
- ## w <- valRes
- ## w <- w/sum(w)
- ## }
+## ## ## dates: new prior taken from obtained posterior
+## ## if(length(valRes)==0) { # if no simul are retained
+## ## warning(paste("No simulation was retained at the given threshold at step",i))
+## ## } else {
+## ## if(optim=="min"){ # define weights for further samplings
+## ## w <- max(valRes,na.rm=TRUE) - valRes
+## ## w <- w/sum(w)
+## ## } else {
+## ## w <- valRes
+## ## w <- w/sum(w)
+## ## }
- ## newDates <- apply(date, 1, function(vec) # used a weighted sampling
- ## sample(vec, size=step.size, replace=TRUE, prob=w))
+## ## newDates <- apply(date, 1, function(vec) # used a weighted sampling
+## ## sample(vec, size=step.size, replace=TRUE, prob=w))
- } # end if(!any(isMissDate))
+## } # end if(!any(isMissDate))
- ## ## OTHER CASE: HANDLE MISSING DATES
- ## if(any(isMissDate)){
+## ## ## OTHER CASE: HANDLE MISSING DATES
+## ## if(any(isMissDate)){
- ## ## Handle distribution and its parameters ##
- ## argList <- list(...)
+## ## ## Handle distribution and its parameters ##
+## ## argList <- list(...)
- ## if(is.null(argList$dateMin) & identical(rMissDate, .rUnifDate)){ # earliest date
- ## argList$dateMin <- min(x.dates,na.rm=TRUE)
- ## } else {
- ## argList$dateMin[is.na(argList$dateMin)] <- min(x.dates,na.rm=TRUE)
- ## }
- ## if(is.null(argList$dateMax) & identical(rMissDate, .rUnifDate)){ # latest date
- ## argList$dateMax <- max(x.dates,na.rm=TRUE)
- ## } else {
- ## argList$dateMax[is.na(argList$dateMax)] <- max(x.dates,na.rm=TRUE)
- ## }
+## ## if(is.null(argList$dateMin) & identical(rMissDate, .rUnifDate)){ # earliest date
+## ## argList$dateMin <- min(x.dates,na.rm=TRUE)
+## ## } else {
+## ## argList$dateMin[is.na(argList$dateMin)] <- min(x.dates,na.rm=TRUE)
+## ## }
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/adegenet -r 597
More information about the adegenet-commits
mailing list