En este taller vamos a trabajar con una implementación en Python del modelo de difusión espacial de Torsten Hägerstrand. Este modelo es una simulación tipo Monte Carlo del proceso de difusión de innovaciones concebido originalmente por Hagerstrand en 1953. Para este taller vamos a trabajar con la versión más simple del modelo, las suposiciones básicas son las siguientes:
En escencia, lo que hacemos es considerar al espacio geográfico como una retícula en la que cada elemento tiene el mismo número de habitantes. A partir de aquí, la probabilidad de que un habitante de una retícula contacte a otro habitante es función únicamente de la distancia.
La simulación Monte Carlo ocurre en dos etapas, para cada portador del mensaje:
Vamos ahora a jugar un poco con la implementación del modelo de Hagerstrand que está en esta misma carpeta.
In [1]:
%pylab inline
from haggerstrand.diffusion import SimpleDiffusion
En la primera linea importamos todas las librerías que forman el paquete pylab y le decimos que queremos graficar inline, es decir, aquí mismo. La segunda linea importa la clase SimpleDiffusion que es en donde está implementado el algoritmo.
Para correr el modelo, primero tenemos que instanciar la clase:
In [ ]:
s = SimpleDiffusion(50,50,9,20,[(20,20)],0.3,15)
Los parámetros para La clase SimpleDiffusion son los siguientes (en orden):
Para ver la lista de parámetros, sus tipos y una descripción, pueden consultar la ayuda:
In [ ]:
help(SimpleDiffusion)
Ahora que ya tenemos instanciada la clase, podemos correr la simulación:
In [ ]:
s.spatial_diffusion()
El algoritmo guarda los resultados para cada iteración en una matriz de M x N x max_iter, el resultado final lo podemos ver en la última rebanada de la matriz:
In [ ]:
s.result[:,:,s.max_iter-1]
Fíjense en la notación: s.result[:,:,s.max_iter-1] le está pidiendo al array s.result que le dé los dos primeros índices completos y la entrada s.max_iter-1 del último índice. (Noten que el último índice está en max_iter-1 porque los arrays empiezan en cero)
Habrán notado que usamos la notación s.algo para hablar con el objeto s, esta es la manera de hablar con las propiedades de los objetos en Python. La lista de propiedades de SimpleDiffusion las encuentran en la sección attributes de la ayuda.
Claro que ver los números de la matriz de resultado no es la mejor manera de visualizarlos, hagamos una gráfica:
In [ ]:
plt.imshow(s.result[:,:,s.max_iter-1])
In [ ]:
plt.imshow(s.result[:,:,3])
Para entender el proceso de difusión espacial, conviene compararlo con un proceso que suceda aleatoriamente, para eso SimpleDiffusion tiene el método random_diffusion:
In [ ]:
s.random_diffusion()
Esto reescribió la matriz de resultados y la populó con una difusión aleatoria. Veamos el resultado final:
In [ ]:
plt.imshow(s.result[:,:,14])
Como reescribimos la matriz de resultados, necesitamos guardar el resultado y volver a correr el algoritmo para comparar ambos casos:
In [ ]:
random_diff = np.copy(s.result)
s.spatial_diffusion()
Ahora sí podemos ver el resultado final en ambos casos:
In [ ]:
subplot(1,2,1)
plt.imshow(s.result[:,:,s.max_iter-1])
subplot(1,2,2)
plt.imshow(random_diff[:,:,s.max_iter-1])
Como vimos, los resultados son muy diferentes para la disusión espacial y la difusión aleatoria, vamos a ver ahora cómo se comporta la evolución temporal. Para esto, la clase SimpleDiffussion nos provee la varable time_series que es un lista con la cantidad de nuevos adoptantes en cada iteración:
In [ ]:
s.time_series
Claro que la podemos graficar:
In [ ]:
plt.plot(s.time_series)
Como podemos ver, el crecimiento parece ser exponencial (piensen por qué), pero ¿Qué pasa si aumentamos el número de iteraciones?
In [ ]:
s = SimpleDiffusion(50,50,9,20,[(20,20)],0.3,20)
s.spatial_diffusion()
plt.plot(s.time_series)
In [2]:
from JSAnimation import IPython_display
import matplotlib.animation as animation
s = SimpleDiffusion(50,50,9,20,[(20,20)],0.3,18)
s.spatial_diffusion()
ims = []
fig = plt.figure()
for i in range(0,s.max_iter):
im = imshow(s.result[:,:,i])
ims.append([im])
animation.ArtistAnimation(fig, ims, interval=100, blit=True)
Out[2]:
El código para crear la animación parece complicado pero, si se fijan con atención, en realidad es bastante directo: creamos una lista vacía para guardar cada imagen (ims = []) y luego la vamos llenando adentro de un ciclo de control, al final llamamos a la librería matplotlib.animation para que se encargue de hacer la animación con las imágenes que ya construimos (la linea from JSAnimation import IPython_display sólo sirve para desplegar la animación en ésta interfaz, normalmente no es necesaria).
Hagan una animación para el caso de difusión aleatoria
In [ ]: