Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_r_kernel_local -p 8888:8888 -d palmoreck/jupyterlab_r_kernel:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_r_kernel_localDocumentación de la imagen de docker palmoreck/jupyterlab_r_kernel:1.1.0 en liga.
Esta nota utiliza métodos vistos en 1.5.Integracion_numerica
Instalación de microbenchmark:
In [2]:
install.packages("microbenchmark",lib="/usr/local/lib/R/site-library/",
repos="https://cran.itam.mx/",verbose=TRUE)
Documentación de Rcpp:
Rcpp permite conectar C++ a R de forma sencilla al utilizar la API de Rcpp.
¿Por qué usar Rcpp?
Aunque C o C++ requieren más líneas de código, son órdenes de magnitud más rápidos que R. Sacrificamos las ventajas que tiene R como rapidez en programación por velocidad en ejecución.
¿Cuando podríamos usar Rcpp?
En loops que no pueden vectorizarse de forma sencilla. Si tenemos loops en los que una iteración depende de la anterior.
Si hay que llamar una función millones de veces.
Si después de hacer perfilamiento y optimización de código no llegamos a nuestro tiempo objetivo.
Por qué no usamos C?
Sí es posible llamar funciones de C desde R pero resulta en más trabajo por parte de l@s programador@s. Por ejemplo, de acuerdo a H. Wickham:
...R’s C API. Unfortunately this API is not well documented. I’d recommend starting with my notes at R’s C interface. After that, read “The R API” in “Writing R Extensions”. A number of exported functions are not documented, so you’ll also need to read the R source code to figure out the details.
Y como primer acercamiento a la compilación de código desde R es preferible seguir las recomendaciones de H. Wickham en utilizar la API de Rcpp.
También se utiliza el paquete microbenchmark para medir tiempos de forma exacta:
Un microbenchmark es la medición del performance de un bloque pequeño de código. El paquete de R con el mismo nombre devolverá el tiempo medido en miliseconds (ms), microseconds ($\mu s$) o nanoseconds (ns) para el bloque de código dado y se repetirá ésta medición un número definido de veces. Las diferencias al correr varias veces la función de microbenchmark pueden deberse a varias razones tan simples como tener otras tareas corriendo en tu computadora.
En lo que sigue se utiliza el método del rectángulo para aproximar la integral definida de una función.
In [3]:
library(Rcpp)
library(microbenchmark)
La regla del rectángulo en código de R y utilizando vapply (vapply es más rápido que sapply pues se especifica con anterioridad el tipo de output que devuelve) es la siguiente:
In [7]:
Rcf1<-function(f,a,b,n){
#Compute numerical approximation using rectangle or mid-point method in
#an interval.
#Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
#Args:
# f (function): function of integrand
# a (int): left point of interval
# b (int): right point of interval
# n (int): number of subintervals
#Returns:
# Rcf (float)
h_hat<-(b-a)/n
sum_res<-0
x<-vapply(0:(n-1),function(j)a+(j+1/2.0)*h_hat,numeric(1))
for(j in 1:n){
sum_res<-sum_res+f(x[j])
}
h_hat*sum_res
}
In [8]:
f<-function(x)exp(-x^2)
Probaremos esta implementación Rcf1 básica para medir su tiempo de ejecución:
In [9]:
n<-10**6
aprox<-Rcf1(f,0,1,n)
In [10]:
aprox
Recuérdese revisar el error relativo:
In [11]:
err_relativo<-function(aprox,obj)abs(aprox-obj)/abs(obj)
In [12]:
obj<-integrate(Vectorize(f),0,1) #en la documentación de integrate
#se menciona que se utilice Vectorize
In [13]:
err_relativo(aprox,obj$value)
In [14]:
system.time(Rcf1(f,0,1,n))
Una implementación que utiliza la función sum de R es la siguiente:
In [15]:
Rcf2<-function(f,a,b,n){
#Compute numerical approximation using rectangle or mid-point method in
#an interval.
#Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
#Args:
# f (function): function of integrand
# a (int): left point of interval
# b (int): right point of interval
# n (int): number of subintervals
#Returns:
# Rcf (float)
h_hat<-(b-a)/n
x<-vapply(0:(n-1),function(j)a+(j+1/2.0)*h_hat,numeric(1))
h_hat*sum(f(x))
}
In [16]:
aprox<-Rcf2(f,0,1,n)
In [17]:
aprox
In [18]:
err_relativo(aprox,obj$value)
In [19]:
system.time(Rcf2(f,0,1,n))
y se redujo el tiempo de cálculo.
En Rcpp se tiene la función cppFunction que recibe código escrito en C++ para definir una función que puede ser utilizada desde R. Antes de usar tal función, reescribamos la regla del rectángulo de modo que no se utilice vapply:
In [21]:
Rcf3<-function(f,a,b,n){
#Compute numerical approximation using rectangle or mid-point method in
#an interval.
#Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
#Args:
# f (function): function of integrand
# a (int): left point of interval
# b (int): right point of interval
# n (int): number of subintervals
#Returns:
# Rcf (float)
h_hat<-(b-a)/n
sum_res<-0
for(i in 0:(n-1)){
x<-a+(i+1/2.0)*h_hat
sum_res<-sum_res+f(x)
}
h_hat*sum_res
}
In [23]:
n<-10**6
aprox<-Rcf3(f,0,1,n)
In [24]:
aprox
In [25]:
err_relativo(aprox,obj$value)
In [26]:
system.time(Rcf4(f,0,1,n))
Entonces se define el source code escrito en C++ que será el primer parámetro que recibirá cppFunction:
In [28]:
f_str<-'double Rcf_Rcpp(double a, double b, int n){
double h_hat;
double sum_res=0;
int i;
double x;
h_hat=(b-a)/n;
for(i=0;i<=n-1;i++){
x = a+(i+1/2.0)*h_hat;
sum_res+=exp(-pow(x,2));
}
return h_hat*sum_res;
}'
In [29]:
cppFunction(f_str)
Si queremos obtener más información de la ejecución de la línea anterior podemos usar:
In [31]:
cppFunction(f_str, verbose=TRUE, rebuild=TRUE) #también usamos rebuild=TRUE
#para que se vuelva a compilar,
#ligar con la librería en C++
#y todo lo que realiza cppFunction
#detrás del telón
Comentarios:
Al ejecutar la línea de cppFunction, Rcpp compilará el código de C++ y construirá una función de R que se conecta con la función compilada de C++ (este se le llama wrapper).
Si se observa en la salida de arriba se verá que hay un bloque de C y un tipo de dato SEXP que de acuerdo a H. Wickham:
...functions that talk to R must use the SEXP type for both inputs and outputs. SEXP, short for S expression, is the C struct used to represent every type of object in R. A C function typically starts by converting SEXPs to atomic C objects, and ends by converting C objects back to a SEXP. (The R API is designed so that these conversions often don’t require copying.)
Revisemos el tiempo de esta función:
In [32]:
aprox_rcpp<-Rcf_Rcpp(0,1,n)
In [33]:
err_relativo(aprox_rcpp,obj$value)
In [34]:
system.time(Rcf_Rcpp(0,1,n))
Y utilizando microbenchmark:
In [35]:
mbk<-microbenchmark(
Rcf1(f,0,1,n),
Rcf2(f,0,1,n),
Rcf3(f,0,1,n),
Rcf_Rcpp(0,1,n),
times=10
)
In [36]:
print(mbk)
Se observa que la función compilada Rcf_Rcpp es dos órdenes de magnitud más rápida que Rcf1 y un orden de magnitud más rápida que Rcf2 y Rcf3.
NumericVector
En Rcpp se tienen clases que se relacionan con los tipos de dato en R para vectores. Entre éstas se encuentran NumericVector, IntegerVector, CharacterVector y LogicalVector que se relacionan con vectores tipo numeric, integer, character y logical. Por ejemplo, para el caso de NumericVector se tiene:
In [37]:
f_str <-'NumericVector el(NumericVector x){
return exp(log(x));
}'
In [38]:
cppFunction(f_str)
In [40]:
print(el(seq(0,1,by=.1)))
Para el caso de la regla de integración del rectángulo, podríamos pensar en R en una implementación como la siguiente:
In [41]:
Rcf3b<-function(f,a,b,n){
#Compute numerical approximation using rectangle or mid-point method in
#an interval.
#Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
#Args:
# f (function): function of integrand
# a (int): left point of interval
# b (int): right point of interval
# n (int): number of subintervals
#Returns:
# Rcf (float)
h_hat<-(b-a)/n
fx<-f(vapply(0:(n-1),function(j)a+(j+1/2.0)*h_hat,numeric(1))) #evaluate f
h_hat*sum(fx)
}
In [42]:
aprox<-Rcf3b(f,0,1,n)
In [43]:
err_relativo(aprox,obj$value)
In [44]:
system.time(Rcf3(f,0,1,n))
Y para poner un ejemplo de NumericVector para esta regla, podemos primero calcular los nodos y evaluar f en ellos:
In [46]:
a<-0
b<-1
h_hat<-(b-a)/n
In [47]:
fx<-f(vapply(0:(n-1),function(j)a+(j+1/2.0)*h_hat,numeric(1)))
In [49]:
print(tail(fx))
In [50]:
f_str<-'
double Rcf_Rcpp2(NumericVector f_x,double h_hat){
double sum_res=0;
int i;
int n = f_x.size();
for(i=0;i<=n-1;i++){
sum_res+=f_x[i];
}
return h_hat*sum_res;
}'
In [51]:
cppFunction(f_str,rebuild=TRUE)
In [52]:
system.time(Rcf_Rcpp2(fx,h_hat))
Revisamos el error relativo:
In [53]:
aprox_rcpp2<-Rcf_Rcpp2(fx,h_hat)
In [54]:
err_relativo(aprox_rcpp2,obj$value)
Y constrastamos con microbenchmark:
In [55]:
mbk<-microbenchmark(
Rcf1(f,0,1,n),
Rcf2(f,0,1,n),
Rcf3(f,0,1,n),
Rcf3b(f,0,1,n),
Rcf_Rcpp(0,1,n),
Rcf_Rcpp2(fx,h_hat),
times=10
)
In [56]:
print(mbk)
Comentarios:
Obsérvese que está utilizando el método .size() que regresa un integer.
No estamos midiendo en condiciones iguales pues las otras funciones construían los nodos... por ejemplo es súper rápida la ejecución de Rcf_Rcpp2 y no tanto la siguiente:
In [57]:
system.time(fx<-f(vapply(0:(n-1),function(j)a+(j+1/2.0)*h_hat,numeric(1))))
Entonces debimos de haber medido como:
In [58]:
mbk<-microbenchmark(
Rcf1(f,0,1,n),
Rcf2(f,0,1,n),
Rcf3(f,0,1,n),
Rcf3b(f,0,1,n),
Rcf_Rcpp(0,1,n),
f(vapply(0:(n-1),function(j)a+(j+1/2.0)*h_hat,numeric(1))),
times=10
)
In [59]:
print(mbk)
NumericVector por ejemplo para crear los nodos:
In [60]:
f_str<-'NumericVector Nodos(double a, double b, int n){
double h_hat=(b-a)/n;
int i;
NumericVector x(n);
for(i=0;i<n;i++)
x[i]=a+(i+1/2.0)*h_hat;
return x;
}'
In [61]:
cppFunction(f_str,rebuild=TRUE)
In [66]:
print(Nodos(0,1,2))
También en Rcpp es posible llamar funciones definidas en el ambiente global, por ejemplo:
In [67]:
f_str='RObject fun(double x){
Environment env = Environment::global_env();
Function f=env["f"];
return f(x);
}
'
In [68]:
cppFunction(f_str,rebuild=TRUE)
In [69]:
fun(1)
In [70]:
f(1)
In [71]:
fun
.Call es una función base para llamar funciones de C desde R:.
There are two ways to call C functions from R: .C() and .Call(). .C() is a quick and dirty way to call an C function that doesn’t know anything about R because .C() automatically converts between R vectors and the corresponding C types. .Call() is more flexible, but more work: your C function needs to use the R API to convert its inputs to standard C data types. H. Wickham.
In [72]:
f
Ejercicios:
Revisar rcpp-sugar y Rcpp syntactic sugar y proponer programas que utilicen sugar.
Reescribe las reglas del trapecio y de Simpson de la nota 1.5.Integracion_numerica con Rcpp.
Referencias
Referencias para tutoriales de C++:
Otras referencias de Rcpp:
Rcpp Note (no es documentación oficial).