Слишком медленная векторизация python? Или это четырехъядерный? Или это мой код? - PullRequest
Новая
       11

Слишком медленная векторизация python? Или это четырехъядерный? Или это мой код?

0 голосов
/ 26 февраля 2020

Я новичок в сообществе и прошу прощения, если я не предоставил информацию, как предполагалось.

Я пытаюсь узнать python, пришедший из Matlab. У меня есть очень простой код для начальных целей: 

from numpy import vectorize
from scipy import integrate
from scipy.special import j1
from math import sqrt, exp, pi, log
import matplotlib.pyplot as plt
import numpy as np



def plot_hc(radius, pd):
    q = np.linspace(0.008, 1.0, num=500)
    y = hc(q, radius, pd)
    plt.loglog(q, y)
    plt.show()

def hc_formfactor(q, radius):
    y = (1.0 / q) * (radius * j1(q * radius))
    y = y ** 2
    return y

def g_distribution(z, radius, pd):
    return (1 / (sqrt(2 * pi) * pd)) * exp(
        -((z - radius) / (sqrt(
            2) * pd)) ** 2)


def ln_distribution(z, radius, pd):
    return (1 / (sqrt(2 * pi) * pd * z / radius)) * exp(
        -(log(z / radius) / (sqrt(2) * pd)) ** 2)

# Dist=1(for G_Distribution)
# Dist=2(for LN Distribution)
Dist = 1


@vectorize
def hc(x, radius, pd):
    global d
    if Dist == 1:
        nmpts = 4
        va = radius - nmpts * pd
        vb = radius + nmpts * pd
        if va < 0:
             va = 0
        d = integrate.quad(lambda z: g_distribution(z, radius, pd), va, vb)
    elif Dist == 2:
        nmpts = 4
        va = radius - nmpts * pd
        vb = radius + nmpts * pd
        if va < 0:
            va = 0
        d = integrate.quad(lambda z: ln_distribution(z, radius, pd), va, vb)
    else:
       d = 1

   def fun(z, x, radius, pd):
        if Dist == 1:
            return hc_formfactor(x, z) * g_distribution(z, radius, pd)
        elif Dist == 2:
            return hc_formfactor(x, z) * ln_distribution(z, radius, pd)
        else:
            return hc_formfactor(x, z)

    y = integrate.quad(lambda z: fun(z, x, radius, pd), va, vb)[0]
    return y/d[0]

if __name__ == '__main__':
plot_hc(radius=40, pd=0.5)

Как и некоторые, я должен использовать для l oop, но это еще больше снизило скорость. Код, используемый для l oop, выглядит следующим образом: 

from numpy import vectorize
from scipy import integrate
from scipy.special import j1
from math import sqrt, exp, pi, log
import matplotlib.pyplot as plt
import numpy as np



def plot_hc(radius, pd):
    q = np.linspace(0.008, 1.0, num=500)
    y = hc(q, radius, pd)
    plt.loglog(q, y)
    plt.show()

def hc_formfactor(q, radius):
    y = (1.0 / q) * (radius * j1(q * radius))
    y = y ** 2
    return y

def g_distribution(z, radius, pd):
    return (1 / (sqrt(2 * pi) * pd)) * exp(
        -((z - radius) / (sqrt(
            2) * pd)) ** 2)


def ln_distribution(z, radius, pd):
    return (1 / (sqrt(2 * pi) * pd * z / radius)) * exp(
        -(log(z / radius) / (sqrt(2) * pd)) ** 2)

# Dist=1(for G_Distribution)
# Dist=2(for LN Distribution)
Dist = 1


def hc(q, radius, pd):
    if Dist == 1:
        nmpts = 4
        va = radius - nmpts * pd
        vb = radius + nmpts * pd
        if va < 0:
            va = 0
        d = integrate.quad(lambda z: g_distribution(z, radius, pd), va,vb)
    elif Dist == 2:
        nmpts = 4
        va = radius - nmpts * pd
        vb = radius + nmpts * pd
        if va < 0:
            va = 0
        d = integrate.quad(lambda z: ln_distribution(z, radius, pd), va, vb)
    else:
        d = 1

    def fun(z, q, radius, pd):
        if Dist == 1:
            return hc_formfactor(q, z) * g_distribution(z, radius, pd)
        elif Dist == 2:
            return hc_formfactor(q, z) * ln_distribution(z, radius, pd)
        else:
            return hc_formfactor(q, z)

    y = empty([len(q)])
    for n in range(len(q)):
        y[n] = integrate.quad(lambda z: fun(z, q[n], radius, pd), va, vb)[0]

    return y / d[0]
if __name__ == '__main__':
plot_hc(radius=40, pd=0.5)

Если я запускаю ту же программу для тех же значений в Matlab, это очень быстро. Я не знаю, какую ошибку я сделал здесь. Пожалуйста, помогите:)

1 Ответ

3 голосов
/ 26 февраля 2020

Примечания

Функция vectorize предназначена в первую очередь для удобства, а не для производительности. Реализация по сути для l oop.

MATLAB имеет jit способностей к компиляции, numpy нет (numba предоставляет некоторые из них). Чтобы получить лучшую numpy скорость, вам нужно мыслить в терминах операций с целым массивом, как мы привыкли делать в старые времена MATLAB.

...