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

Douglas Bates bates at stat.wisc.edu
Tue Nov 16 17:15:16 CET 2010


On Tue, Nov 16, 2010 at 9:43 AM, Rainer M Krug <Rainer at krugs.de> wrote:
> 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?

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/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
>
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
>


More information about the Rcpp-devel mailing list