On Tue, Nov 16, 2010 at 9:43 AM, Rainer M Krug <[email protected]> wrote: > On 11/16/2010 03:16 PM, Douglas Bates wrote: >> On Tue, Nov 16, 2010 at 4:31 AM, Rainer M Krug >> <[email protected]> wrote: >> Hi >> >> I am new to C, and my question might be basic, but I am struggling with it: >> >> I have about 10 lines of code in a simulation, which take up about 50% >> of the whole simulation, As the simulation, takes several hours to run, >> I want to make the code faster. After trying different approaches in R, >> with different degrees of success, I decided to try out inline and Rcpp >> for this. >> >> I did some progress (using sugar), and the net-result is so far, that it >> is as fast / slow as R; but as I said, this is a first step and I am >> happy so far. My problem now is, in the R code, I have the following >> (the complete code is below): >> >> s <- seeds[x:(x+dx2), y:(y+dy2)] >> dispSeeds[x,y] <- sum( >> rbinom( >> s, >> s, >> sd2D >> ), >> na.rm = TRUE >> ) >> >>> I think it unlikely that you will be able to write C/C++ code that >>> runs much faster than that expression in R. The sum and rbinom >>> functions are vectorized and other than a small interpretation >>> overhead not much can be saved by going to compiled code. > > The main reason why I am trying to use C are the loops around the code > above. See attached code example. > >> >> As the rbinom part is only using the vector representation of s, I have >> already translated it into C. >> >> So my problem / question is: can I use a similar syntax for s in inline >> using sugar? >> >> As the code above sits in two large loops, I would like to translate the >> whole loops into C. >> >> Below the code (R and R / inline): >> >> maxx, maxy, dx, dy: scalar integers >> seeds: integer matrix >> dispSeeds: integer matrix >> >> The pure R code: >> >> #+BEGIN_SRC R >> dispRparFor <- function( >> maxx = maxx, >> maxy = maxy, >> seeds = seeds, >> sd2D = sd2D, >> dx = dx, >> dy = dy, >> dispSeeds = dispSeeds # the return value >> ) { >> dx2 <- 2*dx >> dy2 <- 2*dy >> for ( y in 1:maxy ) { >> cat(y, " ") >> for (x in 1:maxx) { >> s <- seeds[x:(x+dx2), y:(y+dy2)] >> if (all(is.na(s))) { >> dispSeeds[x,y] <- NA >> } else { >> dispSeeds[x,y] <- sum( >> rbinom( >> s, >> s, >> sd2D >> ), >> na.rm = TRUE >> ) >> } >> } >> } >> return(dispSeeds) >> } >> #+END_SRC >> >> >> How far I got with using inline and Rcpp / sugar: >> >> >> #+BEGIN_SRC R >> dispRparForSugar <- function( >> maxx = maxx, >> maxy = maxy, >> seeds = seeds, >> sd2D = sd2D, >> dx = dx, >> dy = dy, >> dispSeeds = dispSeeds # the return value >> ) { >> library(inline) >> library(Rcpp) >> dx2 <- 2*dx >> dy2 <- 2*dy >> >> >> fx2 <- cxxfunction( >> sig = signature(N = "integer", SIZE = "integer", >> PROB = "double"), >> body = ' >> Rcpp::IntegerVector n (N); >> Rcpp::IntegerVector size (SIZE); >> Rcpp::NumericVector prob (PROB); >> >> int res; >> res = 0; >> >> for( int i=0; i<n.size(); i++){ >> if (size[i]>0) { >> res += rbinom( 1, size[i], >> prob[i] )[0]; >> } >> } >> >> return wrap( res ); >> ', >> plugin = "Rcpp", >> verbose = TRUE >> ) >> >> >> >> for ( y in 1:maxy ) { >> cat(y, " ") >> for (x in 1:maxx) { >> s <- seeds[x:(x+dx2), y:(y+dy2)] >> if (all(is.na(s))) { >> dispSeeds[x,y] <- NA >> } else { >> dispSeeds[x,y] <- fx2(s, s, sd2D) >> } >> } >> } >> return(dispSeeds) >> } >> #+END_SRC >> >>> It is peculiar to call cxxfunction within another function if the code >>> being compiled is fixed. cxxfunction is a function generator and one >>> usually regards a call to cxxfunction as the equivalent of entering a >>> function definition. > > I changed that - yes, it was peculiar but based on my experimentation. > >> >>> It doesn't appear that you are using n in your C++ code (other than >>> taking its length) and your C++ code seems to be evaluating >>> sum(rbinom(1, s, sd2D)) rather than sum(rbinom(s, s, sd2D)). Is that >>> intentional? > > Left over - changed as well. Thanks. > >> >>> And if your fx2 is always to be called as fx2(s, s, sd2D) then why not >>> just pass s and sd2D? > > Done (partly). > > > To add some more strange things: > > I attach the code in a working example as an R script file which will > show something strange (for me): > > the function overwrites a variable, namely sd2D. Am I doing something > wrong or do I underestimate the dangers of using Rcpp and inline and C > in general?
Yes, it will. You are assigning the NumericVector's s and sd2D to use the storage of the SD2D argument so when you assign s[indS] in the first loop you are overwriting the contents of SD2D. I am testing a rewrite now. > In addition, the code is not doing what it is supposed to be doing, but > that is a different question at the moment. > > Cheers, > > Rainer > >> >> >> >> Ultimately, I would like to have the complete dispRparFor function in C. >> >> >> For the ones interested, I could send a working example with data. >> >> Cheers, >> >> Rainer >> > _______________________________________________ > Rcpp-devel mailing list > Rcpp-devel-Z+qqJ2/[email protected] > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel >>> > > > -- > Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation > Biology, UCT), Dipl. Phys. (Germany) > > Centre of Excellence for Invasion Biology > Natural Sciences Building > Office Suite 2039 > Stellenbosch University > Main Campus, Merriman Avenue > Stellenbosch > South Africa > > Tel: +33 - (0)9 53 10 27 44 > Cell: +27 - (0)8 39 47 90 42 > Fax (SA): +27 - (0)8 65 16 27 82 > Fax (D) : +49 - (0)3 21 21 25 22 44 > Fax (FR): +33 - (0)9 58 10 27 44 > email: [email protected] > > Skype: RMkrug > > _______________________________________________ > Rcpp-devel mailing list > [email protected] > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel > > _______________________________________________ Rcpp-devel mailing list [email protected] https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
