[Rcpp-devel] Question concerning use of Raster in C vi Rcpp plugin and inline

Rainer M Krug Rainer at krugs.de
Tue Nov 16 16:43:06 CET 2010


On 11/16/2010 03:16 PM, Douglas Bates wrote:
> On Tue, Nov 16, 2010 at 4:31 AM, Rainer M Krug <r.m.krug-Re5JQEeQqe8AvxtiuMwx3w at public.gmane.org> 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?

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/841dDXCDGMXqaGq2UG9VpUWMKQH7oEaQurus at public.gmane.org
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:      Rainer at krugs.de

Skype:      RMkrug
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: dispersal.org
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20101116/4b8c5707/attachment-0001.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: dispersal.R
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20101116/4b8c5707/attachment-0001.asc>


More information about the Rcpp-devel mailing list