Reunion Usuarios R Madrid - 10 de diciembre de 2015

Paquetes con Clase

R y C++?

  • Generalmente C++ como "acelerador" de funciones de R: recodificación en un lenguaje más veloz
  • Lo realmente potente: Extensión de R e integración de código externo

En el mundo real

Interfaz al instante

Ingredientes:

  • Librería con objetos de C++ que queremos integrar
  • Función en R que use los objetos anteriores como queremos.
  • R, Rcpp y RStudio. (thx Hadley, thx Dirk)
  • (un poquito de C++ programming skills)

https://github.com/bichopaulinho/PaquetesConClase

Nuestro caso: Regresión Múltiple usando Gradient Descent

  • Integraremos el algoritmo de optimización Gradient Descent (implementación sencilla: librería Purple, ligeramente adaptada a nuestro problema)
  • Aplicación al problema de regresion lineal multiple
  • Resultado: función contenida en un paquete (GradDesc)
Rcpp::Rcpp.package.skeleton(name="GradDesc",path="/tmp")
## GradDesc
## ├── DESCRIPTION
## ├── GradDesc.Rproj
## ├── inst
## ├── man
## ├── NAMESPACE
## ├── R
## └── src
## 
## 4 directories, 3 files

Estructura del directorio

## GradDesc
## ├── DESCRIPTION
## ├── inst
## │   └── include
## │       ├── ObjectiveFunction
## │       ├── OptimizationAlgorithm
## │       └── Utilities
## ├── man
## ├── NAMESPACE
## ├── R
## │   ├── GradDesc.R
## │   └── RcppExports.R
## └── src
##     ├── GradientDescent.cpp
##     ├── Makevars
##     ├── Makevars.win
##     ├── ObjectiveFunction.cpp
##     ├── OptimizationAlgorithm.cpp
##     ├── RcppExports.cpp
##     ├── RcppGradientDescent.cpp
##     ├── RcppModuloSumaCuadrados.cpp
##     ├── Rcpp_purple_wrappers.cpp
##     └── SumaCuadrados.cpp
## 
## 8 directories, 14 files

Funcion objetivo. Esquema

RcppGradientDescent.cpp (Rcpp)

// [[Rcpp::export]]
Rcpp::NumericVector GradientDescentSumaCuadrados(SEXP A,SEXP b) {
    Purple::Vector<double> M_b = Rcpp::as<Purple::Vector<double>>(b);
    Purple::Matrix<double> M_A = Rcpp::as<Purple::Matrix<double>>(A);
    //....
    return(Rcpp::wrap<Purple::Vector<double>>(minimalArgument));
}
  • Argumentos de entrada: SEXP A,SEXP b. Objetos de R (Matriz y vector)
  • Se transforman en objetos de la clase Purple: Purple::Vector y Purple::Matrix: Especializacion de la plantilla de Rcpp::as. Se define en purple_wrappers.cpp

Funcion objetivo. Esquema (II)

GradientDescent.cpp (Rcpp)

// [[Rcpp::export]]
Rcpp::NumericVector GradientDescentSumaCuadrados(SEXP A,SEXP b) {

    Purple::Vector<double> M_b = Rcpp::as<Purple::Vector<double>>(b);
    Purple::Matrix<double> M_A = Rcpp::as<Purple::Matrix<double>>(A);
    
    //....
    return(Rcpp::wrap<Purple::Vector<double>>(minimalArgument));
}
  • Creación objeto función objetivo, algoritmo de optimización, parámetros, ejecución, etc.
  • Salida: parámetros optimizados (objeto R) mediante especialización a medida de Rcpp::wrap. Se define en purple_wrappers.cpp

Compilación

Fuentes de Purple incluidos en src/

  • Ventaja: el paquete es "autocontenido". Es la forma habitual de distribución.
  • Incluir todos los cabeceros de la clase en un directorio "include" (dentro del directorio del paquete)
  • Las clases de la librería se compilan al instalar el paquete.
  • Compilación del paquete con Build tools (devtools, RStudio). Utiliza fichero Makevars (Makevars.win):
## -*- mode: makefile; -*-

PKG_CXXFLAGS=-I../inst/include
## With R 3.1.0 or later, you can uncomment the following line to tell R to
## enable compilation with C++11 (where available)
CXX_STD = CXX11

Cómo funciona Rcpp:

RcppGradientDescent.cpp ([[Rcpp::export]]) -> src/RcppExports.cpp -> R/RcppExports.R (Función en R)

devtools::install(".", args="--no-multiarch")

La hora de la verdad

Estimación de coeficientes de regresión por Gradient Descent vs OLS

library(GradDesc)
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
A <- as.matrix(iris[,2:4])
y <- iris[,1]

La hora de la verdad (II)

Estimación de coeficientes de regresión por Gradient Descent vs OLS

coefs <- GradientDescentSumaCuadrados(A,y)
names(coefs) <- colnames(A)
# Usando función lm (quitamos el término independiente)
lm(data=iris[,-5], formula=Sepal.Length~.-1)$coefficients
##  Sepal.Width Petal.Length  Petal.Width 
##    1.1210617    0.9235289   -0.8956758
coefs
##  Sepal.Width Petal.Length  Petal.Width 
##    1.1210617    0.9235288   -0.8956758

