8 PIEP y POA en R

En esta sección se hará una revisión de los elementos básicos de programación imperativa estructurada procedimental (PIEP) y de programación orientada a arreglos (POA) de una implementación en R (obviamente, R como lenguaje de programación y NO como herramienta de cálculo o de análisis de datos).

En esta revisión se hace una pequeña introducción a R, a sus elementos básicos, a sus principales estructuras de control, a la creación de subprogramas, y al manejo de sus tipos de datos compuestos.

Se espera que al finalizar las actividades de esta sección, el estudiante entienda y tenga clara la manera en que un algoritmo, diseñado bajo los paradigmas de programación imperativa estructurada procedimental y orientado a arreglos, se puede implementar mediante el uso del lenguaje de programación R.

Preparación de clase

Para la siguiente sección, lea todo y ejecute todo el código que allí se incluye, haciendo todas las pruebas, cambios y experimentos que se les puedan ocurrir sobre dicho código.

En sus propias palabras, explique lo que le transmitió y lo que le enseño cada parte de lo que leyó, ejecutó, probó y experimentó; incluya su discusión, reflexiones y conclusiones al respecto; exponga lo que no entendió e intente encontrar por su cuenta respuestas a las preguntas que le surgieron, para poder compartirlas en clase.

8.1 Elementos básicos de PIEP y POA en R

Cuaderno computacional en Google Colaboratory:
Elementos básicos, estructuras de control, creación de subprogramas y tipos de datos compuestos en R

8.2 Ejemplos

Ejemplo 8.1 Haga un adecuado análisis, diseño e implementación en R de un subprograma con una solución de PIEP y uno con una solución de POA (array-oriented) que reciban un número entero y que devuelvan el factorial del número recibido. Implemente una prueba rápida (una serie de instrucciones) que le permita probar los subprogramas requeridos para varios valores distintos. Por último, implemente una serie de instrucciones que permita realiar una comparación de tiempos entre las dos soluciones solicitadas y la función factorial() de R.

Análisis

El factorial de un número entero positivo k, denotado k! se puede definir como el producto de todos los números enteros positivos menores o iguales que k: k! = (2)(3) \cdots (k-2)(k-1)(k) o lo que es lo mismo, k! = (k)(k-1)(k-2) \cdots (3)(2)

Además, es necesario tener en cuenta que, 0! = 1 y es buena idea tener en cuenta que, 1! = 1 \qquad 2! = 2

Implementación en R (solución de PIEP)

# Solamente sirve para un parámetro vector con un solo elemento:
mi_fact_klength1_PIEP <- function(klength1){
  if(length(klength1) == 1){
    if(klength1 > 2){
      res <- klength1
      while(klength1 > 2){
        klength1 <- klength1 - 1
        res <- res * klength1
      }
      return(res)
    }else{
      if(klength1 >= 0){
        return(ifelse(klength1 == 0, 1, klength1))
      }
    }
  }
  return(NaN)
}

# Factorial elemento a elemento de un vector dado:
mi_fact_PIEP <- function(k){
  n <- length(k)
  res <- NULL
  for(i in 1:n){
    res[i] = mi_fact_klength1_PIEP(k[i])
  }
  return(res)
}

Implementación en R (solución de POA)

# Solamente sirve para un parámetro vector con un solo elemento:
mi_fact_klength1_POA <- function(klength1){
  if(length(klength1)==1){
    if(klength1 > 2){
      return(prod(2:klength1))
    }else{
      if(klength1 >= 0){
        return(ifelse(klength1 == 0, 1, klength1))
      }
    }
  }
  return(NaN)
}

# `Vectorize()` toma una función y devuelve una función capaz de devolver 
# el mismo resultado pero elemento a elemento de un vector dado:
mi_fact_POA_Vectorize <- Vectorize(mi_fact_klength1_POA)

# Implementación propia (con la ayuda de `sapply()`) capaz de devolver el 
# factorial elemento a elemento de un vector dado:
mi_fact_POA_sapply <- function(k){
  sapply(k, mi_fact_klength1_POA)
}

