0

Buenas, ando liado con una matriz que se crea con un bucle for múltiple que a su vez dentro del bucle se integra a través de la función quad el tamaño de la matriz depende de unos valores de m y n. Para unos valores de m y n de 8 tarda aproximadamente unos 180 segundos (3 minutos) en crear la matriz. Algún consejo para optimizar el llenado de la matriz.

Un saludo.

M=N=8

def I00(x,i,m,a):
    I00=X(x,i,a)*X(x,m,a)
    return I00
K=np.zeros(shape=(m*n,m*n))
M=np.zeros(shape=(m*n,m*n))

mmm=np.linspace(0,m,num=m,endpoint=False,dtype='int')
ii=np.linspace(0,m,num=m,endpoint=False,dtype='int')
nnn=np.linspace(0,n,num=n,endpoint=False,dtype='int')
jj=np.linspace(0,n,num=n,endpoint=False,dtype='int')
mii=np.linspace(0,n*m,num=m,endpoint=False,dtype='int')
njj=np.linspace(0,n*m,num=n,endpoint=False,dtype='int')

for i,mi in zip(ii,mii):
    for j in jj:
        for mm,nj in zip(mmm,njj):
            for nn in nnn:

                K[mi+j][nj+nn]=(D[0][0]*integrate.quad(I22,0,a,args=(mm+1,i+1,a))[0]
                                       *integrate.quad(I00,0,b,args=(nn+1,j+1,b))[0]) +
                               (D[1][1]*integrate.quad(I00,0,a,args=(mm+1,i+1,a))[0]
                                       *integrate.quad(I22,0,b,args=(nn+1,j+1,b))[0]) +
                               (D[0][1]*((integrate.quad(I20,0,a,args=(mm+1,i+1,a))[0]
                                         *integrate.quad(I02,0,b,args=(nn+1,j+1,b))[0]) +
                                         (integrate.quad(I02,0,a,args=(mm+1,i+1,a))[0]
                                         *integrate.quad(I20,0,b,args=(nn+1,j+1,b))[0]))) +
                               (4*D[2][2]*integrate.quad(I11,0,a,args=(mm+1,i+1,a))[0]
                                         *integrate.quad(I11,0,b,args=(nn+1,j+1,b))[0]) +
                               (2*D[0][2]*((integrate.quad(I21,0,a,args=(mm+1,i+1,a))[0]
                                           *integrate.quad(I01,0,b,args=(nn+1,j+1,b))[0]) +
                                           (integrate.quad(I12,0,a,args=(mm+1,i+1,a))[0]
                                           *integrate.quad(I10,0,b,args=(nn+1,j+1,b))[0]))) +
                               (2*D[1][2]*((integrate.quad(I10,0,a,args=(mm+1,i+1,a))[0]
                                           *integrate.quad(I12,0,b,args=(nn+1,j+1,b))[0]) +
                                           (integrate.quad(I01,0,a,args=(mm+1,i+1,a))[0]
                                           *integrate.quad(I21,0,b,args=(nn+1,j+1,b))[0])))

                M[mi+j][nj+nn]=rho*h*integrate.quad(I00, 0, a, args=(mm+1,i+1,a))[0]
                                    *integrate.quad(I00, 0, b, args=(nn+1,j+1,b))[0]
1
  • Las funciones distintas de I00 son casi idénticas a estas, y todas son integrales sinoidales.
    – user21428
    Commented el 7 nov. 2016 a las 11:51

1 respuesta 1

2

Estás calculando reiteradamente las mismas integrales, una vez y otra. Hazlas por separado, sólo una vez, y guarda su valor para cuando lo necesites luego.

Por ejemplo, un esbozo sería:

M1 = { (i,mm): integrate.quad(I00, 0, a, args=(mm+1,i+1,a))[0]
           for i in ii
               for mm in mmm }
M2 = { (j,nn): integrate.quad(I00, 0, b, args=(nn+1,j+1,b))[0]
           for j in jj
               for nn in nnn }

for i,mi in zip(ii,mii):
    for j in jj:
        for mm,nj in zip(mmm,njj):
            for nn in nnn:

                ... # cálculo de K

                M[mi+j][nj+nn]=rho*h*M1[(i,mm)]*M2[(j,nn)]

