[adegenet-commits] r305 - in pkg: . R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri May 8 16:49:57 CEST 2009


Author: jombart
Date: 2009-05-08 16:49:57 +0200 (Fri, 08 May 2009)
New Revision: 305

Modified:
   pkg/R/seqTrack.R
   pkg/TODO
Log:
 
sasd


Modified: pkg/R/seqTrack.R
===================================================================
--- pkg/R/seqTrack.R	2009-04-30 21:05:38 UTC (rev 304)
+++ pkg/R/seqTrack.R	2009-05-08 14:49:57 UTC (rev 305)
@@ -1,9 +1,11 @@
-#############
-## seqTrack
-#############
-seqTrack <- function(seq.names, seq.dates, D, k=5, lag=3, ...){
+## #############
+## ## seqTrack
+## #############
+seqTrack <- function(seq.names, seq.dates, D, ...){
+    ## CHECKS ##
+    if(!require(graph)) stop("Package graph is not installed")
+    if(!require(RBGL)) stop("Package RBGL is not installed")
 
-    ## CHECKS ##
     if(length(seq.names) != length(seq.dates)){
         stop("inconsistent length for seq.dates")
     }
@@ -14,127 +16,194 @@
         stop("inconsistent dimension for D")
     }
 
+     N <- length(seq.names)
+     id <- 1:N
+    std.time <- as.numeric(seq.dates) # turn dates into numeric
+    std.time <- std.time/sum(std.time) # set max to 1
 
-    ## ASSIGNEMENTS ##
-    N <- length(seq.names)
-    STARTPOINTS <- 1:N # global variable, modified thoughout
-    id <- 1:N
-    INPATH <- NULL
-    CURRENTPATH <- list() # current path
-    CURRENTPATHDIST <- list() # current path distances
-    listPaths <- list() # pre-final output
-    listPathsDist <- list() # pre-final output
-    rownames(D) <- id
-    colnames(D) <- id
 
+    ## BUILD MATRIX OF WEIGHTS ##
+    dist.time <- as.matrix(dist(std.time))
+    W <- D * dist.time
 
-    ## AUXILIARY FUNCTIONS ##
-    ## choose a starting sequence, as recent as possible
-    chooseStartPoint <- function(){
-        return(STARTPOINTS[which.max(seq.dates[STARTPOINTS])])
+
+    ## CONSTRUCTION OF THE GRAPH ##
+    ## build connections for one point
+    prepGraphOneNode <- function(idx){
+        res <- list(edges, weights)
+        res$edges <- which(seq.dates < seq.dates[idx])
+        res$weights <- W[idx, res$edges]
+        return(res)
     }
 
+    ## build the edge list + weights
+    edL <- lapply(id, prepGraphOneNode)
+    names(edL) <- seq.names
 
-    ## find k closest ancestors: returns a list(id=[id of the ancestors], d=[distances to ancestors])
-    findAncestors <- function(currentPoint){
-        candidates <- id[seq.dates < seq.dates[currentPoint]]
-        if(length(candidates)==0) return(NULL) # this will indicate the end of the path
-        res <- D[currentPoint,candidates]
-        if(length(res) <= k) return(list(id=as.integer(names(res)), d=res)) # if less than k candidates
-        res <- sort(res)[1:k] # if more than k, take k closest candidates
-        return(list(id=as.integer(names(res)), d=res))
-    }
+    myGraph <- new("graphNEL", nodes=seq.names, edgeL=edL, edgemode="directed")
 
 
-    ## discard paths stemming from the worst ancestors a step ...
-    discardPath <- function(step){
-        tempDist <- sapply(CURRENTPATHDIST, sum)
-        pathFac <- sapply(CURRENTPATH, function(e) e[step])
-        pathLengths <- tapply(tempDist, pathFac, sum)
-        toKeep <- names(pathLengths)[which.min(temp)]
-        toKeep <- as.integer(toKeep)
-        CURRENTPATH <<- CURRENTPATH[pathFac==toKeep]
-        CURRENTPATHDIST <<- CURRENTPATHDIST[pathFac==toKeep]
-        return()
-    } # end discardPath
+    ## RUN DIJKSTRA ##
+    tempRes <- lapply(seq.names, function(node) sp.between(myGraph, node, seq.names[which.min(seq.dates)]) )
 
 
-    ## update id in INPATH
-    INPATH.up <- function(){
-        temp <- unique(unlist(listPaths))
-        INPATH <<- union(INPATH,temp)
-        return()
-    } # end INPATH.up
+    ## RETURN THE RESULT ##
 
 
-    ## retrieve an already known path
-    findExistingPath <- function(id){
-        if(id %in% INPATH){
-            temp <- sapply(listPaths, function(e) id %in% e) # find in which path
-            res <- listPaths[[which(temp)[1]]] # retrieve the path
-            res <- res[which(res==id):length(res)] # cut the path
-            return(res)
-        } else return(NULL)
-    } # end findExistingPath
+}
 
 
-    ## expand one path
-    expandOnePath <- function(onePath){ # onePath has a single path (a vector)
-        newPoints <- findAncestors(onePath[length(onePath)]) # get new points
-        if(is.null(newPoints)) return(onePath)
-        res <- lapply(1:length(newPoints), function(i) c(onePath, newPoints)) # duplicate path and add new pts
-        return(res)
-    } # end expandOnePath
 
 
-    ## check if CURRENTPATH should stop (TRUE if yes, FALSE otherwise)
-    checkEndCurrentPath <- function(){
-        isEnded <- function(onePath){
-            oldest <- onePath[length(onePath)]
-            temp <- setdiff(id,oldest)
-            if(all(seq.dates[temp] < seq.dates[oldest])) return(TRUE)
-            return(FALSE)
-        }
 
-        temp <- sapply(CURRENTPATH, checkEndCurrentPath)
-        if(all(temp)) return(TRUE) # if all paths are ended
-        return(FALSE)
-    } # end checkEndCurrentPath
 
 
 
-    ## FIND ONE PATH ##
-    findPath <- function(id){
-        temp <- findExistingPath(id) # search for an already existing path
-        if(!is.null(temp)) return(temp) # return existing path if needed
 
-        ## WHILE LOOP ##
-        keepSearching <- TRUE
-        while(keepSearching){
-            for(i in length(CURRENTPATH)){
-                temp <- CURRENPATH[[1]]
 
-            }
 
 
+#########################
+## OLD ADD-HOC VERSION
+#########################
 
-            temp <- findExistingPath(id) # search for an already existing path
-            if(!is.null(temp))
-            keepSearching <- !checkEndCurrentPath() # stop searching if all paths are ended
+## #############
+## ## seqTrack
+## #############
+## seqTrack <- function(seq.names, seq.dates, D, k=5, lag=3, ...){
 
+##     ## CHECKS ##
+##     if(length(seq.names) != length(seq.dates)){
+##         stop("inconsistent length for seq.dates")
+##     }
 
+##     D <- as.matrix(D)
 
+##     if(length(seq.names) != nrow(D)){
+##         stop("inconsistent dimension for D")
+##     }
 
 
-        }
+##     ## ASSIGNEMENTS ##
+##     N <- length(seq.names)
+##     STARTPOINTS <- 1:N # global variable, modified thoughout
+##     id <- 1:N
+##     INPATH <- NULL
+##     CURRENTPATH <- list() # current path
+##     CURRENTPATHDIST <- list() # current path distances
+##     listPaths <- list() # pre-final output
+##     listPathsDist <- list() # pre-final output
+##     rownames(D) <- id
+##     colnames(D) <- id
 
-        INPATH.up() # update id with known paths
 
-        path.lengths <- sapply(CURRENTPATHDIST, sum)# compute path length
-        toKeep <- which.min(path.lengths) # keep best (shortest) one
-        res <- list(id=CURRENTPATH[[toKeep]], d=CURRENTPATHDIST[[toKeep]]) # return best path
-        return(res)
-    }
+##     ## AUXILIARY FUNCTIONS ##
+##     ## choose a starting sequence, as recent as possible
+##     chooseStartPoint <- function(){
+##         return(STARTPOINTS[which.max(seq.dates[STARTPOINTS])])
+##     }
 
 
-} # end seqTrack
+##     ## find k closest ancestors: returns a list(id=[id of the ancestors], d=[distances to ancestors])
+##     findAncestors <- function(currentPoint){
+##         candidates <- id[seq.dates < seq.dates[currentPoint]]
+##         if(length(candidates)==0) return(NULL) # this will indicate the end of the path
+##         res <- D[currentPoint,candidates]
+##         if(length(res) <= k) return(list(id=as.integer(names(res)), d=res)) # if less than k candidates
+##         res <- sort(res)[1:k] # if more than k, take k closest candidates
+##         return(list(id=as.integer(names(res)), d=res))
+##     }
+
+
+##     ## discard paths stemming from the worst ancestors a step ...
+##     discardPath <- function(step){
+##         tempDist <- sapply(CURRENTPATHDIST, sum)
+##         pathFac <- sapply(CURRENTPATH, function(e) e[step])
+##         pathLengths <- tapply(tempDist, pathFac, sum)
+##         toKeep <- names(pathLengths)[which.min(temp)]
+##         toKeep <- as.integer(toKeep)
+##         CURRENTPATH <<- CURRENTPATH[pathFac==toKeep]
+##         CURRENTPATHDIST <<- CURRENTPATHDIST[pathFac==toKeep]
+##         return()
+##     } # end discardPath
+
+
+##     ## update id in INPATH
+##     INPATH.up <- function(){
+##         temp <- unique(unlist(listPaths))
+##         INPATH <<- union(INPATH,temp)
+##         return()
+##     } # end INPATH.up
+
+
+##     ## retrieve an already known path
+##     findExistingPath <- function(id){
+##         if(id %in% INPATH){
+##             temp <- sapply(listPaths, function(e) id %in% e) # find in which path
+##             res <- listPaths[[which(temp)[1]]] # retrieve the path
+##             res <- res[which(res==id):length(res)] # cut the path
+##             return(res)
+##         } else return(NULL)
+##     } # end findExistingPath
+
+
+##     ## expand one path
+##     expandOnePath <- function(onePath){ # onePath has a single path (a vector)
+##         newPoints <- findAncestors(onePath[length(onePath)]) # get new points
+##         if(is.null(newPoints)) return(onePath)
+##         res <- lapply(1:length(newPoints), function(i) c(onePath, newPoints)) # duplicate path and add new pts
+##         return(res)
+##     } # end expandOnePath
+
+
+##     ## check if CURRENTPATH should stop (TRUE if yes, FALSE otherwise)
+##     checkEndCurrentPath <- function(){
+##         isEnded <- function(onePath){
+##             oldest <- onePath[length(onePath)]
+##             temp <- setdiff(id,oldest)
+##             if(all(seq.dates[temp] < seq.dates[oldest])) return(TRUE)
+##             return(FALSE)
+##         }
+
+##         temp <- sapply(CURRENTPATH, checkEndCurrentPath)
+##         if(all(temp)) return(TRUE) # if all paths are ended
+##         return(FALSE)
+##     } # end checkEndCurrentPath
+
+
+
+##     ## FIND ONE PATH ##
+##     findPath <- function(id){
+##         temp <- findExistingPath(id) # search for an already existing path
+##         if(!is.null(temp)) return(temp) # return existing path if needed
+
+##         ## WHILE LOOP ##
+##         keepSearching <- TRUE
+##         while(keepSearching){
+##             for(i in length(CURRENTPATH)){
+##                 temp <- CURRENPATH[[1]]
+
+##             }
+
+
+
+##             temp <- findExistingPath(id) # search for an already existing path
+##             if(!is.null(temp))
+##             keepSearching <- !checkEndCurrentPath() # stop searching if all paths are ended
+
+
+
+
+
+##         }
+
+##         INPATH.up() # update id with known paths
+
+##         path.lengths <- sapply(CURRENTPATHDIST, sum)# compute path length
+##         toKeep <- which.min(path.lengths) # keep best (shortest) one
+##         res <- list(id=CURRENTPATH[[toKeep]], d=CURRENTPATHDIST[[toKeep]]) # return best path
+##         return(res)
+##     }
+
+
+## } # end seqTrack

Modified: pkg/TODO
===================================================================
--- pkg/TODO	2009-04-30 21:05:38 UTC (rev 304)
+++ pkg/TODO	2009-05-08 14:49:57 UTC (rev 305)
@@ -30,8 +30,8 @@
 # CODE ISSUES:
 ==============
 * fix request 1.2-2.04 (implement adjusted heretozygosity in summary)
+* use packageStartupMessage() instead of cat in .First.lib
 
-
 # DOCUMENTATION ISSUES:
 =======================
 * explain new changes inside the tutorial (handling of AFLP/RAPD...)



More information about the adegenet-commits mailing list