Prueba rápida

v <- c(-14, -8, 0, 1, 2, 8, 14)
for(i in v){
  cat("(", sprintf(i, fmt='%3.0f'), ")!\t", 
      sprintf(mi_fact_klength1_PIEP(i), fmt='%12.0f'), "\t", 
      sprintf(mi_fact_klength1_POA(i), fmt='%12.0f'), "\t", 
      sprintf(factorial(i), fmt='%12.0f'), "\n", sep = "")
}

(-14)!           NaN             NaN             NaN
( -8)!           NaN             NaN             NaN
(  0)!             1               1               1
(  1)!             1               1               1
(  2)!             2               2               2
(  8)!         40320           40320           40320
( 14)!   87178291200     87178291200     87178291200

v <- c(-14, -8, 0, 1, 2, 8, 14)
resul <- c(mi_fact_PIEP(v), mi_fact_POA_Vectorize(v), 
           mi_fact_POA_sapply(v), factorial(v))
resul <- matrix(resul, ncol = 4) 
colnames(resul) <- c("mi_fact_PIEP(v)", "mi_fact_POA_Vectorize(v)", 
                     "mi_fact_POA_sapply(v)", "factorial(v)")
rownames(resul) <- paste0("(", v, ")!")
resul

       mi_fact_PIEP(v) mi_fact_POA_Vectorize(v) mi_fact_POA_sapply(v)
(-14)!             NaN                      NaN                   NaN
(-8)!              NaN                      NaN                   NaN
(0)!                 1                        1                     1
(1)!                 1                        1                     1
(2)!                 2                        2                     2
(8)!             40320                    40320                 40320
(14)!      87178291200              87178291200           87178291200
       factorial(v)
(-14)!          NaN
(-8)!           NaN
(0)!              1
(1)!              1
(2)!              2
(8)!          40320
(14)!   87178291200

Comparación de tiempos

# `require()` intenta cargar la librería, si no la puede cargar devuelve FALSE
if (!require(microbenchmark)){
  install.packages("microbenchmark") # Instala librería
  library(microbenchmark) # Carga librería
}

Loading required package: microbenchmark

set.seed(10310506)
k <- sample(8:14, 1)
cat("k =", k)

k = 13

mbm <- microbenchmark("mi_fact_klen1_PIEP"={mi_fact_klength1_PIEP(k)},
                      "mi_fact_klen1_POA"={mi_fact_klength1_POA(k)},
                      "mi_fact_PIEP"={mi_fact_PIEP(k)},
                      "mi_fact_POA_Vectorize"={mi_fact_POA_Vectorize(k)},
                      "mi_fact_POA_sapply"={mi_fact_POA_sapply(k)},
                      "factorial"={factorial(k)},
                      times=1e3)
mbm # Tabla

Unit: nanoseconds
                  expr   min      lq      mean  median      uq     max neval
    mi_fact_klen1_PIEP   943  1121.0  1405.383  1245.5  1412.0   11591  1000
     mi_fact_klen1_POA   747   978.0  1241.781  1118.5  1281.0   11097  1000
          mi_fact_PIEP  1877  2590.0  3387.142  3172.0  3537.5   13772  1000
 mi_fact_POA_Vectorize 31051 32824.5 39474.619 33529.5 34621.0 2351249  1000
    mi_fact_POA_sapply 11767 13025.0 16296.346 13661.5 14343.5 1337000  1000
             factorial   371   510.0   658.587   583.0   681.5   11373  1000

boxplot(mbm, bty="n") # Gráfico boxplot

if (!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}

Loading required package: ggplot2

autoplot(mbm) # gráfico parecido al boxplot

# `require()` intenta cargar la librería, si no la puede cargar devuelve FALSE
if (!require(microbenchmark)){
  install.packages("microbenchmark") # Instala librería
  library(microbenchmark) # Carga librería
}
set.seed(10310506)
k <- sample(0:15, 25, replace=TRUE)
print(k)

 [1]  5  0  6 14 15  7 10 13  6  4 15 13  4  2 15 13 13  5 12  8  7  2  2  9  1

