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   958  1127.5  1311.191  1237.0  1373.0   14896  1000
     mi_fact_klen1_POA   746  1002.0  1187.126  1147.0  1291.5    6041  1000
          mi_fact_PIEP  1891  2594.0  3144.303  3128.0  3406.0   20533  1000
 mi_fact_POA_Vectorize 30691 33200.0 35104.018 33988.5 35106.0  363835  1000
    mi_fact_POA_sapply 11398 12944.0 15274.747 13668.0 14341.5 1285978  1000
             factorial   372   499.0   598.559   568.5   644.0    9912  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 34.162 39.3040 45.107212 41.2065 42.9405 2767.120  1000
 mi_fact_POA_Vectorize 66.424 71.5705 77.771241 74.5260 77.8560  171.114  1000
    mi_fact_POA_sapply 43.561 48.1915 54.999933 50.1345 52.4745 3083.187  1000
             factorial  1.574  1.9230  2.218272  2.1460  2.3370   12.964  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}.