UP | HOME |

Creación de una biblioteca de enlace dinámico con R

Creación de una biblioteca de enlace dinámico con R

Sin necesidad de desarrollar un paquete completo, construyamos una biblioteca de enlace dinámico con R implementando una alternativa a la función expand.grid.

Las hilanderas o la fábula de Aracne, Diego Rodríguez de Silva y Velázquez.

El problema de las permutaciones con repetición

La función expand.grid crea una tabla de datos (data.frame) con todas las combinaciones posibles de los elementos de una lista. Para una tarea necesitamos generar una retícula grid en formato matriz con unas determinadas características: debe obtener todas las permutaciones con repetición de un conjunto x de números enteros o reales (distintos o no), donde cada elemento se repite hasta n veces. Además requiere incluir unas columnas adicionales para almacenar posteriormente información sobre cada permutación.

El siguiente ejemplo detalla la forma de proceder con expand.grid. Calcula las permutaciones con repetición de dos elementos tomados de tres en tres, obteniendo ocho permutaciones, y añade otras cuatro columnas. 

x <- c(1:2)
n <- 3L
ntotal <- 7L

listainicial <- vector("list", ntotal)
for(i in 1:ntotal){
    if(i <= n) {
        listainicial[[i]]  <- x
    } else {
        listainicial[[i]]  <- NA
    }
}
mpointsR <-as.matrix(expand.grid(listainicial, KEEP.OUT.ATTRS = FALSE))
dimnames(mpointsR) <- NULL
print(mpointsR)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    1    1   NA   NA   NA   NA
[2,]    2    1    1   NA   NA   NA   NA
[3,]    1    2    1   NA   NA   NA   NA
[4,]    2    2    1   NA   NA   NA   NA
[5,]    1    1    2   NA   NA   NA   NA
[6,]    2    1    2   NA   NA   NA   NA
[7,]    1    2    2   NA   NA   NA   NA
[8,]    2    2    2   NA   NA   NA   NA

Los parámetros x, n y ntotal varían y cuando alcanzan determinado tamaño, la función expand.grid ralentiza considerablemente todo el proceso de computación.

Compilación de una biblioteca de enlace dinámico

La solución propuesta para superar esta limitación consiste en crear una versión especializada de expand.grid, sin necesidad de considerar todas las casuísticas que incluye esta función. Primero se hizo una versión simplificada en lenguaje R, y como aún así el tiempo de ejecución también sobrepasaba lo razonable, se optó por diseñar una función en C/C++.

Para ello descargamos el siguiente material:

  • El fichero paquete.h, que incluye las cabeceras usuales.
  • El fichero expandgrid.cpp, que contiene una versión especializada de la función expand.grid.
  • El fichero test.r, que compara los resultados de ambos procedimientos.
  • El fichero Makefile automatiza el proceso: crea una biblioteca de enlace dinámico y ejecuta el fichero test.r cuya salida se recoge en el test.r.Rout. El fichero Makefile debe borrarse si se incluye esta carpeta en el src de un paquete. En ese caso se utiliza el fichero Makevars.

Descomprimimos el material, y en ese directorio ejecutamos la instrucción para crear una biblioteca dinámica. 

$ R CMD SHLIB expandgrid.cpp -o paquete.so

Si todo funciona correctamente se crea una biblioteca llamada paquete.so en Linux (en Windows utilice el nombre paquete.dll).

El fichero test.r carga la biblioteca, ejecuta la función expand_grid_cpp definida en paquete.so y comprueba que los resultados son iguales con ambos procedimientos. 

## Load library. But first, we have to unload it if it is loaded.
## Windows: paquete.dll
## Unix: paquete.so
sharedlibrary <- paste0("paquete",.Platform$dynlib.ext)
tryCatch({
    dyn.unload(file.path(".",sharedlibrary)); 
}, warning = function(w) {
    cat("Library is not loaded.\n")
}, error = function(e) {
    cat("Library is not loaded.\n")
})
dyn.load(file.path(".",sharedlibrary), local = TRUE, now = TRUE) 

mpointsCpp <- .Call("expand_grid_cpp",x,  n, ntotal)
print(mpointsCpp)
class(mpointsCpp[1])

all.equal(mpointsR,mpointsCpp) # TRUE

Rendimiento

Si observamos el tiempo de ejecución, la función expand.grid tarda veinte veces más que la diseñada ad hoc. 

