如何减少二维连接域上的集成时间 [英] How to reduce integration time for integration over 2D connected domains

问题描述

我需要在简单连接的域(大部分时间是凸的)上计算许多 2D 积分.我正在使用 python 函数 scipy.integrate.nquad 来做这个集成.但是，与矩形域上的积分相比，此操作所需的时间要长得多.有没有可能更快的实现?

I need to compute many 2D integrations over domains that are simply connected (and convex most of the time). I'm using python function scipy.integrate.nquad to do this integration. However, the time required by this operation is significantly large compared to integration over a rectangular domain. Is there any faster implementation possible?

这是一个例子；我首先在圆形域(使用函数内部的约束)上集成一个常数函数，然后在矩形域(nquad 函数的默认域)上积分.

Here is an example; I integrate a constant function first over a circular domain (using a constraint inside the function) and then on a rectangular domain (default domain of nquad function).

from scipy import integrate
import time

def circular(x,y,a):
  if x**2 + y**2 < a**2/4:
    return 1 
  else:
    return 0

def rectangular(x,y,a):
  return 1

a = 4
start = time.time()
result = integrate.nquad(circular, [[-a/2, a/2],[-a/2, a/2]], args=(a,))
now = time.time()
print(now-start)

start = time.time()
result = integrate.nquad(rectangular, [[-a/2, a/2],[-a/2, a/2]], args=(a,))
now = time.time()
print(now-start)

矩形域只需要 0.00029 秒，而圆形域需要 2.07061 秒来完成.

The rectangular domain takes only 0.00029 seconds, while the circular domain requires 2.07061 seconds to complete.

循环积分也给出以下警告:

Also the circular integration gives the following warning:

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze 
the integrand in order to determine the difficulties.  If the position of a 
local difficulty can be determined (singularity, discontinuity) one will 
probably gain from splitting up the interval and calling the integrator 
on the subranges.  Perhaps a special-purpose integrator should be used.
**opt)

推荐答案

使计算更快的一种方法是使用 numba，一个用于 Python 的即时编译器.

One way to make the calculation faster is to use numba, a just-in-time compiler for Python.

Numba 提供了一个 @jit 装饰器 编译一些 Python 代码并输出可以在多个 CPU 上并行运行的优化机器代码.Jitting 被积函数只需很少的努力，并且可以节省一些时间，因为代码经过优化以运行得更快.人们甚至不必担心类型，Numba 会在幕后完成所有这些工作.

Numba provides a @jit decorator to compile some Python code and output optimized machine code that can be run in parallel on several CPU. Jitting the integrand function only takes little effort and will achieve some time saving as the code is optimized to run faster. One doesn't even have to worry with types, Numba does all this under the hood.

from scipy import integrate
from numba import jit

@jit
def circular_jit(x, y, a):
    if x**2 + y**2 < a**2 / 4:
        return 1 
    else:
        return 0

a = 4
result = integrate.nquad(circular_jit, [[-a/2, a/2],[-a/2, a/2]], args=(a,))

这确实运行得更快，在我的机器上计时时，我得到:

This runs indeed faster and when timing it on my machine, I get:

 Original circular function: 1.599048376083374
 Jitted circular function: 0.8280022144317627

这减少了约 50% 的计算时间.

That is a ~50% reduction of computation time.

由于语言的性质，Python 中的函数调用非常耗时.与 C 等编译语言相比，开销有时会使 Python 代码变慢.

Function calls in Python are quite time consuming due to the nature of the language. The overhead can sometimes make Python code slow in comparison to compiled languages like C.

为了缓解这种情况，Scipy 提供了一个 LowLevelCallable 类，可用于提供对低级编译回调函数的访问.通过这种机制，可以绕过 Python 的函数调用开销，进一步节省时间.

In order to mitigate this, Scipy provides a LowLevelCallable class which can be used to provide access to a low-level compiled callback function. Through this mechanism, Python's function call overhead is bypassed and further time saving can be made.

注意，在nquad的情况下，传递给LowerLevelCallable的cfunc的签名必须是以下之一:

Note that in the case of nquad, the signature of the cfunc passed to LowerLevelCallable must be one of:

double func(int n, double *xx)
double func(int n, double *xx, void *user_data)

其中 int 是参数的数量，参数的值在第二个参数中.user_data 用于需要上下文操作的回调.

where the int is the number of arguments and the values for the arguments are in the second argument. user_data is used for callbacks that need context to operate.

因此，我们可以稍微更改 Python 中的循环函数签名以使其兼容.

We can therefore slightly change the circular function signature in Python to make it compatible.

from scipy import integrate, LowLevelCallable
from numba import cfunc
from numba.types import intc, CPointer, float64


@cfunc(float64(intc, CPointer(float64)))
def circular_cfunc(n, args):
    x, y, a = (args[0], args[1], args[2]) # Cannot do `(args[i] for i in range(n))` as `yield` is not supported
    if x**2 + y**2 < a**2/4:
        return 1 
    else:
        return 0

circular_LLC = LowLevelCallable(circular_cfunc.ctypes)

a = 4
result = integrate.nquad(circular_LLC, [[-a/2, a/2],[-a/2, a/2]], args=(a,))

用这个方法我得到

LowLevelCallable circular function: 0.07962369918823242

与原始版本相比减少了 95%，与函数的 jitted 版本相比减少了 90%.

This is a 95% reduction compared to the original and 90% when compared to the jitted version of the function.

为了使代码更整洁并保持被积函数的签名灵活，可以创建一个定制的装饰器函数.它将 jit 被积函数并将其包装成一个 LowLevelCallable 对象，然后可以与 nquad 一起使用.

In order to make the code more tidy and to keep the integrand function's signature flexible, a bespoke decorator function can be created. It will jit the integrand function and wrap it into a LowLevelCallable object that can then be used with nquad.

from scipy import integrate, LowLevelCallable
from numba import cfunc, jit
from numba.types import intc, CPointer, float64

def jit_integrand_function(integrand_function):
    jitted_function = jit(integrand_function, nopython=True)

    @cfunc(float64(intc, CPointer(float64)))
    def wrapped(n, xx):
        return jitted_function(xx[0], xx[1], xx[2])
    return LowLevelCallable(wrapped.ctypes)


@jit_integrand_function
def circular(x, y, a):
    if x**2 + y**2 < a**2 / 4:
        return 1
    else:
        return 0

a = 4
result = integrate.nquad(circular, [[-a/2, a/2],[-a/2, a/2]], args=(a,))

任意数量的参数

如果参数数量未知，那么我们可以使用方便的carray 函数，用于将 CPointer(float64) 转换为 Numpy 数组.

Arbitrary number of arguments

If the number of arguments is unknown, then we can use the convenient carray function provided by Numba to convert the CPointer(float64) to a Numpy array.

import numpy as np
from scipy import integrate, LowLevelCallable
from numba import cfunc, carray, jit
from numba.types import intc, CPointer, float64

def jit_integrand_function(integrand_function):
    jitted_function = jit(integrand_function, nopython=True)

    @cfunc(float64(intc, CPointer(float64)))
    def wrapped(n, xx):
        ar = carray(xx, n)
        return jitted_function(ar[0], ar[1], ar[2:])
    return LowLevelCallable(wrapped.ctypes)


@jit_integrand_function
def circular(x, y, a):
    if x**2 + y**2 < a[-1]**2 / 4:
        return 1
    else:
        return 0

ar = np.array([1, 2, 3, 4])
a = ar[-1]
result = integrate.nquad(circular, [[-a/2, a/2],[-a/2, a/2]], args=ar)

这篇关于如何减少二维连接域上的集成时间的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！