Ускорение кода с помощью numpy

Question

Я новичок в python, и у меня есть этот код для вычисления потенциала внутри поля 1x1 с использованием ряда Fourier, но часть его идет слишком медленно (отмечено в коде ниже) ,

Если кто-то может мне помочь, я подозреваю, что мог бы что-то сделать с библиотекой numpy, но я не знаком с этим.

import matplotlib.pyplot as plt 
import pylab 
import sys 
from matplotlib import rc 
rc('text', usetex=False) 
rc('font', family = 'serif') 

#One of the boundary conditions for the potential. 

def func1(x,n): 
    V_c = 1 
    V_0 = V_c * np.sin(n*np.pi*x) 
    return V_0*np.sin(n*np.pi*x) 

#To calculate the potential inside a box: 
def v(x,y): 
    n = 1; 
    sum = 0; 
    nmax = 20; 

    while n < nmax: 
     [C_n, err] = quad(func1, 0, 1, args=(n),); 
     sum = sum + 2*(C_n/np.sinh(np.pi*n)*np.sin(n*np.pi*x)*np.sinh(n*np.pi*y)); 
     n = n + 1; 


    return sum; 

def main(argv): 
    x_axis = np.linspace(0,1,100) 
    y_axis = np.linspace(0,1,100) 
    V_0 = np.zeros(100) 
    V_1 = np.zeros(100) 

    n = 4; 

    #Plotter for V0 = v_c * sin() x 

    for i in range(100): 
     V_0[i] = V_0_1(i/100, n) 

    plt.plot(x_axis, V_0) 
    plt.xlabel('x/L') 
    plt.ylabel('V_0') 
    plt.title('V_0(x) = sin(m*pi*x/L), n = 4') 
    plt.show() 

    #Plot for V_0 = V_c(1-(x-1/2)^4) 

    for i in range(100): 
     V_1[i] = V_0_2(i/100) 


    plt.figure() 

    plt.plot(x_axis, V_1) 
    plt.xlabel('x/L') 
    plt.ylabel('V_0') 
    plt.title('V_0(x) = 1- (x/L - 1/2)^4)') 
    #plt.legend() 
    plt.show() 

    #Plot V(x/L,y/L) on the boundary: 

    V_0_Y = np.zeros(100) 
    V_1_Y = np.zeros(100) 
    V_X_0 = np.zeros(100) 
    V_X_1 = np.zeros(100) 

    for i in range(100): 
     V_0_Y[i] = v(0, i/100) 
     V_1_Y[i] = v(1, i/100) 
     V_X_0[i] = v(i/100, 0) 
     V_X_1[i] = v(i/100, 1) 

    # V(x/L = 0, y/L): 

    plt.figure() 
    plt.plot(x_axis, V_0_Y) 
    plt.title('V(x/L = 0, y/L)') 
    plt.show() 

    # V(x/L = 1, y/L): 

    plt.figure() 
    plt.plot(x_axis, V_1_Y) 
    plt.title('V(x/L = 1, y/L)') 
    plt.show() 

    # V(x/L, y/L = 0): 

    plt.figure() 
    plt.plot(x_axis, V_X_0) 
    plt.title('V(x/L, y/L = 0)') 
    plt.show() 

    # V(x/L, y/L = 1): 

    plt.figure() 
    plt.plot(x_axis, V_X_1) 
    plt.title('V(x/L, y/L = 1)') 
    plt.show() 

    #Plot V(x,y) 
####### 
# This is where the code is way too slow, it takes like 10 minutes when n in v(x,y) is 20. 
####### 

    V = np.zeros(10000).reshape((100,100)) 
    for i in range(100): 
     for j in range(100): 
      V[i,j] = v(j/100, i/100) 

    plt.figure() 
    plt.contour(x_axis, y_axis, V, 50) 
    plt.savefig('V_1') 
    plt.show() 

if __name__ == "__main__": 
    main(sys.argv[1:]) 

Answer 1

Вы можете найти как использовать FFT/ДПФ в этом документе:

Discretized continuous Fourier transform with numpy

Кроме того, в отношении вашей V матрицы, есть много способов, чтобы улучшить скорость выполнения. Один из них - убедиться, что вы используете Python 3 или xrange() вместо range(), если вы все еще находитесь на Python 2. . Я обычно ставлю эти строки в моем коде Python, чтобы позволить ему работать равномерно кастрированный баран Я использую Python 3. или 2. *