library(microbenchmark)
microbenchmark(
    cpp=.Call("expand_grid_cpp",x,  n, ntotal),
    R=expand.grid(listainicial, KEEP.OUT.ATTRS = FALSE)
)


Unit: nanoseconds
 expr   min      lq     mean  median      uq   max neval
  cpp   735   831.0  1446.15   917.0  1918.5 11941   100
    R 19065 19577.5 20709.61 20016.5 20706.5 60480   100

Lea el manual de tortura

Si ha llegado hasta aquí, debe avanzar un poco más y leer el manual Writing R Extensions, en particular la sección de problemas de memoria. Dado que supone una tortura, aquí puede encontrar una forma de proceder con valgrind, que hemos de aplicar siempre a nuestras bibliotecas.

Coda

Este ejemplo detalla cómo crear una biblioteca de enlace dinámico y cómo se utiliza en un programa de R sin necesidad de disponer de un paquete completo. Sirve para desarrollar poco a poco las funciones necesarias.

Referencias

Código en C/C++ de la versión alternativa a expand.grid

Aquí se presenta el código del fichero expandgrid.cpp y el de la cabecera paquete.h.

El código de expandgrid.cpp consta de estas instrucciones: 

#include "paquete.h"

/*

Create a grid of length(x)^n, with the values of x
and add (ntotal - n) columns.
Return a matrix with dimensions (length(x)^n , ntotal)
Particular case of expand.grid

x <- 1:3
n <- 4L
ntotal <- 7L

mpoints <- .Call("expand_grid_cpp",x,  n, ntotal)

listainicial <- vector("list", ntotal)
for(i in 1:ntotal){
  if(i <= n) {
 listainicial[[i]]  <- x
                     } else {
 listainicial[[i]]  <- NA
                     }
 }
mpointsR <-as.matrix(expand.grid(listainicial, KEEP.OUT.ATTRS = FALSE))
dimnames(mpointsR) <- NULL

all.equal(mpointsR,mpoints) # TRUE

head(mpoints)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    1    1    1   NA   NA   NA
[2,]    2    1    1    1   NA   NA   NA
[3,]    3    1    1    1   NA   NA   NA
[4,]    1    2    1    1   NA   NA   NA
[5,]    2    2    1    1   NA   NA   NA
[6,]    3    2    1    1   NA   NA   NA
...

*/


