<jupyter_start><jupyter_text>**Notas para contenedor de docker:** Comando de docker para ejecución de la nota de forma local:nota: cambiar `` por la ruta de directorio que se desea mapear a `/datos` dentro del contenedor de docker.```docker run --rm -v :/datos --name jupyterlab_r_kernel_local -p 8888:8888 -d palmoreck/jupyterlab_r_kernel:1.1.0```password para jupyterlab: `<PASSWORD>`Detener el contenedor de docker:```docker stop jupyterlab_r_kernel_local``` Documentación de la imagen de docker `palmoreck/jupyterlab_r_kernel:1.1.0` en [liga](https://github.com/palmoreck/dockerfiles/tree/master/jupyterlab/r_kernel). --- Esta nota utiliza métodos vistos en [1.5.Integracion_numerica](https://github.com/ITAM-DS/analisis-numerico-computo-cientifico/blob/master/temas/I.computo_cientifico/1.5.Integracion_numerica.ipynb) **Instalación de microbenchmark:**<jupyter_code>install.packages("microbenchmark",lib="/usr/local/lib/R/site-library/",
                repos="https://cran.itam.mx/",verbose=TRUE)<jupyter_output>system (cmd0): /usr/lib/R/bin/R CMD INSTALL

foundpkgs: microbenchmark, /tmp/Rtmpu90LwL/downloaded_packages/microbenchmark_1.4-7.tar.gz

files: /tmp/Rtmpu90LwL/downloaded_packages/microbenchmark_1.4-7.tar.gz

1): succeeded '/usr/lib/R/bin/R CMD INSTALL -l '/usr/local/lib/R/site-library' /tmp/Rtmpu90LwL/downloaded_packages/microbenchmark_1.4-7.tar.gz'<jupyter_text>Rcpp Documentación de Rcpp:* [rcpp por <NAME>](http://dirk.eddelbuettel.com/code/rcpp.html)* [rcpp.org](http://www.rcpp.org/) **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 <NAME>:*...R’s C API. Unfortunately this API is not well documented. I’d recommend starting with my notes at [R’s C interface](http://adv-r.had.co.nz/C-interface.html). After that, read “[The R API](http://cran.rstudio.com/doc/manuals/r-devel/R-exts.htmlThe-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](https://github.com/wch/r-source) to figure out the details.*Y como primer acercamiento a la compilación de código desde `R` es preferible seguir las recomendaciones de <NAME> en utilizar la API de `Rcpp`. **También se utiliza el paquete [microbenchmark](https://www.rdocumentation.org/packages/microbenchmark/versions/1.4-7/topics/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.**<jupyter_code>library(Rcpp)
library(microbenchmark)<jupyter_output><empty_output><jupyter_text>La regla del rectángulo en código de `R` y utilizando [vapply](https://www.rdocumentation.org/packages/functools/versions/0.2.0/topics/Vapply) (`vapply` es más rápido que `sapply` pues se especifica con anterioridad el tipo de `output` que devuelve) es la siguiente:<jupyter_code>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
}
f<-function(x)exp(-x^2)<jupyter_output><empty_output><jupyter_text>Probaremos esta implementación `Rcf1` básica para medir su tiempo de ejecución:<jupyter_code>n<-10**6
aprox<-Rcf1(f,0,1,n)
aprox<jupyter_output><empty_output><jupyter_text>**Recuérdese** revisar el error relativo:<jupyter_code>err_relativo<-function(aprox,obj)abs(aprox-obj)/abs(obj)
obj<-integrate(Vectorize(f),0,1) #en la documentación de integrate
                                 #se menciona que se utilice Vectorize
err_relativo(aprox,obj$value)
system.time(Rcf1(f,0,1,n))<jupyter_output><empty_output><jupyter_text>Una implementación que utiliza la función `sum` de `R` es la siguiente:<jupyter_code>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))
}
aprox<-Rcf2(f,0,1,n)
aprox
err_relativo(aprox,obj$value)
system.time(Rcf2(f,0,1,n))<jupyter_output><empty_output><jupyter_text>y se redujo el tiempo de cálculo.  Hacia la compilación con Rcpp En `Rcpp` se tiene la función [cppFunction](https://www.rdocumentation.org/packages/Rcpp/versions/1.0.3/topics/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`:<jupyter_code>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
}
n<-10**6
aprox<-Rcf3(f,0,1,n)
aprox
err_relativo(aprox,obj$value)
system.time(Rcf4(f,0,1,n))<jupyter_output><empty_output><jupyter_text>Entonces se define el `source code` escrito en `C++` que será el primer parámetro que recibirá `cppFunction`:<jupyter_code>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;
        }'