La hora de la verdad (III)

Estimación de coeficientes de regresión por Gradient Descent vs OLS

# Con término independiente: incluimos una columna de unos en la matriz del modelo
A <- cbind(Interc=rep(1,nrow(A)), A)
coefs <- GradientDescentSumaCuadrados(A,y)
names(coefs) <- colnames(A)
lm(data=iris[,-5], formula=Sepal.Length~.)$coefficients
##  (Intercept)  Sepal.Width Petal.Length  Petal.Width 
##    1.8559975    0.6508372    0.7091320   -0.5564827
coefs
##       Interc  Sepal.Width Petal.Length  Petal.Width 
##    1.8559957    0.6508376    0.7091322   -0.5564831

Extra: RCPP_MODULE

  • Macro RCPP_MODULE: permite exportar objetos directamente. Similar a // [[Rcpp::export]]
  • Muy útil para exportar clases completas (constructores, métodos, etc)
  • En src/ModuloSumaCuadrados.cpp definimos la clase GradientDescentSumaCuadradosClase
class GradientDescentSumaCuadradosClase {
    public:
        GradientDescentSumaCuadradosClase(SEXP A, SEXP y){
            //...
        }
        SEXP getMinimalArgument(){
            return Rcpp::wrap<Purple::Vector<double>>(gradient_descent.getMinimalArgument());
        }
        void setInitialArgument(SEXP vect_inicial){
            gradient_descent.setInitialArgument(Rcpp::as<Purple::Vector<double>>(vect_inicial));
        }
        //...
    private:
        Purple::GradientDescent gradient_descent;
        //...
    
};

Extra: RCPP_MODULE (II)

  • La exportamos como módulo. En src/ModuloSumaCuadrados.cpp:
class GradientDescentSumaCuadradosClase{
    //...
};
RCPP_MODULE(ModuloSumaCuadrados){
    using namespace Rcpp;

    class_<GradientDescentSumaCuadradosClase>( "GradientDescentSumaCuadradosClase" )

    .constructor<SEXP,SEXP>()

    .method( "solve", &GradientDescentSumaCuadradosClase::getMinimalArgument )
    .method( "setGradientNormGoal", &GradientDescentSumaCuadradosClase::setGradientNormGoal )
    .method( "setEvaluationGoal", &GradientDescentSumaCuadradosClase::setEvaluationGoal )
    .method( "starting_point", &GradientDescentSumaCuadradosClase::setInitialArgument )
    ;
}
  • Hay que añadir en DESCRIPTION : RcppModules: ModuloSumaCuadrados
  • Hay que establecer un punto de carga en un archivo .R: R/GradDesc.R: loadModule("ModuloSumaCuadrados",TRUE)

Extra: RCPP_MODULE (III)

RCPP_MODULE crea una clase S4 de R.

obj = new(GradientDescentSumaCuadradosClase,A,y)

punto_inicial_alternativo = c(1,-1,-2,1)
obj$starting_point(punto_inicial_alternativo)

coefs2 = obj$solve()
names(coefs2)<-colnames(A)
coefs2
##       Interc  Sepal.Width Petal.Length  Petal.Width 
##    1.9291884    0.6326327    0.6991565   -0.5391704
lm(data=iris[,-5], formula=Sepal.Length~.)$coefficients
##  (Intercept)  Sepal.Width Petal.Length  Petal.Width 
##    1.8559975    0.6508372    0.7091320   -0.5564827

Extra: Otra alternativa de Compilación

Purple como librería externa

  • Ventaja: obtenemos la librería externa de su ubicación original (integración continua)
  • Se enlaza al compilar el paquete referenciandola en en el Makevars (PKG_LIBS)
## -*- mode: makefile; -*-
PKG_CXXFLAGS=-I../../Purple
PKG_LIBS=-L../../lib -lPurple
## With R 3.1.0 or later, you can uncomment the following line to tell R to
## enable compilation with C++11 (where available)
CXX_STD = CXX11
  • Previamente hay que compilar Purple. Es necesario? Recomendable para evitar problemas de compatibilidad con la tool-chain que compila los paquetes en R (sobre todo en Windows)
  • Con ayuda del fichero Makefile (Makefile.win): make -f Makefile
  • Disponible en VersionCompilacionExterna/

Problemas Frecuentes

  • A veces no recompila los archivos. Usar: devtools::clean_dll(); devtools::clean_source()
  • Borrar los archivos R/RcppExport.R y src/RcppExport.cpp
  • Si usamos roxygen2 revisar el NAMESPACE. Tiene que tener:
## # Generated by roxygen2: do not edit by hand
## 
## export(GradientDescentSumaCuadrados)
## export(GradientDescentSumaCuadradosClase)
## import(Rcpp)
## useDynLib(GradDesc)
  • Si tenemos algún memory leak: R -f script_que_peta.R -d valgrind ref
  • Si queremos debbugear con gdb: R --debugger=gdb ref
  • Pedir socorro: Stackoverflow, lista Rcpp-devel, lista R-es.

Referencias

MUCHAS GRACIAS!