From 6f6c22587d809972b1f9c2ba65f6319b357bb1bd Mon Sep 17 00:00:00 2001
From: Sasiri Juliana Vargas Urbano <vargass@jupyterMiLAB>
Date: Tue, 18 May 2021 05:32:35 -0500
Subject: [PATCH] primer commit
---
main.py | 427 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 427 insertions(+)
create mode 100644 main.py
diff --git a/main.py b/main.py
new file mode 100644
index 0000000..7d9050b
--- /dev/null
+++ b/main.py
@@ -0,0 +1,427 @@
+import random as rm
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+
+from matplotlib.animation import FuncAnimation #Para las animaciones
+import moviepy.editor as mp #Para transformar el gif en mp4
+from IPython.display import Video #Para reproducir el video mp4 notebook
+
+#------------
+# DINÃMICA
+#------------
+def Dinamica(vec, l, Pr, Pj, Ps):
+
+ #Vector que almacena las concentraciones de las partÃculas A y B
+ concentracion=np.array(vec)
+ #Vector que almacena la distribución de espaciamientos
+ Pesp=np.zeros(len(concentracion))
+
+ # Números aleatorios:
+ #--------------------
+
+ #Se tienen 4 números aleatorios. 3 de ellos representan la posibilidad de tener o no,
+ #reacción, incerción y salto. El restante es el sitio en el que ocurrirá la dinámica
+
+ Escogida = rm.randint(0,2*l) #Num aleatorio que escoge el sitio de red para efectuar la dinpamica
+ particula = concentracion[Escogida] #Tipo de partÃcula escogida
+
+ p_r = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de reacción
+ p_j = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de incerción
+ p_s = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de salto
+
+
+ # Saltos
+ # -------
+
+ #Para los saltos se tiene en cuenta lo siguiente:
+ #Si 0 < p_s < Ps -> ocurre un salto hacia la derecha
+ #Si Ps < p_s < 2*Ps -> ocurre un salto hacia la izquierda
+ #(Siempre y cuando el sitio al que saltará la partÃcula este vacÃo)
+ #Esto se realiza para que podamos tener información tanto de la posibilidad de salto,
+ #como de la posibilidad de no saltar
+
+ #Ejemplo: Si Ps=0.5, entonces:
+ #Salto a la derecha : 0 < p_s < 0.5 --> 50% de las veces
+ #Salto a la izquierda: 0.5 < p_s < 1 --> 50% de las veces
+ #No salta --> 0% de las veces
+
+ if p_s < 2*Ps:
+
+ #Salto general (Se pueden escoger todos los sitios, excepto los extremos)
+ if all((Escogida != 0, Escogida != 2*l, concentracion[Escogida] != 0)):
+ #Salto derecha
+ if all((p_s < Ps, concentracion[Escogida + 1] == 0)):
+ concentracion[Escogida + 1] = particula
+ concentracion[Escogida] = 0
+ #Salto izquierda
+ elif all((Ps < p_s < 2*Ps, concentracion[Escogida - 1] == 0)):
+ concentracion[Escogida - 1] = particula
+ concentracion[Escogida] = 0
+
+ #Extremo izquierdo (dominio de la partÃcula A)
+ elif Escogida == 0:
+ #Incercion: Si p_j < Pj, entonces se efectúa la incerción
+ if all((p_j < Pj, concentracion[Escogida] == 0)):
+ concentracion[Escogida] = -1
+ #Salto: Se considera que hay salto a la derecha
+ elif all((p_s < Ps, concentracion[0] != 0, concentracion[1] == 0)):
+ concentracion[1] = particula
+ concentracion[0] = 0
+
+ #Extremo derecho (dominio de la partÃcula B)
+ elif Escogida == 2*l:
+ #Incersion: Si p_j < Pj, entonces se efectúa la incerción
+ if all((p_j < Pj, concentracion[Escogida] == 0)):
+ concentracion[Escogida] = 1
+ #Salto: Se considera que hay salto a la izquierda
+ elif all((p_s < Ps, concentracion[2*l] != 0, concentracion[2*l-1] == 0)):
+ concentracion[2*l-1] = particula
+ concentracion[2*l] = 0
+
+ # Reaccion
+ #----------
+
+ #La reacción se efectúa si p_r < Pr y si la distancia entre rma y lmb es igual a 1 (si A y B son vecinos)
+
+ try:
+ #El try y el except es puesto porque hay ocasiones en las que se tiene por ejemplo:
+ #[-1,-1,-1,0,0,0]
+ #Entonces en este caso no es posible encontrar un lmb
+
+ #rma
+ where1=np.where(concentracion == -1) #devuelve una tupla
+ rmav=where1[0] #Arreglo de las ubicaciones de las partÃculas A
+ rma=rmav[-1] #PartÃcula más a la derecha
+
+ #lmb
+ where2=np.where(concentracion == 1)
+ lmbv=where2[0] #Arreglo de las ubicaciones de las partÃculas B
+ lmb=lmbv[0] #PartÃcula más a la izquierda
+
+ if (p_r < Pr) and (abs(rma-lmb) == 1): #Las partÃculas se aniquilan
+ concentracion[rma] = 0
+ concentracion[lmb] = 0
+
+ #Se encuentran los nuevos rma y lmb para encontrar la distr. de espaciamientos
+ #rma
+ where1=np.where(concentracion == -1)
+ rmav=where1[0]
+ rma=rmav[-1]
+
+ #lmb
+ where2=np.where(concentracion == 1)
+ lmbv=where2[0]
+ lmb=lmbv[0]
+
+ #Función que devuelve un arreglo de la distribución de espaciamientos
+ Pesp= P_i_xt(rma,lmb,l)
+
+ except:
+ pass
+
+ return concentracion, Pesp # Devuelve 2 arreglos
+
+# ---------
+# CICLOS
+# ---------
+
+#Repite la función 'Dinamica' una cantidad M de veces, determinada por 'time'
+
+def ciclos(c_ini, l, Pr, Pj, Ps, time):
+
+ resultado=np.copy(c_ini) # copia de la condición inicial
+ salida_conc_A=np.zeros(len(resultado)) # arreglo de salida para la concentración de A
+ salida_conc_B=np.zeros(len(resultado)) # arreglo de salida para la concentración de B
+ salida_esp=np.zeros(len(resultado)) # arreglo de salida para la dist. de espaciamientos
+
+ for i in range(time+1):
+ resultado, espaciamiento = Dinamica(resultado, l, Pr, Pj, Ps) #Dinamica para el tiempo i
+ A, B = concentracionA_B(resultado) # Separación de arreglos
+
+ return A, B, espaciamiento # Devuelve 3 arreglos
+
+
+#------------------------
+# FUNCIONES ADICIONALES
+#------------------------
+
+#------------------------------------------
+# Distribución de espaciamientos discreta
+#------------------------------------------
+
+#La distribución de espaciamiento discreta se encuentra con la función P_i_xt(r_,l_,L_),
+#la cual utiliza P_i(r,l,i).
+
+#Devuelve un arreglo por ejemplo: [0,0,1,0,0,0]
+#Indicando que la distancia entre rma y lmb es 2
+
+def P_i(r,l,i):
+ salida=0
+ if (abs(r-l)==i):
+ salida=1
+ return salida
+
+def P_i_xt(r_,l_,L_):
+ return [P_i(r_,l_,j) for j in range(2*L_+1)]
+
+
+#------------------------------------
+# Separación de concentraciones
+#------------------------------------
+
+#Entra un arreglo del tipo: [-1,-1, 0, 1, 1]
+# Y devuelve dos arreglos: [ 1, 1, 0, 0, 0] --> Para las partÃculas A
+# [ 0, 0, 0, 1, 1] --> Para las partÃculas B
+
+def concentracionA_B(entrada):
+ salida_A = np.zeros(len(entrada))
+ salida_B = np.zeros(len(entrada))
+
+ for i in range(len(entrada)):
+
+ if entrada[i] == -1:
+ salida_A[i] = 1
+
+ elif entrada[i] == 1:
+ salida_B[i] = 1
+
+ return salida_A, salida_B
+
+#--------------
+# ANALISIS
+#--------------
+def ANALISIS(L=5, dx=1, dt=1, D=0.5, j=1, K=1, tiempo=2000, rep=10, opcion='graficos', nombre = 'Data'):
+
+ #PARÃMETROS
+
+ beta= K/dx # Tasa de reacción
+ lamb = D/dx**2 #Tasa de salto
+ P_S = lamb * dt # Probabilidad de salto
+ P_J= j #Probabilidad de incersion
+ P_R= beta*dt #Probabilidad de reacción
+
+ #CONDICION INICIAL
+
+ # 0 : VacÃo
+ # -1 : PartÃcula A
+ # 1 : Particula B
+
+ condicion_inicial=[-1 if i<L else 1 for i in range(2*L+1)] #Ilustracion para L=3: [-1,-1,-1,1,1,1]
+
+#-----------------------------------OPCIONES------------------------------------
+
+ if opcion=='graficos':
+
+ x = np.linspace(-L, L, 2*L+1) #Posiciones
+ z = np.linspace(0, 2*L, 2*L+1) #Distancias dist. de espaciamientos
+
+ salida_A=np.zeros(len(x))
+ salida_B=np.zeros(len(x))
+ salida_Desp=np.zeros(len(z))
+
+ # Se almacena la información de TODA la dinámica en 3 arreglos de salida
+ # Se suman cada iteración de tiempo para posteriormente hacer su respectiva normalización y promedio
+ # y encontrar la densidad local contÃnua y dist. de espaciamientos contÃnua
+
+ for i in range(0,rep):
+ A, B, d_esp = ciclos(condicion_inicial, L, P_R, P_J, P_S, tiempo)
+
+ salida_A = salida_A + A
+ salida_B = salida_B + B
+ salida_Desp = salida_Desp + d_esp
+
+ # Cada 1000 iteraciones sacamos un print en pantalla para conocer el estado de la simulacion
+ if not(i % 1000):
+ print('Working in iteration number {} out of {}.'.format(i, rep))
+
+ salida_A = salida_A/(rep*2*L)
+ salida_B = salida_B/(rep*2*L)
+ salida_Desp = salida_Desp/salida_Desp.sum()
+
+
+ #--------------------------------------------------------
+ #------------------ARCHIVOS------------------------------
+ #--------------------------------------------------------
+
+ #Escribe 3 archivos txt en la carpeta Archi:
+ # - Info para A
+ # - Info para B
+ # - Info para dist. de espaciamientos
+
+ print(nombre)
+
+ nombre1 = 'Archi/'+str(nombre)+'A.txt'
+ nombre2 = 'Archi/'+str(nombre)+'B.txt'
+ nombre3 = 'Archi/'+str(nombre)+'Desp.txt'
+
+ File_A = open(nombre1, 'w')
+ File_B = open(nombre2, 'w')
+ File_P_i = open(nombre3, 'w')
+
+ for i in range(len(x)):
+ File_A.write(str(x[i])+' '+str(salida_A[i])+'\n')
+ File_B.write(str(x[i])+' '+str(salida_B[i])+'\n')
+ File_P_i.write(str(z[i])+' '+str(salida_Desp[i])+'\n')
+
+ File_A.close()
+ File_B.close()
+ File_P_i.close()
+
+
+ #----------------------------------------------------------
+ #----------------------GRAFICOS----------------------------
+ #----------------------------------------------------------
+
+ #-------------------PERFILES DE DENSIDAD-------------------------
+ fig1=plt.figure(1)
+ plt.plot(x,salida_A,c='b')
+ plt.plot(x,salida_B,c='r')
+ plt.legend(['Particulas A', 'Particulas B'],loc='upper right',bbox_to_anchor=(1.35, 1))
+ plt.xlabel('x')
+ plt.ylabel('Perfil de densidad')
+ plt.show()
+
+ #---------------DISTRIBUCIÓN DE ESPACIAMIENTOS--------------------
+ fig2=plt.figure(2)
+ plt.plot(z,salida_Desp,c='b', alpha=0.5)
+ plt.xlabel('x')
+ plt.ylabel('Distribución de espaciamientos')
+ plt.show()
+
+ #-------------------------------------------------------------------
+
+ #---------------
+ # Dinamica
+ #---------------
+
+ #Imprime en cada paso de tiempo, el vector concentración
+ #Con esta opción se puede visualizar lo que ocurre en cada paso de tiempo
+
+ elif opcion=='dinamica':
+
+ salida_din=np.copy(condicion_inicial)
+
+ for i in range(0,tiempo):
+ print('tiempo ', i)
+ salida_din, otro = Dinamica(salida_din, L, P_R, P_J, P_S)
+ print(salida_din)
+ print(' ')
+
+ #--------------------------------------------------------------------
+
+ #---------------
+ # Estacionario
+ #---------------
+
+ #Comprueba si se llega al estado estacionario en el tiempo ingresado. Esto se hace comparando
+ #el sistema en el tiempo = 'tiempo' y tiempo = 'tiempo' + 1000
+
+ elif opcion=='estacionario':
+
+ pruebas=[]
+
+ #Evolución hasta el tiempo = 'tiempo'
+ A_est, B_est, Desp_est = ciclos(condicion_inicial, L, P_R, P_J, P_S, tiempo)
+
+ estacionario = B_est - A_est
+
+ pruebas.append(estacionario)
+ print('Estado en la iteración '+ str(tiempo)+':')
+ print(estacionario)
+ print(' ')
+
+ #Evolución hasta el tiempo = 'tiempo + 1000'
+ A_est2, B_est2, Desp_est2 = ciclos(estacionario, L, P_R, P_J, P_S, 1000)
+
+ estacionario2 = B_est2 - A_est2
+
+ pruebas.append(estacionario2)
+ print('Estado en la iteración '+ str(tiempo+1000)+':')
+ print(estacionario2)
+ print(' ')
+
+ salida_estacionario = pruebas[0] == pruebas[1]
+
+ num_verdadero = list(salida_estacionario).count(True)
+
+ print('Estado estacionario: ' + str(len(salida_estacionario) == num_verdadero))
+
+ return len(salida_estacionario) == num_verdadero
+
+
+#-----------CREACION DEL CANVAS--------------
+
+L = 5 # 1/2 de la longitud del sistema
+x = np.linspace(-L, L, 2*L+1) #Representacion del arreglo
+y = np.zeros(2*L+1) #Para representar el arreglo unidimensional
+
+fig = plt.figure()
+ax = plt.axes(xlim=(-L-1, L+1), ylim=(-0.005,0.005))
+scatt = ax.scatter([],[],s=70)
+plt.close()
+
+
+#------------FUNCIONES PARA FUNCANIMATION--------------------
+
+# def init():
+# scatt.set_offsets([])
+# return scatt
+
+cmap = matplotlib.cm.bwr
+norm = matplotlib.colors.Normalize(vmin=-1,vmax=1)
+
+def animate(i, *fargs):
+
+ concentracion = fargs[i]
+
+ scatt.set_offsets(np.c_[x,y]) #Para correr los datos
+ scatt.set_color(cmap(norm(concentracion))) #Para los colores de la gráfica
+ scatt.set_edgecolor('black')
+
+ return scatt,
+
+
+#-----------FUNCIÓN PARA LAS ANIMACIONES (CONTIENE FUNCANIMATION)----------
+
+def animacion(time, Pr, Pj, Ps):
+
+ #Longitud del arreglo predefinida:
+ l=5
+
+ # Se genera la dinámica del sistema utilizando la función 'Dinamica'
+ # hasta llegar al tiempo 'time', los resultados obtenidos se almacenan
+ # en la lista vec_animate
+
+ vec_animate = []
+
+ condicion_inicial=[-1 if i<l else 1 for i in range(2*l+1)]
+
+ salida = np.copy(condicion_inicial)
+
+ for i in range(time):
+ salida, otro = Dinamica(salida, l, Pr, Pj, Ps)
+ vec_animate.append(salida)
+
+
+ # Se hace uso de FuncAnimation
+
+ anim = FuncAnimation(fig, animate, frames=time - 1, interval=200, blit=True, repeat=False, fargs = vec_animate)
+
+
+ # Se guarda como gif
+
+ nombre_gif = 'anim.gif'
+
+ anim.save(nombre_gif, writer=matplotlib.animation.PillowWriter())
+
+
+ # Se transforma el gif en mp4
+
+ nombre_video = 'video_dinamica.mp4'
+
+ clip = mp.VideoFileClip(nombre_gif)
+ clip.write_videofile(nombre_video)
+
+ return nombre_video
+
--
GitLab