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

[Rainer M Krug]

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

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


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