From eff682665500cbbf0a971b056afb33bc614f23f6 Mon Sep 17 00:00:00 2001
From: Sasiri Juliana Vargas Urbano <sasiri.vargas@correounivalle.edu.co>
Date: Tue, 25 May 2021 18:42:40 +0000
Subject: [PATCH] Upload New File
---
Numba.ipynb | 483 ++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 483 insertions(+)
create mode 100644 Numba.ipynb
diff --git a/Numba.ipynb b/Numba.ipynb
new file mode 100644
index 0000000..1c26fa2
--- /dev/null
+++ b/Numba.ipynb
@@ -0,0 +1,483 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import random as rm\n",
+ "import numpy as np\n",
+ "import matplotlib\n",
+ "import matplotlib.pyplot as plt\n",
+ "import time\n",
+ "from numba import jit, njit, prange\n",
+ "import multiprocessing as mp"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#------------------------\n",
+ "# FUNCIONES ADICIONALES\n",
+ "#------------------------\n",
+ "\n",
+ "#------------------------------------------\n",
+ "# Distribución de espaciamientos discreta\n",
+ "#------------------------------------------\n",
+ "\n",
+ "#La distribución de espaciamiento discreta se encuentra con la función P_i_xt(r_,l_,L_),\n",
+ "#la cual utiliza P_i(r,l,i).\n",
+ "\n",
+ "#Devuelve un arreglo por ejemplo: [0,0,1,0,0,0]\n",
+ "#Indicando que la distancia entre rma y lmb es 2\n",
+ "\n",
+ "@njit\n",
+ "def P_i(r,l,i):\n",
+ " salida=0.\n",
+ " if (abs(r-l)==i):\n",
+ " salida=1.\n",
+ " return salida\n",
+ "\n",
+ "@njit\n",
+ "def P_i_xt(r_,l_,L_): \n",
+ " salida = np.zeros(2*L_+1) \n",
+ " for j in range(0,2*L_+1):\n",
+ " salida[j] = P_i(r_,l_,j)\n",
+ " return salida"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#------------------------------------\n",
+ "# Separación de concentraciones\n",
+ "#------------------------------------\n",
+ "\n",
+ "#Entra un arreglo del tipo: [-1,-1, 0, 1, 1]\n",
+ "# Y devuelve dos arreglos: [ 1, 1, 0, 0, 0] --> Para las partÃculas A \n",
+ "# [ 0, 0, 0, 1, 1] --> Para las partÃculas B\n",
+ "@njit\n",
+ "def concentracionA_B(entrada):\n",
+ " salida_A = np.zeros(len(entrada))\n",
+ " salida_B = np.zeros(len(entrada))\n",
+ "\n",
+ " for i in range(len(entrada)):\n",
+ "\n",
+ " if entrada[i] == -1:\n",
+ " salida_A[i] = 1\n",
+ "\n",
+ " elif entrada[i] == 1:\n",
+ " salida_B[i] = 1\n",
+ "\n",
+ " return salida_A, salida_B"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#------------\n",
+ "# DINÃMICA\n",
+ "#------------\n",
+ "\n",
+ "@njit\n",
+ "def Dinamica(vec, l, Pr, Pj, Ps):\n",
+ "\n",
+ " #Vector que almacena las concentraciones de las partÃculas A y B\n",
+ " concentracion=np.copy(vec)\n",
+ " #Vector que almacena la distribución de espaciamientos\n",
+ " Pesp=np.zeros(len(concentracion)) \n",
+ "\n",
+ " # Números aleatorios:\n",
+ " #--------------------\n",
+ "\n",
+ " #Se tienen 4 números aleatorios. 3 de ellos representan la posibilidad de tener o no,\n",
+ " #reacción, incerción y salto. El restante es el sitio en el que ocurrirá la dinámica\n",
+ "\n",
+ " Escogida = rm.randint(0,2*l) #Num aleatorio que escoge el sitio de red para efectuar la dinpamica\n",
+ " particula = concentracion[Escogida] #Tipo de partÃcula escogida\n",
+ "\n",
+ " p_r = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de reacción\n",
+ " p_j = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de incerción\n",
+ " p_s = rm.uniform(0,1e7)/1e7 #Num aleatorio para la prob de salto\n",
+ "\n",
+ "\n",
+ " # Saltos\n",
+ " # -------\n",
+ "\n",
+ " #Para los saltos se tiene en cuenta lo siguiente:\n",
+ " #Si 0 < p_s < Ps -> ocurre un salto hacia la derecha\n",
+ " #Si Ps < p_s < 2*Ps -> ocurre un salto hacia la izquierda\n",
+ " #(Siempre y cuando el sitio al que saltará la partÃcula este vacÃo)\n",
+ " #Esto se realiza para que podamos tener información tanto de la posibilidad de salto,\n",
+ " #como de la posibilidad de no saltar\n",
+ "\n",
+ " #Ejemplo: Si Ps=0.5, entonces:\n",
+ " #Salto a la derecha : 0 < p_s < 0.5 --> 50% de las veces \n",
+ " #Salto a la izquierda: 0.5 < p_s < 1 --> 50% de las veces\n",
+ " #No salta --> 0% de las veces\n",
+ "\n",
+ " if p_s < 2*Ps:\n",
+ "\n",
+ " #Salto general (Se pueden escoger todos los sitios, excepto los extremos)\n",
+ " if Escogida != 0 and Escogida != 2*l and concentracion[Escogida] != 0:\n",
+ " #Salto derecha\n",
+ " if p_s < Ps and concentracion[Escogida + 1] == 0:\n",
+ " concentracion[Escogida + 1] = particula\n",
+ " concentracion[Escogida] = 0\n",
+ " #Salto izquierda\n",
+ " elif Ps < p_s < 2*Ps and concentracion[Escogida - 1] == 0:\n",
+ " concentracion[Escogida - 1] = particula\n",
+ " concentracion[Escogida] = 0\n",
+ "\n",
+ " #Extremo izquierdo (dominio de la partÃcula A) \n",
+ " elif Escogida == 0:\n",
+ " #Incercion: Si p_j < Pj, entonces se efectúa la incerción\n",
+ " if p_j < Pj and concentracion[Escogida] == 0:\n",
+ " concentracion[Escogida] = -1\n",
+ " #Salto: Se considera que hay salto a la derecha\n",
+ " elif p_s < Ps and concentracion[0] != 0 and concentracion[1] == 0:\n",
+ " concentracion[1] = particula\n",
+ " concentracion[0] = 0\n",
+ "\n",
+ " #Extremo derecho (dominio de la partÃcula B)\n",
+ " elif Escogida == 2*l:\n",
+ " #Incersion: Si p_j < Pj, entonces se efectúa la incerción\n",
+ " if p_j < Pj and concentracion[Escogida] == 0:\n",
+ " concentracion[Escogida] = 1\n",
+ " #Salto: Se considera que hay salto a la izquierda\n",
+ " elif p_s < Ps and concentracion[2*l] != 0 and concentracion[2*l-1] == 0:\n",
+ " concentracion[2*l-1] = particula\n",
+ " concentracion[2*l] = 0\n",
+ "\n",
+ " # Reaccion\n",
+ " #----------\n",
+ "\n",
+ " #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)\n",
+ "\n",
+ " try: \n",
+ " #El try y el except es puesto porque hay ocasiones en las que se tiene por ejemplo:\n",
+ " #[-1,-1,-1,0,0,0]\n",
+ " #Entonces en este caso no es posible encontrar un lmb\n",
+ "\n",
+ " #rma\n",
+ " where1=np.where(concentracion == -1) #devuelve una tupla\n",
+ " rmav=where1[0] #Arreglo de las ubicaciones de las partÃculas A\n",
+ " rma=rmav[-1] #PartÃcula más a la derecha\n",
+ "\n",
+ " #lmb\n",
+ " where2=np.where(concentracion == 1) \n",
+ " lmbv=where2[0] #Arreglo de las ubicaciones de las partÃculas B\n",
+ " lmb=lmbv[0] #PartÃcula más a la izquierda\n",
+ "\n",
+ " if (p_r < Pr) and (abs(rma-lmb) == 1): #Las partÃculas se aniquilan\n",
+ " concentracion[rma] = 0 \n",
+ " concentracion[lmb] = 0\n",
+ "\n",
+ " #Se encuentran los nuevos rma y lmb para encontrar la distr. de espaciamientos\n",
+ " #rma\n",
+ " where1=np.where(concentracion == -1)\n",
+ " rmav=where1[0]\n",
+ " rma=rmav[-1]\n",
+ "\n",
+ " #lmb\n",
+ " where2=np.where(concentracion == 1)\n",
+ " lmbv=where2[0]\n",
+ " lmb=lmbv[0]\n",
+ "\n",
+ " #Función que devuelve un arreglo de la distribución de espaciamientos\n",
+ " Pesp= P_i_xt(rma,lmb,l)\n",
+ "\n",
+ " except: \n",
+ " pass\n",
+ "\n",
+ " return concentracion, Pesp # Devuelve 2 arreglos"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# ---------\n",
+ "# CICLOS\n",
+ "# ---------\n",
+ "\n",
+ "#Repite la función 'Dinamica' una cantidad M de veces, determinada por 'time'\n",
+ "\n",
+ "@njit\n",
+ "def ciclos(c_ini, l, Pr, Pj, Ps, time):\n",
+ "\n",
+ " resultado=np.copy(c_ini) # copia de la condición inicial\n",
+ " salida_conc_A=np.zeros(len(resultado)) # arreglo de salida para la concentración de A\n",
+ " salida_conc_B=np.zeros(len(resultado)) # arreglo de salida para la concentración de B\n",
+ " salida_esp=np.zeros(len(resultado)) # arreglo de salida para la dist. de espaciamientos\n",
+ "\n",
+ " for i in range(time+1):\n",
+ " resultado, espaciamiento = Dinamica(resultado, l, Pr, Pj, Ps) #Dinamica para el tiempo i\n",
+ " A, B = concentracionA_B(resultado) # Separación de arreglos\n",
+ "\n",
+ " return A, B, espaciamiento # Devuelve 3 arreglos"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# ----------------\n",
+ "# REPETICIONES\n",
+ "# ----------------\n",
+ "\n",
+ "#Repite la función 'ciclos' una cantidad determinada (rep) de repeticiones \n",
+ "\n",
+ "pool = mp.Pool(mp.cpu_count())\n",
+ "\n",
+ "global results\n",
+ "results=[]\n",
+ "\n",
+ "def get_result(result):\n",
+ " results.append(result)\n",
+ "\n",
+ "#Uso de paralelización\n",
+ "def repeticiones(rep, c_ini, l, Pr, Pj, Ps, time):\n",
+ "\n",
+ " for i in prange(rep):\n",
+ " pool.apply_async(ciclos, args=(c_ini, l, Pr, Pj, Ps, time), callback=get_result)\n",
+ " \n",
+ " # Cada 1000 iteraciones sacamos un print en pantalla para conocer el estado de la simulacion\n",
+ " if not(i % 100):\n",
+ " print('Working in iteration number {} out of {}.'.format(i, rep))\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#--------------\n",
+ "# ANALISIS\n",
+ "#--------------\n",
+ "\n",
+ "def ANALISIS(L=5, dx=1., dt=1., D=0.5, j=1., K=1., tiempo=2000, rep=10, opcion='graficos', nombre = 'Data'):\n",
+ " \n",
+ " #PARÃMETROS\n",
+ " \n",
+ " beta= K/dx # Tasa de reacción\n",
+ " lamb = D/dx**2 #Tasa de salto\n",
+ " P_S = lamb * dt # Probabilidad de salto\n",
+ " P_J= j #Probabilidad de incersion\n",
+ " P_R= beta*dt #Probabilidad de reacción\n",
+ " \n",
+ " #CONDICION INICIAL\n",
+ " \n",
+ " # 0 : VacÃo \n",
+ " # -1 : PartÃcula A\n",
+ " # 1 : Particula B \n",
+ " \n",
+ " condicion_inicial = np.zeros(2*L+1)\n",
+ " \n",
+ " for i in range(2*L+1): #Ilustracion para L=3: [-1,-1,-1,1,1,1]\n",
+ " if i < L:\n",
+ " condicion_inicial[i] = -1\n",
+ " else:\n",
+ " condicion_inicial[i] = 1\n",
+ "\n",
+ "\n",
+ "#-----------------------------------OPCIONES------------------------------------\n",
+ " \n",
+ " if opcion=='graficos':\n",
+ "\n",
+ " x = np.linspace(-L, L, 2*L+1) #Posiciones\n",
+ " z = np.linspace(0, 2*L, 2*L+1) #Distancias dist. de espaciamientos\n",
+ "\n",
+ " salida_A=np.zeros(len(x))\n",
+ " salida_B=np.zeros(len(x))\n",
+ " salida_Desp=np.zeros(len(z))\n",
+ " \n",
+ " # Se almacena la información de TODA la dinámica en 3 arreglos de salida\n",
+ " # Se suman cada iteración de tiempo para posteriormente hacer su respectiva normalización y promedio\n",
+ " # y encontrar la densidad local contÃnua y dist. de espaciamientos contÃnua\n",
+ " \n",
+ " #Función que aplica la función ciclos para 'rep' repeticiones\n",
+ " repeticiones(rep, condicion_inicial, L, P_R, P_J, P_S, tiempo)\n",
+ " \n",
+ " #ESTO FALTARIA POR OPTIMIZAR!!!!!!!! ( HACE LENTO EL CÓDIGO )\n",
+ " pool.close()\n",
+ " pool.join()\n",
+ " \n",
+ " #La variable global results es la que almacena los resultados\n",
+ " for item in results:\n",
+ " A, B, D = item\n",
+ "\n",
+ " salida_A += A\n",
+ " salida_B += B\n",
+ " salida_Desp += D\n",
+ " \n",
+ " results.clear() #Se limpia la variable global\n",
+ "\n",
+ " #Normalización y cálculo de la densidad de concentraciones\n",
+ " salida_A = salida_A/(rep*2*L)\n",
+ " salida_B = salida_B/(rep*2*L)\n",
+ " salida_Desp = salida_Desp/salida_Desp.sum()\n",
+ "\n",
+ "\n",
+ " #--------------------------------------------------------\n",
+ " #------------------ARCHIVOS------------------------------\n",
+ " #--------------------------------------------------------\n",
+ "\n",
+ " #Escribe 3 archivos txt en la carpeta Archi:\n",
+ " # - Info para A\n",
+ " # - Info para B\n",
+ " # - Info para dist. de espaciamientos\n",
+ " \n",
+ " print('TERMINA LA SIMULACIÓN')\n",
+ " \n",
+ " #-----------------------------------------------\n",
+ " \n",
+ " info = 'Archi/'+str(nombre)+'_info.txt'\n",
+ " File_info = open(info, 'w')\n",
+ " \n",
+ " File_info.write('L: '+str(L)+'\\n')\n",
+ " File_info.write('D: '+str(D)+'\\n')\n",
+ " File_info.write('j: '+str(j)+'\\n')\n",
+ " File_info.write('K: '+str(K)+'\\n')\n",
+ " File_info.write('tiempo: '+str(tiempo)+'\\n')\n",
+ " File_info.write('Repetición: '+str(rep)+'\\n')\n",
+ " \n",
+ " #-----------------------------------------------\n",
+ "\n",
+ " nombre1 = 'Archi/'+str(nombre)+'A.txt'\n",
+ " nombre2 = 'Archi/'+str(nombre)+'B.txt'\n",
+ " nombre3 = 'Archi/'+str(nombre)+'Desp.txt'\n",
+ "\n",
+ " File_A = open(nombre1, 'w')\n",
+ " File_B = open(nombre2, 'w')\n",
+ " File_P_i = open(nombre3, 'w')\n",
+ "\n",
+ " for i in range(len(x)):\n",
+ " File_A.write(str(x[i])+' '+str(salida_A[i])+'\\n')\n",
+ " File_B.write(str(x[i])+' '+str(salida_B[i])+'\\n')\n",
+ " File_P_i.write(str(z[i])+' '+str(salida_Desp[i])+'\\n')\n",
+ "\n",
+ " File_A.close()\n",
+ " File_B.close()\n",
+ " File_P_i.close()\n",
+ "\n",
+ "\n",
+ " #----------------------------------------------------------\n",
+ " #----------------------GRAFICOS----------------------------\n",
+ " #----------------------------------------------------------\n",
+ "\n",
+ " #-------------------PERFILES DE DENSIDAD-------------------------\n",
+ " fig1=plt.figure(1)\n",
+ " plt.plot(x,salida_A,c='b')\n",
+ " plt.plot(x,salida_B,c='r')\n",
+ " plt.legend(['Particulas A', 'Particulas B'],loc='upper right',bbox_to_anchor=(1.35, 1))\n",
+ " plt.xlabel('x')\n",
+ " plt.ylabel('Perfil de densidad')\n",
+ " plt.show()\n",
+ "\n",
+ " #---------------DISTRIBUCIÓN DE ESPACIAMIENTOS--------------------\n",
+ " fig2=plt.figure(2)\n",
+ " plt.plot(z,salida_Desp,c='b', alpha=0.5)\n",
+ " plt.xlabel('x')\n",
+ " plt.ylabel('Distribución de espaciamientos')\n",
+ " plt.show()\n",
+ "\n",
+ " #-------------------------------------------------------------------\n",
+ " \n",
+ " #---------------\n",
+ " # Dinamica\n",
+ " #---------------\n",
+ " \n",
+ " #Imprime en cada paso de tiempo, el vector concentración\n",
+ " #Con esta opción se puede visualizar lo que ocurre en cada paso de tiempo\n",
+ " \n",
+ " elif opcion=='dinamica':\n",
+ " \n",
+ " salida_din=np.copy(condicion_inicial)\n",
+ " \n",
+ " for i in range(0,tiempo):\n",
+ " print('tiempo ', i)\n",
+ " salida_din, otro = Dinamica(salida_din, L, P_R, P_J, P_S)\n",
+ " print(salida_din)\n",
+ " print(' ')\n",
+ " \n",
+ " #--------------------------------------------------------------------\n",
+ " \n",
+ " #---------------\n",
+ " # Estacionario\n",
+ " #---------------\n",
+ " \n",
+ " #Comprueba si se llega al estado estacionario en el tiempo ingresado. Esto se hace comparando\n",
+ " #el sistema en el tiempo = 'tiempo' y tiempo = 'tiempo' + 1000\n",
+ " \n",
+ " elif opcion=='estacionario':\n",
+ "\n",
+ " pruebas=[]\n",
+ " \n",
+ " #Evolución hasta el tiempo = 'tiempo'\n",
+ " A_est, B_est, Desp_est = ciclos(condicion_inicial, L, P_R, P_J, P_S, tiempo)\n",
+ " \n",
+ " estacionario = B_est - A_est\n",
+ "\n",
+ " pruebas.append(estacionario)\n",
+ " print('Estado en la iteración '+ str(tiempo)+':')\n",
+ " print(estacionario)\n",
+ " print(' ')\n",
+ " \n",
+ " #Evolución hasta el tiempo = 'tiempo + 1000'\n",
+ " A_est2, B_est2, Desp_est2 = ciclos(estacionario, L, P_R, P_J, P_S, 1000)\n",
+ " \n",
+ " estacionario2 = B_est2 - A_est2\n",
+ " \n",
+ " pruebas.append(estacionario2)\n",
+ " print('Estado en la iteración '+ str(tiempo+1000)+':')\n",
+ " print(estacionario2)\n",
+ " print(' ')\n",
+ "\n",
+ " salida_estacionario = pruebas[0] == pruebas[1]\n",
+ " \n",
+ " num_verdadero = list(salida_estacionario).count(True) \n",
+ " \n",
+ " print('Estado estacionario: ' + str(len(salida_estacionario) == num_verdadero))\n",
+ "\n",
+ " return len(salida_estacionario) == num_verdadero"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.8.5"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
--
GitLab