Это пример кода, который я выполняю, который генерирует матрицу измерения (size x size). Матрица отправляется в FFT (через несколько итераций), и их норма является требуемым результатом. Поскольку это тестовый прогон, я установил size = 256 и итерации (zaxis) в 3. В настоящее время для обработки матрицы требуется 1-2 минуты.
Фактический производственный цикл требует: матрица 512 x 512, 1024 x 1024 (или, может быть, больше) с примерно 25 итерациями каждая, и мне интересно, смогу ли я ускорить этот скрипт на python.
Короче говоря , я генерирую сложную матрицу => присваивать ненулевые значения элемент за элементом в цикле => отправлять в FFT => вычислить norm = > сохранить норму в массиве. Работает отлично!
Тяжелая работа выполняется в следующем фрагменте кода, где ненулевые значения вычисляются как val. Здесь двумерный интеграл вычисляется отдельно для вещественной и мнимой составляющих. В идеале я должен иметь возможность выполнить это на нескольких ядрах . (* хотя я думаю , если назначение различных ненулевых матричных элементов может быть полностью выгружено на несколько ядер, это было бы очень эффективно. У меня нет опыта в многопроцессорной обработке. Спецификация системы: 1700X AMD, 8 ядер, 32 ГБ ОЗУ под управлением Python3, Win 10; в качестве альтернативы доступна система Ubuntu с 12 ядрами, 64 ГБ ОЗУ)
if (r < (dim/2) ):
c1 = fy(a,r,x,y)
c2 = r*complexAmp(a,b,x,y)
re = I_real ( r ) # double integral, real
im = I_imag ( r ) # double integral, imaginary
val=c1 * c2 *(re[0]+1j*im[0])
Итак, мой вопрос. Есть ли хорошие способы повысить скорость таких операций (и, надеюсь, я смогу узнать больше об эффективном программировании на python). Сейчас я проверяю ray и multiprocessing .
Ниже приведен полный сценарий ввода . Вывод показан внизу.
import time
import math
import cmath as cm
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import dblquad
#-----------------------------------------------------------
# DEFINE FUNCTIONS
def fy(a,b,x,y):
return (a*b**3+(x/2.5)+y)/50
def complexAmp(a,b,x,y):
return ( (cm.exp(-1j*x*y*cm.sqrt(a+b)))/ cm.sqrt( a ) ) *b
def wrap(r, rho, phi):
return cm.cos(phi)*cm.exp(-1j*2*math.pi*cm.sqrt(rho**2 \
+ r**2))/cm.sqrt(rho**2 + r**2)
def wrap_real(r, rho, phi):
res = cm.cos(phi)*cm.exp(-1j*2*math.pi*cm.sqrt(rho**2 +\
r**2))/cm.sqrt(rho**2 + r**2)
return res.real
def wrap_imag(r, rho, phi):
res = cm.cos(phi)*cm.exp(-1j*2*math.pi*cm.sqrt(rho**2 + \
r**2))/cm.sqrt(rho**2 + r**2)
return res.imag
rMax = 5
def I_real (value ):
return dblquad(lambda rho, phi: wrap_real (value, rho, phi) \
, 0, rMax, lambda x: 0, lambda x: 2*math.pi)
def I_imag (value ):
return dblquad(lambda rho, phi: wrap_imag (value, rho, phi) ,\
0, rMax, lambda x: 0, lambda x: 2*math.pi)
#-----------------------------------------------------------
# TEST INTEGRATION
print("\n-----------COMPLEX INTEGRATION RESULT ----------")
print (I_real ( 6 ), I_imag ( 6 ))
print("--------------------------------------------------")
# parameters governing size and step of grid
size=256
depth=10
step=1.75 # step of grid
n2 = 1.45
theta = math.asin(1.00025/n2)
# complex matrix to keep data
inp = np.zeros((size,size) , dtype=complex )
zaxis = np.arange(-60, -10, 20)
result = np.zeros(zaxis.shape[0])
n2 = 1.454
theta = math.asin(1.00025/n2) # update theta
dim = 16000
# The main program -----------------------------------------------
for z in range(zaxis.shape[0]):
print ("In the loop {0}".format(z))
start = time.time()
for i in range(inp.shape[0]):
for j in range(inp.shape[1]):
x = step*(i-(size/2))
y = step*(j-(size/2))
r = x**2 + y**2
#print(i, (i-(size/2)),j, (j-(size/2)) )
b = r*( math.sin( 1.00025 /n2)) *math.sqrt(2*r**2)
a = 250/abs(zaxis[z]-r)
rMax = abs(zaxis[z])*math.tan(theta)
val=0
if (r < (dim/2) ):
c1 = fy(a,r,x,y)
c2 = r*complexAmp(a,b,x,y)
re = I_real ( r ) # double integral, real
im = I_imag ( r ) # double integral, imaginary
val=c1 * c2 *(re[0]+1j*im[0])
inp [i][j] = val # substitue the value to matrix
end = time.time()
print("Time taken : {0} sec \n" . format( round(end - start,7 )))
b = np.fft.fft2(inp)
result [z] = np.linalg.norm(b)
Выход:
-----------COMPLEX INTEGRATION RESULT ----------
(-0.0003079405888916291, 1.0879642638692853e-17) (-0.0007321233659418995, 2.5866160149768244e-17)
--------------------------------------------------
In the loop 0
Time taken : 138.8842542 sec
[plot]
In the loop 1
Time taken : 134.3815458 sec
[plot]
In the loop 2
Time taken : 56.848331 sec
[plot]
[plot]