Как я могу сделать scipy.odeint быстрее? - PullRequest
Новая
       68

Как я могу сделать scipy.odeint быстрее?

0 голосов
/ 05 мая 2020

В настоящее время я решаю интегрированную систему из 559 нелинейных дифференциальных уравнений. Мне нужно подогнать полученные решения к некоторым экспериментальным данным, варьируя константы c1, c2 b и g. Я использую scipy.odeint, и я хотел бы знать, есть ли способ ускорить мою программу, поскольку для ее выполнения требуется много времени.

Код такой: 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import random as rd
from numba import jit

L=np.loadtxt('C:/Users/Pablo/Desktop/TFG/Probas/matriz_L_Pablo.txt')  
I=np.loadtxt('C:/Users/Pablo/Desktop/TFG/Probas/vector_I_Pablo.txt')



k=np.diag(L)
n=len(k)  #Contamos o numero de nodos

u=np.zeros(n)
for i in range (n):
    u[i]=rd.random()


M=np.zeros((n,n))
derivs=np.zeros(n)

c1=100 ; c2=10000 ; b=0.01 ; g=1


@jit
def f(y,t,params):
    suma=0
    c1,c2,b,g=params
    for i in range(n):
        for j in range(n):
            if i==j:
                M[i,i]=(1-y[i]/b)+g*(1-y[i])+c2*I[i]*(1/n-1)
            if i!=j:
                M[i,j]=(1/n)*(c1*L[i,j]+c2*I[i])
            out=(M[i,j]*y[j])
            suma=suma+out
        derivs[i]=suma
        suma=0
    return derivs





#Condicions iniciais

y0=u

#lista cos parametros
params=[c1,c2,b,g]

#tempos de int
tf=1
deltat=0.001
t=np.arange(0,tf,deltat)

#solucion

sol= odeint(f, y0,t, args=(params,))

(Извините, если не очень понятно, я здесь впервые)

1 Ответ

0 голосов
/ 07 мая 2020

Вы можете попробовать векторизовать свой код. Функция f делает 2 вещи: сначала создает матрицу M, а затем выполняет умножение $$ M y $$. Умножение $$ M y $$ легко векторизовать, потому что все, что нам нужно сделать, это использовать функцию matmul numpy. 

def f(y,t,params):
    suma=0
    c1,c2,b,g=params
    for i in range(n):
        for j in range(n):
            if i==j:
                M[i,i]=(1-y[i]/b)+g*(1-y[i])+c2*I[i]*(1/n-1)
            if i!=j:
                M[i,j]=(1/n)*(c1*L[i,j]+c2*I[i])

    return np.matmul(M,y)

Это должно немного помочь во время выполнения. Но наиболее трудоемкой частью является тот факт, что вся матрица M формируется каждый раз, когда вызывается f, и что она формируется по одному элементу за раз. Единственные части M, которые необходимо изменить при вызове f, - это части, которые зависят от y. Таким образом, все недиагональные записи в M могут быть заполнены до вызова решателя од. Таким образом, если M имеет размер (569x569), вместо того, чтобы вычислять все ~ 250000 + элементов M каждый раз, когда вызывается f, вам нужно будет вычислить только 569 элементов M по диагонали. Остальные 250000 записей M не зависят от y и указываются перед вызовом решателя од. Внесение этой модификации должно привести к огромному ускорению, так как это кажется основным узким местом в вашем коде. Наконец, вы также можете векторизовать заполнение диагонали M, используя что-то вроде numpy .diag_indices.

...