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

[Douglas Bates]

On Tue, Nov 16, 2010 at 4:31 AM, Rainer M Krug <r.m.krug at gmail.com> wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> 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.

> 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.

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?

And if your fx2 is always to be called as fx2(s, s, sd2D) then why not
just pass s and sd2D?


>
> 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
>
> - --
> 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
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.10 (GNU/Linux)
> Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
>
> iEYEARECAAYFAkziXZ0ACgkQoYgNqgF2egoL7wCfR5RapMcddSmM2f8ACBzjVG6p
> GaQAni+XQJs2wzcJQNXpMxwFY/FqRa/W
> =PmVZ
> -----END PGP SIGNATURE-----
> _______________________________________________
> 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
>