cppFunction(f_str)<jupyter_output><empty_output><jupyter_text>Si queremos obtener más información de la ejecución de la línea anterior podemos usar:<jupyter_code>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<jupyter_output>Generated code for function definition: 
--------------------------------------------------------

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
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;
        }

Generated extern "C" functions 
--------------------------------------------------------


#include <Rcpp.h>
// Rcf_Rcpp
double Rcf_Rcpp(double a, double b, int n);
RcppExport SEXP sourceCpp_1_Rcf_Rcpp(SEXP aSEXP, SEXP bSEXP, SEXP nSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< double >::type a(aSEXP);
    Rcpp::traits::input_parameter< double >::type b(bSEXP);
    Rcpp::traits::input_parameter< int[...]<jupyter_text>**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 <NAME>:*...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:<jupyter_code>aprox_rcpp<-Rcf_Rcpp(0,1,n)
err_relativo(aprox_rcpp,obj$value)
system.time(Rcf_Rcpp(0,1,n))<jupyter_output><empty_output><jupyter_text>Y utilizando `microbenchmark`:<jupyter_code>mbk<-microbenchmark(
    Rcf1(f,0,1,n),
    Rcf2(f,0,1,n),
    Rcf3(f,0,1,n),
    Rcf_Rcpp(0,1,n), 
    times=10
    )
print(mbk)<jupyter_output>Unit: milliseconds
              expr       min         lq       mean     median         uq
  Rcf1(f, 0, 1, n) 1134.9580 1143.60462 1286.19387 1207.79665 1217.92329
  Rcf2(f, 0, 1, n)  668.5134  687.85826  746.14824  731.34272  751.19857
  Rcf3(f, 0, 1, n)  524.9488  536.67870  545.12018  539.86892  552.13084
 Rcf_Rcpp(0, 1, n)   16.5403   17.25606   18.64566   17.97957   19.13085
        max neval
 1723.17939    10
  947.80422    10
  577.35055    10
   24.64291    10<jupyter_text>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:<jupyter_code>f_str <-'NumericVector el(NumericVector x){
        return exp(log(x));
}'
cppFunction(f_str)
print(el(seq(0,1,by=.1)))<jupyter_output>[1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0<jupyter_text>Para el caso de la regla de integración del rectángulo, podríamos pensar en `R` en una implementación como la siguiente:<jupyter_code>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)
}
aprox<-Rcf3b(f,0,1,n)
err_relativo(aprox,obj$value)
system.time(Rcf3(f,0,1,n))<jupyter_output><empty_output><jupyter_text>Y para poner un ejemplo de `NumericVector` para esta regla, podemos primero calcular los nodos y evaluar `f` en ellos:<jupyter_code>a<-0
b<-1
h_hat<-(b-a)/n
fx<-f(vapply(0:(n-1),function(j)a+(j+1/2.0)*h_hat,numeric(1)))
print(tail(fx))
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;
        }'
