10

I would like to learn to use Rcpp. I went through the docs on the package's CRAN website, but i feel working on a practical example (second practical, considering convolve3) would be more helpfull .

I propose the following code, from the the robustbase package because it's neither too long nor too short, uses a combination of R types and R functions and has one of those small arithmetic iterations that are too slow in R. How would you go about Rcpp-ing it ?

scaleTau2<-function (x, c1 = 4.5, c2 = 3, consistency = TRUE, mu.too = FALSE){
n <- length(x)
medx <- median(x)
x. <- abs(x - medx)
sigma0 <- median(x.)
mu <- if (c1 > 0) {
    x. <- x./(sigma0 * c1)
    w <- 1 - x. * x.
    w <- ((abs(w) + w)/2)^2
    sum(x * w)/sum(w)
}
else medx
x <- (x - mu)/sigma0
rho <- x^2
rho[rho > c2^2] <- c2^2
if (!identical(consistency, FALSE)) {
    Erho <- function(b) 2*((1-b^2)*pnorm(b)-b*dnorm(b)+b^2)-1
    Es2 <- function(c2) Erho(c2*qnorm(3/4))
    nEs2 <-ifelse(consistency == "finiteSample",n-2,n)*Es2(c2)
}
else nEs2 <- n
c(if (mu.too) mu, sigma0 * sqrt(sum(rho)/nEs2))
}

Please explain as much as you can.

EDIT It's really the idea of a step by step explanation of how you would go about converting a well written (and documented) R code (so at least the foundations are okay) unto an efficient implementation. The choice of the code is arguable a bit random but i think it reflects the arch-typical script on our blueprints (calls R functions that one doesn't want to translate, uses arithmetic loops....).

EDIT2 from the comments i realize this may actually be a big work to do in C++ (i didn't realize it when posting the code). In regard of this, using individual pieces as pedagogical devices is ok. I'll eventually parse the pieces together by editing the question.

chl
  • 50,972
  • 18
  • 205
  • 364
user603
  • 21,225
  • 3
  • 71
  • 135
  • 3
    I *strongly disagree* with your edited title and added/edited question. You are simply mistaken if you consider Rcpp to be a code compiler, or when asking us to rewrite code for you. – Dirk Eddelbuettel Oct 27 '10 at 14:41
  • @Dirk:> sure, what would a suggested title be (the older one?). I don't really care about this particular function I'm interested in learning ways to make my codes run faster. If you have another example, please post it. I'll happily close this one. – user603 Oct 27 '10 at 14:49

1 Answers1

14

Interesting question, but quite possibly too challenging to be discussed briefly:

  • You would need a C++-side implementations of median()

  • The cited code from package robustbase is highly 'R-optimised' which may not be the best starting point.

  • Rcpp is not an 'R compiler' that you toss any such function at to 'make it faster'. It is more about connecting existing C++ code, or writing new C++ code.

  • Of course the above can be translated (Turing-equivalence and all that) but that may not be the best way to learn about using Rcpp. I think we have simpler examples on the mailing list.

Lastly, isn't this a programming question for SO? ;-)

Dirk Eddelbuettel
  • 8,362
  • 2
  • 28
  • 43
  • @Dirk:> a) there are plenty of implementation of median() (say 'pull' in package pcaPP) so it's fair game. b) you mean one won't notice a sizable increase in running times ? c) okay, but i think the issue with this code is not really the translation to C++, rather the idea of calling some R functions [pnorm, dnorm,...] in C++ (of course i can be really wrong) d) can you provide the link to your mailing list ? – user603 Oct 27 '10 at 13:21
  • Can we please split the sub-questions off one by one? A) you can call R function from C++ -- for convenience but not necessarily speed. See the examples/ in Rcpp. B) I said no such thing. C) That is a all easy since Rcpp 0.8.7, see the 'Rcpp sugar' docs, posts on Rcpp-devel and our recent presentations. D) It hangs off the R-forge page; just google for 'rcpp-devel'. – Dirk Eddelbuettel Oct 27 '10 at 13:25
  • @Dirk to A)B)C)D->thanks i'll read more. Would you know of a directory/listings of examples of simple R functions together with there translation in C++ (i know there are C++ funciton in any library but i would prefer modern examples, that is, ones that use as much as possible the features of Rcpp). – user603 Oct 27 '10 at 13:36
  • Did you see examples/ within the Rcpp package? Also, some of the HPC tutorials are by now 'too old' for current Rcpp practice. Make sure you lick a recent one -- say from useR this summer. – Dirk Eddelbuettel Oct 27 '10 at 13:41
  • The HPC tutorial i saw was over a year ago (btw thanks for these and the rest of your contribs), so maybe there is more there. I saw the examples in '/Rcpp/inst/examples/functionCallback/' but that's really not much. I realize there may be more spread around on your mailing list (and on your HPC slides) but that's too sparse. Would you have learned regression if it were taught that way (couple of contrived examples on slides and ML)? – user603 Oct 27 '10 at 13:56
  • 1
    1) Start at http://dirk.eddelbuettel.com/presentations.html and work your way down. 2) There are six subdirectories to examples/ so I am unsure why you focus on one. 3) There are 770+ unit tests that double as examples if you care to look closely enough. 4) There are eight (8) vignettes in the Rcpp package. 5) We authored a few other packages that use Rcpp, you could look at those too. 6) Lastly, CRAN lists fifteen packages depending upon Rcpp -- these all are examples too. – Dirk Eddelbuettel Oct 27 '10 at 14:00
  • 1
    Dude: There is a mailing list for the project you are interested in. *All* our documentation suggests to ask on the mailing list. So why-oh-why do you keep piling on here? Can we **please** stop that now. Lastly, your 'too superficial' would require some backing up. I will gladly review patches, just don't post them *here*. Ok? – Dirk Eddelbuettel Oct 27 '10 at 14:38
  • 2
    @kwak: Responding to "It's something that should be outsourced to the community": I look forward to seeing your contributions as you work through these examples yourself. – Joshua Ulrich Oct 27 '10 at 14:45
  • @Joshua Ultrich:> Sure: my question was my initial contribution to the process. Let's see if we can keep this ball rolling. – user603 Oct 27 '10 at 14:52