es.davy.ai

Preguntas y respuestas de programación confiables

¿Tienes una pregunta?

Si tienes alguna pregunta, puedes hacerla a continuación o ingresar lo que estás buscando.

  • es.davy.ai/
  • Programación /
  • El código Python tiene un gran cuello de botella, pero no tengo suficiente experiencia para identificar dónde se encuentra.

El código Python tiene un gran cuello de botella, pero no tengo suficiente experiencia para identificar dónde se encuentra.

26 Jul, 2023 Programación

Mi código pretende modelar la energía promedio de la desintegración alfa, funciona pero es muy lento.

import numpy as np
from numpy import sin, cos, arccos, pi, arange, fromiter
import matplotlib.pyplot as plt
from random import choices

r_cell, d, r, R, N = 5.5, 15.8, 7.9, 20, arange(1,10000, 50)

def total_decay(N):
    theta = 2*pi*np.random.rand(2,N)
    phi = arccos(2*np.random.rand(2,N)-1)

x = fromiter((r*sin(phi[0][i])*cos(theta[0][i]) for i in range(N)),float, count=-1)
dx = fromiter((x[i] + R*sin(phi[1][i])*cos(theta[1][i]) for i in range(N)), float,count=-1)
y = fromiter((r*sin(phi[0][i])*sin(theta[0][i]) for i in range(N)),float, count=-1)
dy = fromiter((y[i] + R*sin(phi[1][i])*sin(theta[1][i]) for i in range(N)),float,count=-1)
z = fromiter((r*cos(phi[0][i]) for i in range(N)),float, count=-1)
dz = fromiter((z[i] + R*cos(phi[1][i]) for i in range(N)),float, count=-1)

return x, y, z, dx, dy, dz

def inter(x,y,z,dx,dy,dz, N):
    intersections = 0

for i in range(N): #Checks to see if a line between two points intersects with the target cell
    a = (dx[i] - x[i])*(dx[i] - x[i]) + (dy[i] - y[i])*(dy[i] - y[i]) + (dz[i] - z[i])*(dz[i] - z[i])
    b = 2*((dx[i] - x[i])*(x[i]-d) + (dy[i] - y[i])*(y[i])+(dz[i] - z[i])*(z[i]))
    c = d*d + x[i]*x[i] + y[i]*y[i] + z[i]*z[i] - 2*(d*x[i]) - r_cell*r_cell
    if b*b - 4*a*c >= 0:
        intersections += 1
return intersections

def hits(N):
    I = []
    for i in range(len(N)):
        decay = total_decay(N[i])
        I.append(inter(decay[0],decay[1],decay[2],decay[3],decay[4],decay[5],N[i]))
    return I

def AE(I,N): 
    p1, p2 = 52.4 / (52.4 + 18.9), 18.9 / (52.4 + 18.9)
    E = [choices([5829.6, 5793.1], cum_weights=(p1,p2),k=1)[0] for _ in range(I)]
    return sum(E)/N

def list_AE(I,N):
    E = [AE(I[i],N[i]) for i in range(len(N))]
    return E

plt.plot(N, list_AE(hits(N),N))
plt.title('Energía promedio por dosis con respecto al número de desintegraciones')
plt.xlabel('Número de desintegraciones [N]')
plt.ylabel('Energía promedio [keV]')
plt.show()

¿Alguien con experiencia puede señalar dónde se produce el cuello de botella, explicar por qué sucede y cómo optimizarlo? Gracias de antemano.

Tags:  , ,

Answer

  1. Avatar for davy.ai
    davy.ai

    El cuello de botella del código ocurre en la función inter, donde un bucle itera sobre todos los N decaimientos y calcula la intersección de cada uno con la celda objetivo. Esta operación es O(N) y se realiza para cada valor de N en el rango especificado, lo que resulta en una desaceleración significativa a medida que N aumenta.

    Para optimizar esto, podemos aprovechar el hecho de que el cálculo de la intersección es independiente para cada decaimiento y utilizar operaciones vectorizadas para realizar el cálculo en lotes. Podemos extraer las coordenadas relevantes de x, y, y z para todos los decaimientos a la vez y calcular los elementos de la fórmula de intersección en forma vectorial utilizando Numpy. Luego, podemos vectorizar la decisión de si ocurre una intersección verificando si el discriminante es no negativo para cada decaimiento utilizando las operaciones de matriz de Numpy.

    Aquí está la función inter optimizada:

    def inter(x,y,z,dx,dy,dz, N):
    a = np.sum((dx - x)**2 + (dy - y)**2 + (dz - z)**2, axis=0)
    b = 2*((dx - x)*(x - d) + (dy - y)*y + (dz - z)*z)
    c = d*d + x*x + y*y + z*z - 2*d*x - r_cell*r_cell
    discriminant = b*b - 4*a*c
    intersecciones = np.sum(discriminante >= 0)
    return intersecciones

    Esta implementación aprovecha las operaciones vectorizadas en Numpy y reduce significativamente el tiempo de cálculo.

    Otra posible optimización es precalcular valores que no dependen de N, como 2*pi*np.random.rand(2,N) y arccos(2*np.random.rand(2,N)-1) fuera de la función total_decay, y pasarlos como argumentos. Esto evita los gastos generados repetidamente para cada valor de N.

    def total_decay(N, theta, phi):
    x = r*np.sin(phi[0])*np.cos(theta[0])
    dx = x + R*np.sin(phi[1])*np.cos(theta[1])
    y = r*np.sin(phi[0])*np.sin(theta[0])
    dy = y + R*np.sin(phi[1])*np.sin(theta[1])
    z = r*np.cos(phi[0])
    dz = z + R*np.cos(phi[1])
    return x, y, z, dx, dy, dz

theta = 2*pi*np.random.rand(2,10000)
phi = arccos(2*np.random.rand(2,10000)-1)

def hits(N):
    decay = total_decay(N, theta[:,:N], phi[:,:N])
    I = inter(decay[0],decay[1],decay[2],decay[3],decay[4],decay[5],N)
    return I

def list_AE(I,N):
    p1, p2 = 52.4 / (52.4 + 18.9), 18.9 / (52.4 + 18.9)
    E = [choices([5829.6, 5793.1], cum_weights=(p1,p2),k=1)[0] for _ in range(I)]
    return sum(E)/N

plt.plot(N, list_AE(hits(N),N))
plt.title('Energía promedio por dosis en función del número de decaimientos')
plt.xlabel('Número de decaimientos [N]')
plt.ylabel('Energía promedio [keV]')
plt.show()

    Esta implementación solo genera theta y phi una vez, fuera del bucle sobre N, lo que resulta en una mayor aceleración.

    Estas optimizaciones juntas deberían mejorar significativamente el rendimiento del código.

Comments are closed.