R/C functions in R code

From Wiki1

< R
Jump to: navigation, search

This is a quick way to write a loop as a C function in R. This description assumes you are using Windows. A few small changes are required for other OSs.

First, install Rtools

The following R code illustrates a C function, its compilation and its inclusion in an R function.

'benchmark' is used to compare the speed of of the raw R function, the compiled R function, the C function using a 'double' accumulator, and the C function using a 'long double' accumulator.

The speed advantage of using the C function compared with the raw R function is in the order of 100. The compiled R function runs roughly twice as fast as the raw function.

install.packages('rbenchmark')
library(rbenchmark)
library(compiler)

# Speeding up a function with a loop

Rfun <- function(t, N = 100000) {
  tmp1 <- 0
  tmp2 <- 0
  tmp3 <- 0
  for (j in N:1) {
    tmp1 <- tmp1 + 0.5*log(j*(j+1)/(j*(j+1)-2*t))
    tmp2 <- tmp2 + 1/(j*(j+1)-2*t)
    tmp3 <- tmp3 + 2/((j*(j+1)-2*t))^2
  }
  c(tmp1, tmp2, tmp3)
}

# Using the R compiler

Rfun.comp <- cmpfun(Rfun)

# The function in C

cat('
    #include <R.h>
    #include <math.h>
    void fun(double *t, double *N, double *ret)
{
    double tmp1, tmp2, tmp3, j;
    tmp1 = tmp2 = tmp3 = 0;
    for (j = *N; j > 0.9 ; j--) {
      tmp1 = tmp1 + 0.5*log(j*(j+1)/(j*(j+1)-2.0*t[0]));
      tmp2 = tmp2 + 1.0/(j*(j+1)-2.0*t[0]);
      tmp3 = tmp3 + 2.0/pow(((j*(j+1)-2.0*t[0])),2.0);
    };
    ret[0] = tmp1;
    ret[1] = tmp2;
    ret[2] = tmp3;
    ret[3] = j;
}
    ', file = 'fun.c')
system( 'rm fun.dll')
system( "R CMD SHLIB fun.c" ) 
dyn.load("fun.dll")

Cfun <- function(t, N = 100000){
  ret <- double(4)
  ret2 <- .C('fun', t = as.double(t), N = as.double(N),
             ret = ret)
  ret2$ret 
}  

# Using long double as an accumulator to see what effect that has

cat('
    #include <R.h>
    #include <math.h>
    void funl(double *t, double *N, double *ret)
{
    long double tmp1, tmp2, tmp3, j, Nl;
    tmp1 = tmp2 = tmp3 = 0;
    Nl = *N;
    for (j = Nl; j > 0.9 ; j--) {
      tmp1 = tmp1 + 0.5*log(j*(j+1)/(j*(j+1)-2.0*t[0]));
      tmp2 = tmp2 + 1.0/(j*(j+1)-2.0*t[0]);
      tmp3 = tmp3 + 2.0/pow(((j*(j+1)-2.0*t[0])),2.0);
    };
    ret[0] = tmp1;
    ret[1] = tmp2;
    ret[2] = tmp3;
    ret[3] = j;
}
    ', file = 'funl.c')
system( 'rm funl.dll')
system( "R CMD SHLIB funl.c" ) 
dyn.load("funl.dll")

Cfunl <- function(t, N = 100000){
  ret <- double(4)
  ret2 <- .C('funl', t = as.double(t), N = as.double(N),
               ret = ret)
    ret2$ret 
}  

# -475 is a plausible argument for this function
Rfun( -475) - Rfun.comp( -475)
Cfun( -475) - Cfunl( -475)
Rfun( -475) - Cfunl(-475)[1:3]


benchmark( Rfun( -475), Rfun.comp(-475),Cfun(-475),Cfunl(-475), replications=10)
Personal tools