mbm <- microbenchmark("mi_fact_PIEP"={mi_fact_PIEP(k)},
                      "mi_fact_POA_Vectorize"={mi_fact_POA_Vectorize(k)},
                      "mi_fact_POA_sapply"={mi_fact_POA_sapply(k)},
                      "factorial"={factorial(k)},
                      times=1e3)
mbm # Tabla

Unit: microseconds
                  expr    min      lq      mean  median      uq      max neval
          mi_fact_PIEP 32.816 39.4255 43.749542 42.0235 45.0275   83.321  1000
 mi_fact_POA_Vectorize 67.148 73.7655 87.369809 76.0095 85.8530 3072.639  1000
    mi_fact_POA_sapply 43.002 48.4200 53.645399 50.4640 55.6535  134.417  1000
             factorial  1.585  1.9185  2.231503  2.1180  2.2895   21.043  1000

if (!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}
autoplot(mbm) # gráfico parecido al boxplot

Ejemplo 8.2 Haga un adecuado análisis, diseño e implementación en R de un subprograma con una solución de PIEP y uno con una solución de POA (array-oriented) que devuelvan una aproximación de la función coseno evaluada en un número real dado entre 0 y \frac{\pi}{4}. Los subprogramas solicitados deben tener en cuenta que:

\sum\limits_{k = 0}^{n} \frac{(-1)^k}{(2k)!} c^{2k} \xrightarrow[n \to \infty]{} \cos(c)
En esta ocasión, fije la cantidad de términos a usar de la serie. Es decir, primero analice y determine, a lo sumo, cuántos términos son mayores que 1 \times 10^{-15} para c entre 0 y \frac{\pi}{4}, y luego en la implementación use siempre esa misma cantidad de términos.

Adicionalmente, implemente una prueba rápida (una serie de instrucciones) que le permita probar los subprogramas requeridos para varios valores distintos. Por último, implemente una serie de instrucciones que permita realizar una comparación de tiempos entre las dos soluciones solicitadas y la función cos() de R.

Análisis

Para c entre 0 y \frac{\pi}{4}, \frac{1}{(2k)!} c^{2k} < \frac{1}{(2k)!} \left(\frac{\pi}{4}\right)^{2k} o sea que si se encuentra k_0 tal que para todo k > k_0 se tiene que, \frac{1}{(2k)!} \left(\frac{\pi}{4}\right)^{2k} < 1 \times 10^{-15} entonces \frac{1}{(2k)!} c^{2k} < 1 \times 10^{-15} para todo k > k_0 y para todo c entre 0 y \frac{\pi}{4}.

Reescribiendo,

\begin{align*} \frac{1}{(2k)!} \left(\frac{\pi}{4}\right)^{2k} &< 1 \times 10^{-15} \\ \left(\frac{\pi}{4}\right)^{2k} &< \frac{(2k)!}{1 \times 10^{15}} \\ \left(\frac{\pi}{4}\right)^{2k} - \frac{(2k)!}{1 \times 10^{15}} &< 0 \end{align*}

Encontrando numéricamente una raíz para la función f(k) = \left(\frac{\pi}{4}\right)^{2k} - \frac{(2k)!}{1 \times 10^{15}}, concluimos que \frac{1}{(2k)!} \left(\frac{\pi}{4}\right)^{2k} < 1 \times 10^{-15} para todo k > 8.000312.

Es decir, para un c entre 0 y \frac{\pi}{4}, los términos k = 0, 1, \dots, 8 de la serie serán los únicos que podrían llegar a ser mayores que 1 \times 10^{-15}. En otras palabras, los términos k = 0, 1, \dots, 8 serán más que suficientes para cualquier c entre 0 y \frac{\pi}{4}, ya que, para todo c entre 0 y \frac{\pi}{4}, los términos k = 9 en adelante serán menores que 1 \times 10^{-15}.