# Don't want to generate huge lists in memory... use standard range for Python 3.* 
range = xrange if isinstance(range(2), 
          list) else range 

Тогда вместо повторного вычисления j/100 и i/100, вы можете предвычисления этих значений и помещаем их в массив; зная, что разделение намного дороже, чем умножение! Что-то вроде:

ratios = np.arange(100)/100 

V = np.zeros(10000).reshape((100,100)) 
j = 0 
while j < 100: 
    i = 0 
    while i < 100: 
     V[i,j] = v(values[j], values[i]) 
     i += 1 
    j += 1 

Ну, во всяком случае, это довольно косметическое и не спасет вашу жизнь; и вам все равно нужно вызвать функцию v() ...

Затем вы можете использовать переплетение:

http://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/weave.html

Или написать всю ваш чистый код вычисления/петлю в C, скомпилировать и сгенерировать модуль который вы можете вызвать из Python.

Answer 2

Вы должны смотреть на broadcasting трюков Numpy в и vectorization (нескольких ссылок, один из первых хороших ссылок, которые выскакивают от Matlab, но это так же, как это применимо к NumPy - может кто-нибудь порекомендовать мне хорошую ссылку Numpy в комментариях, что В будущем я могу указать других пользователей?).

То, что я видел в своем коде (как только вы удалите все ненужные биты, как участки и неиспользуемые функции), является то, что вы в основном делают это:

from __future__ import division 
from scipy.integrate import quad 
import numpy as np 
import matplotlib.pyplot as plt 

def func1(x,n): 
    return 1*np.sin(n*np.pi*x)**2 

def v(x,y): 
    n = 1; 
    sum = 0; 
    nmax = 20; 

    while n < nmax: 
     [C_n, err] = quad(func1, 0, 1, args=(n),); 
     sum = sum + 2*(C_n/np.sinh(np.pi*n)*np.sin(n*np.pi*x)*np.sinh(n*np.pi*y)); 
     n = n + 1; 
    return sum; 

def main(): 
    x_axis = np.linspace(0,1,100) 
    y_axis = np.linspace(0,1,100) 
####### 
# This is where the code is way too slow, it takes like 10 minutes when n in v(x,y) is 20. 
####### 

    V = np.zeros(10000).reshape((100,100)) 
    for i in range(100): 
     for j in range(100): 
      V[i,j] = v(j/100, i/100) 

    plt.figure() 
    plt.contour(x_axis, y_axis, V, 50) 
    plt.show() 

if __name__ == "__main__": 
    main() 

Если вы посмотрите внимательно (можно использовать профайлер тоже), вы увидите, что вы интегрируете свою функцию func1 (которую я переименую в integrand) примерно 20 раз для каждого элемента массива 100x100 V. Однако подынтегральное выражение не меняется! Таким образом, вы уже можете вывести его из своей петли. Если вы сделаете это, и использовать вещания трюки, вы могли бы в конечном итоге с чем-то вроде этого:

import numpy as np 
from scipy.integrate import quad 
import matplotlib.pyplot as plt 

def integrand(x,n): 
    return 1*np.sin(n*np.pi*x)**2 

sine_order = np.arange(1,20).reshape(-1,1,1) # Make an array along the third dimension 
integration_results = np.empty_like(sine_order, dtype=np.float) 
for enu, order in enumerate(sine_order): 
    integration_results[enu] = quad(integrand, 0, 1, args=(order,))[0] 

y,x = np.ogrid[0:1:.01, 0:1:.01] 

term = integration_results/np.sinh(np.pi * sine_order) * np.sin(sine_order * np.pi * x) * np.sinh(sine_order * np.pi * y) 
# This is the key: you have a 3D matrix here and with this summation, 
# you're basically squashing the entire 3D structure into a flat, 2D 
# representation. This 'squashing' is done by means of a sum. 
V = 2*np.sum(term, axis=0) 

x_axis = np.linspace(0,1,100) 
y_axis = np.linspace(0,1,100) 
plt.figure() 
plt.contour(x_axis, y_axis, V, 50) 
plt.show() 

, которая проходит менее чем за одну секунду в моей системе. Трансляция становится более понятной, если вы возьмете ручку & и вытяните векторы, которые вы «вещаете», как если бы вы строили здание, от базовых блоков Тетрис.

Эти две версии функционально одинаковы, но один полностью векторизован, в то время как другой использует python for-loops. Как новый пользователь для python и numpy, я определенно рекомендую прочитать основы вещания.Удачи!