extern "C" SEXP expand_grid_cpp(SEXP x, SEXP n, SEXP ntotal){

  // X is a vector of integer or real numbers.
  // N and NTOTAL are scalar integers such that 0< N <= NTOTAL.
  // Return a matrix of the permutations with n repetitions of x,
  // and add (ntotal -n ) NA columns.

  if (!Rf_isVectorAtomic(x) || !LENGTH(x))
    Rf_error("X must be an atomic vector with length greater than 0 (%d).", LENGTH(x));

  if(!Rf_isInteger(n) || LENGTH(n) != 1)
    Rf_error("N is either not an integer value ('%s') or its length is not 1 (%d).", Rf_type2char(TYPEOF(n)), LENGTH(n));

  if(!Rf_isInteger(ntotal) || LENGTH(ntotal) != 1)
    Rf_error("NTOTAL is either not an integer value ('%s') or its length is not 1 (%d).", Rf_type2char(TYPEOF(ntotal)), LENGTH(n));

  const R_len_t nrepetitions = INTEGER(n)[0];
  const R_len_t ncolumns = INTEGER(ntotal)[0];

  if( nrepetitions < 1 || ncolumns < nrepetitions || nrepetitions == NA_INTEGER || ncolumns == NA_INTEGER )
    Rf_error("N is '%d', lesser than 1 or greater than the number of columns NCOL (%d).", nrepetitions,ncolumns);

  // Number of combinations. It can become a very, very, very big magnitude.
  const R_len_t lx=LENGTH(x);
  const R_len_t nrows = pow(lx, nrepetitions);

  // BE CAREFUL! ISNAN(nrows) does not work here!
  if(nrows < 0 )
    Rf_error("The number of combinations (%.0f) is too big! ('nrows'). Check it within R: Can you create a 'matrix(NA, nrow= %d^%d, ncol= %d)'?",pow(lx,nrepetitions),  lx, nrepetitions, ncolumns);

  // We create the matrix. If it is too big, it gives the error 'Not enough memory for your dreams'.
  int protecti=0;
  SEXP result = PROTECT(Rf_allocMatrix(TYPEOF(x), nrows , ncolumns)); protecti++;

  // vcount is the count/value for each column.
  // vcount[jcol] = 0,1,2,..., lx-1
  R_len_t *vcount = (R_len_t*) (R_alloc(nrepetitions,sizeof(R_len_t)));
  for(R_len_t i = 0; i < nrepetitions ; i++)  vcount[i] = 0;

  const R_len_t maxcount = lx -1; // The maximum count that can appear in a column
  const R_len_t lastcolumn  = nrepetitions -1; // Last column of the repetition submatrix.
  const R_len_t firstcolumn = 0; // First column.

  switch(TYPEOF(x))
    {
    case LGLSXP:
    case INTSXP:
      {

        const int* px = INTEGER(x);
        int* presult = INTEGER(result);

        // for each row
        for(R_len_t irow = 0; irow < nrows ; irow++){

          // for the columns to fill with values of x
          for(R_len_t jcol = 0; jcol < nrepetitions ; jcol++){

            // When we are at the first column,
            // we decide that a count must vary
            // when the first column gets an excessive count.
            if(vcount[firstcolumn] > maxcount && jcol == firstcolumn ) {

              // Numk is the first position/column which its count/value can be increased.
              // Therefore, the previous columns have the maximum count/value.
              // So we reset the count/values for the previous columns
              // and increase the count/values of the position Numk.
              R_len_t numk = 0;
              for( numk =0;  numk < nrepetitions && vcount[numk] >= maxcount; numk++);
              for(R_len_t k =0;  k < numk; k++){
                vcount[k]=0;
              }

              // If Numk is greater than the last column, we have finished the grid:
              // All the columns have the maximum count/value.
              // Otherwise, increase the value of the column numk
              if(numk < nrepetitions)  vcount[numk]++; // numk <= lastcolumn

            }

            // Fill this position with the value of X that corresponds to
            // the count of this column.
            presult[ IDX(irow,jcol,nrows) ] = px[ vcount[jcol]  ] ;

            // When it is the last column, the count of the first column always increases.
            // It is the signal to reset some columns.
            if(jcol == lastcolumn)  vcount[firstcolumn]++;

          } // End of columns of permutations

          // Fill the others columns with NA
          for(R_len_t jcol = nrepetitions; jcol < ncolumns ; jcol++){
            presult[ IDX(irow,jcol,nrows) ] = NA_INTEGER;
          }
        }

        break;
      }
    case REALSXP:
      {

        const double* px = REAL(x);
        double* presult = REAL(result);

        for(R_len_t irow = 0; irow < nrows ; irow++){

          for(R_len_t jcol = 0; jcol < nrepetitions ; jcol++){

            if(vcount[firstcolumn] > maxcount && jcol == firstcolumn ) {

              int numk = 0;
              for( numk =0;  numk < nrepetitions && vcount[numk] >= maxcount; numk++);
              for(R_len_t k =0;  k < numk; k++){
                vcount[k]=0;
              }

              if(numk < nrepetitions){
                vcount[numk ]++;
              }

            }

            presult[ IDX(irow,jcol,nrows) ] = px[ vcount[jcol]  ] ;
            if(jcol == lastcolumn){
              vcount[firstcolumn]++;
            }
          }

          for(R_len_t jcol = nrepetitions; jcol < ncolumns ; jcol++){
            presult[ IDX(irow,jcol,nrows) ] = NA_REAL;
          }
        }
        break;
      }

    default:
      UNPROTECT(protecti);
      Rf_error("Unsupported type for X ('%s')",Rf_type2char(TYPEOF(x)));
      return R_NilValue;
    }

  UNPROTECT(protecti);
  return result;

}



// ============================================================
// ============================================================

El código de la cabecera paquete.h consta de estas instrucciones: 

#ifndef PAQUETE_H
#define PAQUETE_H


#include <R.h>
#define R_NO_REMAP
#include <Rinternals.h>
#include <Rdefines.h>

#ifdef _OPENMP
#include <omp.h>
#include <pthread.h>
#endif

// ============================================================

// Matrix[i,j] = Matrix[ i + j*dim(0) ] i=0,1,2,..,dim0-1; j=0,1,2,3,..
#define IDX(i, j, dim0) ( (i) + (j) * (dim0) )

// ============================================================

// expandgrid.cpp
extern "C" SEXP expand_grid_cpp(SEXP x, SEXP n, SEXP ntotal);

// ============================================================

#endif