[Rcpp-commits] r4296 - in pkg/RcppArmadillo: . inst/include/RcppArmadilloExtensions inst/unitTests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Mar 29 14:37:51 CET 2013


Author: edd
Date: 2013-03-29 14:37:50 +0100 (Fri, 29 Mar 2013)
New Revision: 4296

Modified:
   pkg/RcppArmadillo/ChangeLog
   pkg/RcppArmadillo/inst/include/RcppArmadilloExtensions/sample.h
   pkg/RcppArmadillo/inst/unitTests/runit.sample.R
Log:
applied another patch by Christian Gunning [mailed to rcpp-devel] for sample()


Modified: pkg/RcppArmadillo/ChangeLog
===================================================================
--- pkg/RcppArmadillo/ChangeLog	2013-03-29 07:33:00 UTC (rev 4295)
+++ pkg/RcppArmadillo/ChangeLog	2013-03-29 13:37:50 UTC (rev 4296)
@@ -1,3 +1,9 @@
+2013-03-29  Dirk Eddelbuettel  <edd at debian.org>
+
+	* inst/include/RcppArmadilloExtensions/sample.h: Applied another
+	patch by Christian Gunning
+	* inst/unitTests/runit.sample.R (test.sample): Idem
+
 2013-03-12  Dirk Eddelbuettel  <edd at debian.org>
 
 	* DESCRIPTION: Release 0.3.800.1
@@ -5,8 +11,8 @@
 	* inst/include/*: Upgraded to new release 3.800.1 of Armadillo
 
 	* inst/unitTests/cpp/armadillo.cpp: Factored out of
-	runit.RcppArmadillo.R to accelerate unit test run 
-	* inst/unitTests/runit.RcppArmadillo.R: Deploy refactored cpp code 
+	runit.RcppArmadillo.R to accelerate unit test run
+	* inst/unitTests/runit.RcppArmadillo.R: Deploy refactored cpp code
 
 2013-03-11  Dirk Eddelbuettel  <edd at debian.org>
 

Modified: pkg/RcppArmadillo/inst/include/RcppArmadilloExtensions/sample.h
===================================================================
--- pkg/RcppArmadillo/inst/include/RcppArmadilloExtensions/sample.h	2013-03-29 07:33:00 UTC (rev 4295)
+++ pkg/RcppArmadillo/inst/include/RcppArmadilloExtensions/sample.h	2013-03-29 13:37:50 UTC (rev 4296)
@@ -63,8 +63,14 @@
                 // Copy the given probabilities into an arma vector
                 arma::vec prob(fixprob.begin(), fixprob.size());
                 if (replace) {
-                    // check for walker alias conditions here??
-                    ProbSampleReplace(index, nOrig, size, prob);
+                    // check for walker alias conditions 
+                    int walker_test = sum( (probsize * prob) > 0.1);
+                    if (walker_test < 200) { 
+                        ProbSampleReplace(index, nOrig, size, prob);
+                    } else {
+                        throw std::range_error( "Walker Alias method not implemented. R-core sample() is likely faster for this problem.");
+                        // WalkerProbSampleReplace(index, nOrig, size, prob);
+                    }
                 } else {
                     ProbSampleNoReplace(index, nOrig, size, prob);
                 }

Modified: pkg/RcppArmadillo/inst/unitTests/runit.sample.R
===================================================================
--- pkg/RcppArmadillo/inst/unitTests/runit.sample.R	2013-03-29 07:33:00 UTC (rev 4295)
+++ pkg/RcppArmadillo/inst/unitTests/runit.sample.R	2013-03-29 13:37:50 UTC (rev 4296)
@@ -23,16 +23,6 @@
 }
 
 test.sample <- function() {
-    ## Seed needs to be reset to compare R to C++
-    seed <- 441
-
-    ## Input vectors to sample
-    N <- 100
-
-    ## Number of samples
-    ## works for size == N?!
-    size <- N%/%2
-
     ## set up S3 dispatching,
     ## simplifies lapply(tests, ...) below
     csample <- function(x, ...) UseMethod("csample")
@@ -42,6 +32,14 @@
     csample.character <- csample_character
     csample.logical <- csample_logical
 
+    ## Seed needs to be reset to compare R to C++
+    seed <- 441
+    ## Input vectors to sample
+    N <- 100
+    ## Number of samples
+    ## works for size == N?!
+    size <- N%/%2
+
     ## all atomic vector classes except raw
     ## for each list element, check that sampling works
     ## with and without replacement, with and without prob
@@ -70,18 +68,18 @@
         set.seed(seed)
         r.yes.no <- sample(dat, size, replace=T)
         set.seed(seed)
-        r.no.yes <- sample(dat, size, replace=F, prob=1:N)
+        r.no.yes <- sample(dat, size, replace=F, prob=probs)
         set.seed(seed)
-        r.yes.yes <- sample(dat, size, replace=T, prob=1:N)
+        r.yes.yes <- sample(dat, size, replace=T, prob=probs)
         ## C
         set.seed(seed)
         c.no.no <- csample(dat, size, replace=F)
         set.seed(seed)
         c.yes.no <- csample(dat, size, replace=T)
         set.seed(seed)
-        c.no.yes <- csample(dat, size, replace=F, prob=1:N)
+        c.no.yes <- csample(dat, size, replace=F, prob=probs)
         set.seed(seed)
-        c.yes.yes <- csample(dat, size, replace=T, prob=1:N)
+        c.yes.yes <- csample(dat, size, replace=T, prob=probs)
         ## comparisons
         checkEquals(r.no.no, c.no.no, msg=sprintf("sample.%s.no_replace.no_prob",.class))
         checkEquals(r.yes.no, c.yes.no, msg=sprintf("sample.%s.replace.no_prob",.class))
@@ -89,4 +87,21 @@
         checkEquals(r.no.yes, c.no.yes, msg=sprintf("sample.%s.no_replace.prob",.class))
         checkEquals(r.yes.yes, c.yes.yes, msg=sprintf("sample.%s.replace.prob",.class))
     })
+    ## Walker Alias method test
+    ## With replacement, >200 "nonzero" probabilities
+    ## Not implemented, see below
+    walker.N <- 1e3
+    walker.sample <- (1:walker.N)/10
+    walker.probs <- rep(0.1, walker.N)
+    ## uncomment following 5 lines if/when walker alias method is implemented
+    #set.seed(seed)
+    #r.walker <- sample( walker.sample, walker.N, replace=T, prob=walker.probs)
+    #set.seed(seed)
+    #c.walker <- csample( walker.sample, walker.N, replace=T, prob=walker.probs)
+    #checkEquals(r.walker, c.walker, msg=sprintf("Walker Alias method test"))
+    ## Walker Alias method is not implemented.
+    ## For this problem (replace, >200 non-zero probs) R is much faster
+    ## So throw an error and refuse to proceed
+    walker.error <- try( csample( walker.sample, walker.N, replace=T, prob=walker.probs), TRUE)
+    checkEquals(inherits(walker.error, "try-error"), TRUE, msg=sprintf("Walker Alias method test"))
 }



More information about the Rcpp-commits mailing list