cppFunction(f_str,rebuild=TRUE)
system.time(Rcf_Rcpp2(fx,h_hat))<jupyter_output><empty_output><jupyter_text>**Revisamos** el error relativo:<jupyter_code>aprox_rcpp2<-Rcf_Rcpp2(fx,h_hat)
err_relativo(aprox_rcpp2,obj$value)<jupyter_output><empty_output><jupyter_text>Y constrastamos con `microbenchmark`:<jupyter_code>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
    )
print(mbk)<jupyter_output>Unit: milliseconds
                 expr         min          lq        mean      median
     Rcf1(f, 0, 1, n) 1192.918665 1228.678904 1315.967484 1264.626690
     Rcf2(f, 0, 1, n)  708.988752  721.018438  838.991386  791.609931
     Rcf3(f, 0, 1, n)  533.528447  557.654910  642.007659  599.152741
    Rcf3b(f, 0, 1, n)  688.495578  723.240941  840.585161  743.023979
    Rcf_Rcpp(0, 1, n)   16.944433   17.587898   21.350258   21.209751
 Rcf_Rcpp2(fx, h_hat)    1.047825    1.074875    1.348535    1.126084
          uq         max neval
 1414.261255 1489.786964    10
  855.935157 1190.395839    10
  690.248859  850.867292    10
  942.462597 1213.450117    10
   24.679616   29.429521    10
    1.200255    3.288781    10<jupyter_text>**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:<jupyter_code>system.time(fx<-f(vapply(0:(n-1),function(j)a+(j+1/2.0)*h_hat,numeric(1))))<jupyter_output><empty_output><jupyter_text>Entonces debimos de haber medido como:<jupyter_code>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
    )
print(mbk)<jupyter_output>Unit: milliseconds
                                                                expr       min
                                                    Rcf1(f, 0, 1, n) 1118.9874
                                                    Rcf2(f, 0, 1, n)  661.8217
                                                    Rcf3(f, 0, 1, n)  518.4573
                                                   Rcf3b(f, 0, 1, n)  657.6112
                                                   Rcf_Rcpp(0, 1, n)   16.7331
 f(vapply(0:(n - 1), function(j) a + (j + 1/2) * h_hat, numeric(1)))  708.2641
        lq       mean     median         uq        max neval
 1148.5022 1221.05506 1177.96496 1294.26572 1396.64626    10
  669.4143  712.31013  681.45650  695.29276 1009.91352    10
  531.4959  596.47586  557.93353  652.30480  778.27262    10
  683.1871  735.76753  686.61225  727.13774 1014.09308    10
   17.0036   18.01895   18.27387   18.49182   19.71365    10
  744.7894  824.46560  758.49632  923.38135 1010.55658    10<jupyter_text>* También se pueden devolver vectores de tipo `NumericVector` por ejemplo para crear los nodos:<jupyter_code>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;
}'
cppFunction(f_str,rebuild=TRUE)
print(Nodos(0,1,2))<jupyter_output>[1] 0.25 0.75<jupyter_text>**También en `Rcpp` es posible llamar funciones definidas en el ambiente global, por ejemplo:**<jupyter_code>f_str='RObject fun(double x){
Environment env = Environment::global_env();
Function f=env["f"];
return f(x);
}
'
cppFunction(f_str,rebuild=TRUE)
fun(1)
f(1)
fun<jupyter_output><empty_output><jupyter_text>**.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.* **<NAME>**.<jupyter_code>f<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12914143213395085, 'fraction_numerical': 0.06679729248307802, 'mean_word_length': 3.1271747120803726, 'pattern_counts': {'":': 0, '<': 132, '<?xml version=': 0, '>': 88, 'https://': 7, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '337450', 'license': 'Apache-2.0', 'n_tokens_mistral': 6984, 'n_tokens_neox': 6295, 'n_words': 1598, 'path': 'temas/I.computo_cientifico/1.7.Reescribir_funciones_a_C++_Rcpp.ipynb', 'repo_name': 'Danahirmt/analisis-numerico-computo-cientifico', 'size': '16865'}