[adegenet-commits] r349 - pkg/R

[noreply at r-forge.r-project.org]

Author: jombart
Date: 2009-06-03 17:53:08 +0200 (Wed, 03 Jun 2009)
New Revision: 349

Modified:
   pkg/R/seqTrack.R
Log:
Trying to retaining only a percentage of results. 


Modified: pkg/R/seqTrack.R
===================================================================
--- pkg/R/seqTrack.R	2009-06-03 14:05:43 UTC (rev 348)
+++ pkg/R/seqTrack.R	2009-06-03 15:53:08 UTC (rev 349)
@@ -389,7 +389,7 @@
 ## 3) uncomment, adapt, and test code for missing data
 ##
 optimize.seqTrack <- function(nstep=10, step.size=1e3,
-                              seq.names, seq.dates, W, thres=NULL, optim=c("min","max"),
+                              seq.names, seq.dates, W, thres=0.1, optim=c("min","max"),
                               prox.mat=NULL, mu0, seq.length,
                               rMissDate=.rUnifTimeSeq, ...){
 
@@ -441,11 +441,11 @@
     }
 
 
-    ## SET THRESHOLD IF NEEDED ##
-    if(is.null(thres)){
-        thres <- sum(seqTrack(seq.names=seq.names, seq.dates=seq.dates, W=W,
-                                    optim=optim, prox.mat=prox.mat, ...)$weight, na.rm=TRUE)
-    }
+    ## SET THRESHOLD IF NEEDED ## ## NO LONGER USED
+    ##  if(is.null(thres)){
+    ##         thres <- sum(seqTrack(seq.names=seq.names, seq.dates=seq.dates, W=W,
+    ##                                     optim=optim, prox.mat=prox.mat, ...)$weight, na.rm=TRUE)
+    ##     }
 
 
     ## AUXILIARY FUNCTIONS ##
@@ -490,57 +490,61 @@
 
                 res.new <- seqTrack(seq.names=seq.names, seq.dates=myDates, W=W,
                                     optim=optim, prox.mat=prox.mat, ...)
-                temp <- val.res(res.new)
-                if(ifelse(optim=="min", temp < thres, temp > thres)){
-                    ances <- cbind(ances, res.new$ances)
-                    date <- cbind(date, as.character(res.new$date))
-                    ances.date <- cbind(ances.date, as.character(res.new$ances.date))
-                    valRes <- c(valRes, temp)
-                }
+
+                ##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
 
-            ## 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)
-                ##                 }
+            ## retain thres% of the dates ##
+            toKeep <- valRes < quantile(valRes, thres) ## NOT WORKING FOR optim==max !!!
+            date <- date[toKeep]
+            newDates <- apply(date, 1, function(vec)
+                              sample(vec, size=step.size, replace=TRUE))
+            newDates <- t(newDates)
 
+            ## re-initialize posterior distributions
+            if(i<nstep){
+                ances <- integer(0)
+                date <- character(0)
+                ances.date <- character(0)
+                valRes <- numeric(0)
+            } # end if
+        } # end for i
 
-                ## new dates
-                ## DEBUGING ##
-                temp <- apply(date,1,function(vec) length(unique(vec)))
-                cat("\n i =", i, "j =", j, "Number of dates per sequence:\n")
-                print(temp)
-                cat("\nHead of date:\n")
-                print(head(date))
-                ## cat("\nProba vector:\n")
-                ## print(w)
-                ## END DEBUGING ##
+        ##  ## 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)
-                                  sample(vec, size=step.size, replace=TRUE))
-                newDates <- t(newDates)
 
-                ## re-initialize posterior distributions
-                if(i<nstep){
-                    ances <- integer(0)
-                    date <- character(0)
-                    ances.date <- character(0)
-                    valRes <- numeric(0)
-                }
-            }
-        }
-    }
+        ## new dates
+        ## DEBUGING ##
+        ## temp <- apply(date,1,function(vec) length(unique(vec)))
+        ##                 cat("\n i =", i, "j =", j, "Number of dates per sequence:\n")
+        ##                 print(temp)
+        ##                 cat("\nHead of date:\n")
+        ##                 print(head(date))
+        ## cat("\nProba vector:\n")
+        ## print(w)
+        ## END DEBUGING ##
 
+        ## newDates <- apply(date, 1, function(vec) #  used a weighted sampling
+        ##                                  sample(vec, size=step.size, replace=TRUE, prob=w))
 
+
+    } # end if(!any(isMissDate))
+
+
     ##  ## OTHER CASE: HANDLE MISSING DATES
     ##     if(any(isMissDate)){