实现此2D数值积分的计算上更快的方式是什么? [英] What would be the computationally faster way to implement this 2D numerical integration?

问题描述

我对进行2D数值积分感兴趣.现在我正在使用scipy.integrate.dblquad，但是它非常慢.请参见下面的代码.我需要用完全不同的参数来评估这100次积分.因此，我想使处理尽可能快和高效.代码是:

I am interested in doing a 2D numerical integration. Right now I am using the scipy.integrate.dblquad but it is very slow. Please see the code below. My need is to evaluate this integral 100s of times with completely different parameters. Hence I want to make the processing as fast and efficient as possible. The code is:

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def f(q, z, t):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)


y = np.empty([len(q)])
for n in range(len(q)):
    y[n] = integrate.dblquad(lambda t, z: f(q[n], z, t), 0, 50, lambda z: 10, lambda z: 60)[0]

end = time.time()
print(end - start)

花费的时间是

212.96751403808594

太多了.请提出一种更好的方法来实现我的目标.在来到这里之前，我曾尝试过进行搜索，但没有找到任何解决方案.我读过quadpy可以更好，更快地完成这项工作，但我不知道如何实现相同的目的.请帮忙.

This is too much. Please suggest a better way to achieve what I want to do. I tried to do some search before coming here, but didn't find any solution. I have read quadpy can do this job better and very faster but I have no idea how to implement the same. Please help.

推荐答案

您可以使用Numba或可调用的低级

几乎是您的示例

我只是将函数直接传递给scipy.integrate.dblquad，而不是使用lambda生成函数的方法.

I simply pass function directly to scipy.integrate.dblquad instead of your method using lambdas to generate functions.

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def f(t, z, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):
    return 10.

def upper_inner(z):
    return 60.


y = np.empty(len(q))
for n in range(len(q)):
    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)
#143.73969149589539

这已经快了一点点(143 vs. 151s)，但是唯一的用途是有一个简单的示例进行优化.

This is already a tiny bit faster (143 vs. 151s) but the only use is to have a simple example to optimize.

仅使用Numba编译功能

要使其运行，还需要 Numba 和

To get this to run you need additionally Numba and numba-scipy. The purpose of numba-scipy is to provide wrapped functions from scipy.special.

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
import numba as nb

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

#error_model="numpy" -> Don't check for division by zero
@nb.njit(error_model="numpy",fastmath=True)
def f(t, z, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):
    return 10.

def upper_inner(z):
    return 60.


y = np.empty(len(q))
for n in range(len(q)):
    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)
#8.636585235595703

使用低级别的可呼叫电话

scipy.integrate函数还提供了传递C回调函数而不是Python函数的可能性.这些函数可以用例如C，Cython或Numba编写，我在本示例中使用它们.主要优点是，在函数调用中无需Python解释程序交互.

The scipy.integrate functions also provide the possibility to pass C-callback function instead of a Python function. These functions can be written for example in C, Cython or Numba, which I use in this example. The main advantage is, that no Python interpreter interaction is necessary on function call.

出色的 answer @Jacques Gaudin的显示了一种简单的方法，包括附加参数.

An excellent answer of @Jacques Gaudin shows an easy way to do this including additional arguments.

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
import numba as nb
from numba import cfunc
from numba.types import intc, CPointer, float64
from scipy import LowLevelCallable

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def jit_integrand_function(integrand_function):
    jitted_function = nb.njit(integrand_function, nopython=True)

    #error_model="numpy" -> Don't check for division by zero
    @cfunc(float64(intc, CPointer(float64)),error_model="numpy",fastmath=True)
    def wrapped(n, xx):
        ar = nb.carray(xx, n)
        return jitted_function(ar[0], ar[1], ar[2])
    return LowLevelCallable(wrapped.ctypes)

@jit_integrand_function
def f(t, z, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):
    return 10.

def upper_inner(z):
    return 60.


y = np.empty(len(q))
for n in range(len(q)):
    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)
#3.2645838260650635

这篇关于实现此2D数值积分的计算上更快的方式是什么?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！