Por concisión, sólo he puesto M, pero se puede hacer lo mismo con K expresándolo en función de M1 y M2.


Editado: "Compresiones de diccionarios"

Se pueden definir listas, diccionarios y conjuntos a través de lo que se conoce como "Compresiones" ("Comprehesions").

Por ejemplo, este código:

d = { k: inicia(k) for k in range(10) }

Es equivalente a hacer:

d = {}
for k in range(10):
    d[k] = inicia(k)

Ambos código generan un diccionario de 10 elementos con claves 0, 1, ..., 9 y con valores definido por la llamada a la función inicia.

La ventaja de la comprensión de diccionarios sería:

  • el intérprete sabe explícitamente cuáles son las variables que controlan la creación del diccionario. En lugar de añadir elemento por elemento, puede optar por añadir grupos de elementos y optimizar mejor.
  • se cumple (en general) el principio de Transparencia referencial, o lo que es lo mismo, puedes sustituir la expresión por el diccionario generado y el resto de código debería funcionar igual, lo que debería evitar los "efectos colaterales". Es sencillo aplicar algoritmos de concurrencia para acelerar el proceso.

En la forma clásica de anidar bucles for (que he puesto antes):

  • nada asegura que el diccionario esté vacío, ni siquiera tienes la certeza de que el objeto usado es un diccionario. (La experiencia te dice que es conveniente colocarlo antes de entrar al bucle para tenerlo a la vista).
  • el bucle for obliga a recorrerlo en un orden dado ya que se permite modificar la variable de control. No se cumple (en general) el principio de Transparencia referencial, el intérprete no es consciente de cuál será el producto final. No hay mucho más para optimizar que ir siguiendo el bucle hasta que se termine.

En resumen, al contar con más información sobre cómo se crea la secuencia, el intérprete es capaz de optimizar muchísimo más e, incluso, podría establecer mecanismos para aprovechar el multiproceso (Apache Spark, Hadoop, ...) y acelerar más aún el proceso de inicialización (Dejo este tema para otra ocasión).


Editado: Concurrencia

Creo que será mejor añadirte la idea de cómo inicializar un diccionario mediante Futures, introducido en Python3.5.

La idea es lanzar las integrales para que se ejecuten en un pool de procesos, de modo que luego sea fácil combinar todas las respuestas obtenidas para formar el diccionario.

Como no tengo tu código, meto el cálculo de factorial en la inicialización para que se note cómo arden las CPUs.

from concurrent import futures
from functools import reduce
from operator import mul

def factorial(n):
    return reduce(mul, range(2,n+1), 1)

def ini(k):
    n = factorial(15000)
    return (k, str(k))

if __name__ == "__main__":

    with futures.ProcessPoolExecutor(40) as ex:
        d = {k:v for (k,v) in ex.map(ini, range(1000), chunksize=10)}

    print(d)

Cuando se ejecuta este código desde la línea de comandos, se ejecuta la sección principal if __name__ ..., que creará un contexto de ejecución de 40 workers. A este pool van llegando los procesos para calcular ini() sobre cada elemento del iterable (range(1000)). El chunksize=10 indica que se envíen los procesos de 10 en 10, en lugar de ir uno a uno.

Lo he puesto así por más claridad, pero la creación del diccionario se puede hacer más directa:

with futures.ProcessPoolExecutor(40) as ex:
    d = dict(ex.map(ini, range(1000), chunksize=10))

Si lo adaptas a las integrales, deberías notar bastante mejora de velocidad al aprovechar toda la potencia de tu ordenador.

3
  • No conozco nada de scipy, pero estoy seguro que habrá algo más directo para crearte una sparse matrix con el mismo uso que el diccionario indicado. Commented el 7 nov. 2016 a las 16:19
  • Muchisimas gracias por la respuesta.
    – user21428
    Commented el 7 nov. 2016 a las 19:33
  • Por cierto podrias explicarme brevemente como funciona la estructura {():....for for} ? Un saludo
    – user21428
    Commented el 7 nov. 2016 a las 20:57

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.