El Mapa Logístico: un ejercicio de programación y graficación científica.
Hace unas semanas y por pura casualidad (¿o habrá sido Caos?) me topé con el video de Youtube “This equation will change how you see the world” y me sentí de inmediato inspirado, no solo por la temática y la calidad del video en cuestión, sino porque el material me transportó de vuelta a mis clases de pregrado y a los cursos de física computacional, cuando me topé por primera vez y de manera formal con el concepto del Caos y los Sistemas Dinámicos.
Más aún, el revisitar estos conceptos me hizo recordar que ya en aquel entonces me tocó, como parte de algún proyecto o tarea, graficar el mismo Diagrama de Bifurcación del Mapa Logístico. En ese entonces lo hice en lenguaje C, si mal no recuerdo usando las primitivas gráficas de X11 en Linux. Por supuesto, pasando mucho trabajo y dolores de cabeza. Pero justo tras terminar de ver el excelente video, pensé: ¡apuesto a que con esto se podría hacer un bonito gráfico en R!
Así que, minutos después, estaba ya frente al computador tronando los dedos y preparándome para llevar a cabo un satisfactorio ejercicio de programación, ciencia de datos y graficación científica.
El primer paso, por supuesto, era obtener los datos necesarios para trazar la curva generada por el Mapa Logístico, definido como:
Pero curiosamente mi primera reacción fue levantar Spyder y empecé a escribir el código en Python. Cuando se trata de escribir funciones, generar datos y hacer prototipado rápido de algoritmos, considero que Python suele ser más rápido y sencillo de usar. Así que, por supuesto, lo primero fue importar las librerías necesarias:
import numpy as np
import matplotlib.pyplot as plt
from random import randint
import pandas as pdY luego, implementar una función que permitiera iterar el Mapa Logístico un número N de veces, a partir de un valor de X inicial, así como de un valor de r dado.
Así, tenemos:
def logistic_iteration(N, xinicial, r):
logistic=0
for i in range(N):
logistic = r*xinicial*(1-xinicial)
xinicial = logistic
return logisticEl paso siguiente es simplemente definir una lista vacía para almacenar los valores que tomará la función tras cada grupo de iteraciones, otra lista para almacenar los valores de r seleccionados a lo largo del eje de las abscisas, y definir en un intervalo regular, en este caso entre 0 y 4, las cantidades distintas de valores de dicho parámetro r que se emplearán para trazar la curva:
valor_f = []
eje_r = []
valor_r = list(np.linspace(0, 4, 100000))Luego, solo basta con iterar a lo largo de los elementos de esta última lista, y aplicar la función de iteración seleccionando al azar en un rango dado el valor de N en cada caso, así como el valor inicial de X:
N = 2000
Xo = 0.5for elems in valor_r:
print(elems)
log = 0
log = logistic_iteration(randint(100, N), Xo, elems)
valor_f.append(log)
eje_r.append(elems)
Finalizadas la iteraciones, podemos observar el resultado haciendo un gráfico sencillo:
plt.scatter(eje_r, valor_f)
plt.show()
Como podemos ver, se obtiene exactamente lo que buscamos: el Diagrama de Bifurcación, con sus distintas regiones periódicas, y luego la región caótica.
Ya que estamos satisfechos con la cantidad de puntos y forma de los datos, podemos almacenar esto en la forma de un dataframe, y guardarlo luego como un archivo separado por comas:
datos = {‘r’:eje_r, ‘xn’:valor_f}
data = pd.DataFrame(datos, dtype = float)
data.to_csv(‘dataLogisticMap.csv’, index = False)Y así tenemos listo el dataset que emplearemos en R para construir nuestra gráfica.
Una vez en RStudio, primero que nada se deben cargar las librerías a emplear:
library(tidyverse)
library(viridis)
library(latex2exp)Luego, cargamos el dataset recién creado:
dataset <- read_csv('dataLogisticMap.csv')Y al echarle un ojo a los datos:
dataset## # A tibble: 100,000 x 2
## r xn
## <dbl> <dbl>
## 1 0 0
## 2 0.0000400 0
## 3 0.0000800 0
## 4 0.000120 0
## 5 0.000160 0
## 6 0.000200 0
## 7 0.000240 0
## 8 0.000280 0
## 9 0.000320 0
## 10 0.000360 0
## # ... with 99,990 more rows
Vemos que tenemos un dataframe de 100.000 filas, y solo 2 columnas: el valor del parámetro r de la ecuación, y el valor del mapa logístico iterado en cada caso.
Solo nos basta esto para generar un primer gráfico:
dataset %>% ggplot(aes(x=r, y=xn)) + geom_point()
¡Muy bien! Ya tenemos nuestro Diagrama de Bifurcación en R. Aunque tal vez no con la estética que deseamos. Bien, llegó el momento de empezar a afinar la visualización.
En primer lugar, vamos a asociar el color de los puntos al valor del mapa en cada iteración, es decir, agregemos el parámetro colour a la estética de la función ggplot():
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point()
Mucho mejor, pero no me gusta la leyenda, así que vamos a eliminarla aprovechando los parámetros de la función theme() y que, además, se empleará luego para afinar todavía más parámetros:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point() +
theme(legend.position = 'None')
Perfecto. Ya tenemos el Diagrama de Bifurcación en toda su gloria, pero los colores por defecto no me convencen, así que (como se observó cuando cargué las librerías y como lo tenía pensado incluso antes de escribir la primera línea de código) vamos a emplear la escala de color de la librería viridis:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point() +
scale_color_viridis() +
theme(legend.position = 'None')
¡Ya la gráfica empieza a tomar forma! Lo siguiente es cambiar el color del fondo. Desde el principio sabía que quería que el fondo fuese negro, pues resulta sin duda unos de los mejores colores para hacer contraste a los tonos que ofrece la escala viridis. Así que, hagamos negro el fondo:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point() +
scale_color_viridis() +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000"))
¡Oops! Hace falta también colocar de negro el color del “panel” interior del gráfico, y para eliminar el contorno blanco del fondo, hagamos ‘NA’ el color de la línea:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point() +
scale_color_viridis() +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000", color = NA),
panel.background = element_rect(fill = "#000000"))
Hmmm, vamos poco a poco. Ahora, considero necesario incorporar el título y subtítulo del gráfico a fin de empezar a ver cómo se distribuyen todos los elementos. Por supuesto, es necesario cambiar el color de dicho texto pues de lo contrario no se observará en la gráfica porque, por defecto, es negro:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point() +
scale_color_viridis() +
ggtitle(TeX('Diagrama de Bifurcación del Mapa Logístico $x_{n+1} = rx_{n}(1-x_{n})$')) +
labs(subtitle = 'Borde del Caos r ~ 3.57') + xlab(TeX('$r$')) + ylab(TeX('$x_{n+1}$')) +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000", color = NA),
panel.background = element_rect(fill = "#000000"),
title = element_text(color = "#abc1d3", face = "bold"))
Como se puede observar, en el título y los ejes hago uso de la función TeX() de la librería latex2exp a fin de poder incorporar ecuaciones al texto.
Ahora que tengo esto, observo que la cuadrícula interior del gráfico en color blanco no es de mi agrado, así que el paso siguiente es modificarla tanto de color, como de tipo de línea y tamaño:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point() +
scale_color_viridis() +
ggtitle(TeX('Diagrama de Bifurcación del Mapa Logístico $x_{n+1} = rx_{n}(1-x_{n})$')) +
labs(subtitle = 'Borde del Caos r ~ 3.57') + xlab(TeX('$r$')) + ylab(TeX('$x_{n+1}$')) +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000", color = NA),
panel.background = element_rect(fill = "#000000"),
title = element_text(color = "#abc1d3", face = "bold"),
panel.grid.major.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.major.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"))
Mucho mejor. Ahora, quiero modificar el texto y color tanto de la numeración de los ejes como de los títulos de los ejes:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point() +
scale_color_viridis() +
ggtitle(TeX('Diagrama de Bifurcación del Mapa Logístico $x_{n+1} = rx_{n}(1-x_{n})$')) +
labs(subtitle = 'Borde del Caos r ~ 3.57') + xlab(TeX('$r$')) + ylab(TeX('$x_{n+1}$')) +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000", color = NA),
panel.background = element_rect(fill = "#000000"),
title = element_text(color = "#abc1d3", face = "bold"),
panel.grid.major.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.major.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
axis.title.y = element_text(color = "#ffffff", face = "bold", family = "mono"),
axis.title.x = element_text(color = "#ffffff", face = "bold", family = "mono"),
axis.text.x = element_text(color = "#ffffff", family = "mono"),
axis.text.y = element_text(color = "#ffffff", family = "mono"),
text=element_text(family = "TT Arial"))
También, vamos a agregar un ‘caption’ al gráfico, con la información del autor, así como unos pequeños márgenes:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point() +
scale_color_viridis() +
ggtitle(TeX('Diagrama de Bifurcación del Mapa Logístico $x_{n+1} = rx_{n}(1-x_{n})$')) +
labs(subtitle = 'Borde del Caos r ~ 3.57',
caption = 'Ronald Delgado | linkedin.com/in/ronald-delgado') +
xlab(TeX('$r$')) +
ylab(TeX('$x_{n+1}$')) +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000", color = NA),
panel.background = element_rect(fill = "#000000"),
title = element_text(color = "#abc1d3", face = "bold"),
panel.grid.major.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.major.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
axis.title.y = element_text(color = "#ffffff", face = "bold", family = "mono"),
axis.title.x = element_text(color = "#ffffff", face = "bold", family = "mono"),
axis.text.x = element_text(color = "#ffffff", family = "mono"),
axis.text.y = element_text(color = "#ffffff", family = "mono"),
text=element_text(family = "TT Arial"),
plot.margin=unit(c(1.0,1.0,1.0,1.0),"cm"))
Ya nuestro gráfico está casi listo. Sin embargo, veo que los puntos del diagrama están muy grandes, por lo que vamos a disminuir su tamaño. También, me interesa colocar una linea vertical en el umbral en donde comienza la región caótica de la curva (r ~ 3.57):
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point(size = 0.1) +
geom_vline(xintercept = 3.57, color = 'red') +
scale_color_viridis() +
ggtitle(TeX('Diagrama de Bifurcación del Mapa Logístico $x_{n+1} = rx_{n}(1-x_{n})$')) +
labs(subtitle = 'Borde del Caos r ~ 3.57',
caption = 'Ronald Delgado | linkedin.com/in/ronald-delgado') +
xlab(TeX('$r$')) +
ylab(TeX('$x_{n+1}$')) +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000", color = NA),
panel.background = element_rect(fill = "#000000"),
title = element_text(color = "#abc1d3", face = "bold"),
panel.grid.major.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.major.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
axis.title.y = element_text(color = "#ffffff", face = "bold", family = "mono"),
axis.title.x = element_text(color = "#ffffff", face = "bold", family = "mono"),
axis.text.x = element_text(color = "#ffffff", family = "mono"),
axis.text.y = element_text(color = "#ffffff", family = "mono"),
text=element_text(family = "TT Arial"),
plot.margin=unit(c(1.0,1.0,1.0,1.0),"cm"))
Agregemos, además, una anotación en la región caótica, sólo por pura estética.
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point(size = 0.1) +
geom_vline(xintercept = 3.57, color = 'red') +
geom_text(x = 3.57, y = 0.05, label = "Región Caótica", color = "red", hjust = -0.2) +
scale_color_viridis() +
ggtitle(TeX('Diagrama de Bifurcación del Mapa Logístico $x_{n+1} = rx_{n}(1-x_{n})$')) +
labs(subtitle = 'Borde del Caos r ~ 3.57',
caption = 'Ronald Delgado | linkedin.com/in/ronald-delgado') +
xlab(TeX('$r$')) +
ylab(TeX('$x_{n+1}$')) +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000", color = NA),
panel.background = element_rect(fill = "#000000"),
title = element_text(color = "#abc1d3", face = "bold"),
panel.grid.major.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.major.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
axis.title.y = element_text(color = "#ffffff", face = "bold", family = "mono"),
axis.title.x = element_text(color = "#ffffff", face = "bold", family = "mono"),
axis.text.x = element_text(color = "#ffffff", family = "mono"),
axis.text.y = element_text(color = "#ffffff", family = "mono"),
text=element_text(family = "TT Arial"),
plot.margin=unit(c(1.0,1.0,1.0,1.0),"cm"))En este punto, deberíamos tener a la salida un gráfico casi definitivo que luce como el siguiente:
Sin embargo, todavia podemos jugar un poco más con el tamaño del texto y los demás elementos, hasta obtener lo siguiente:
dataset %>% ggplot(aes(x=r, y=xn, colour=xn)) + geom_point(size = 0.1) +
geom_vline(xintercept = 3.57, color = 'red') +
geom_text(x = 3.57, y = 0.05, label = "Región Caótica", color = "red", hjust = -0.2, size = 4) +
scale_color_viridis() +
ggtitle(TeX('Diagrama de Bifurcación del Mapa Logístico $x_{n+1} = rx_{n}(1-x_{n})$')) +
labs(subtitle = 'Borde del Caos r ~ 3.57',
caption = 'Ronald Delgado | linkedin.com/in/ronald-delgado') +
xlab(TeX('$r$')) +
ylab(TeX('$x_{n+1}$')) +
theme(legend.position = 'None',
plot.background = element_rect(fill = "#000000", color = NA),
panel.background = element_rect(fill = "#000000"),
title = element_text(color = "#abc1d3", face = "bold", size = 22),
panel.grid.major.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.x = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.major.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
panel.grid.minor.y = element_line(size = 0.3,
linetype = 'dashed',
colour = "#3f4d58"),
axis.title.y = element_text(color = "#ffffff", face = "bold", family = "mono", size = 20),
axis.title.x = element_text(color = "#ffffff", face = "bold", family = "mono", size = 20),
axis.text.x = element_text(color = "#ffffff", family = "mono", size = 12),
axis.text.y = element_text(color = "#ffffff", family = "mono", size = 12),
text=element_text(family = "TT Arial"),
plot.margin=unit(c(1.0,1.0,1.0,1.0),"cm"),
plot.caption = element_text(size = 11),
plot.subtitle = element_text(size = 16))¡Hemos terminado nuestro gráfico!
Queda a gusto del usuario o autor afinar los tamaños, colores y demás parámetros del gráfico, y luego guardarlo con las dimensiones adecuadas.
Como se vió, el ejercicio de programar el Mapa Logístico, generar sus datos y luego transformarlos en una bonita visualización del Diagrama de Bifurcación no solo sirvió para obtener una gráfica de gran calidad, sino como práctica tanto de Python, R y ciencia de datos en general.
Veamos… ¿Qué otras gráficas basadas en la Teoría del Caos podemos elaborar?
