python 中的多处理嵌套数值积分

Question

我在 python 中使用嵌套数值积分，其中每一层的限制取决于下一层。我的代码的整体结构看起来像

import numpy as np
import scipy.integrate as si

def func(x1, x2, x3, x4):
    return x1**2 - x2**3+x3*x2 - x4*x3**3  

def int1():
    """integrates `int2` over x1"""
    a1, b1 = -1, 3
    def int2(x1):
        """integrates `func` over x2 at given x1.""" 
        #partial_func1 = lambda x2: func(x1, x2)
        b2 = 1 - np.abs(x1)
        a2 = -np.abs(x1**3)
        def int3(x2):
            a3 = x2
            b3 = -a3
            def int4(x3):
                partial_func = lambda x4: func(x1, x2, x3, x4)
                a4 = 1+np.abs(x3)
                b4 = - a4
                return si.quad(partial_func,a4,b4)[0]
            return si.quad(int4, a3, b3)[0]
        return si.quad(int3, a2, b2)[0]     
    return si.quad(int2, a1, b1)[0]
result = int1()  # -22576720.048151683

在我的完整版代码中，积分和极限比较复杂，需要好几个小时运行，很不方便。每个积分似乎都可以很容易地并行化：似乎我应该能够使用多处理将集成分发到多个 CPU 并加快运行时间。

参考其他关于堆栈溢出的帖子，我尝试了以下方法：

def testfunc(intfunc,fmin,fmax):
    return scint.quad(intfun,fmin,fmax,epsabs=10**-40)[0]

result = pool.map(partial(partial(testfunc, intfunc = int4),fmin = a3),[b3])

但是我得到一个错误，本地对象不能被 pickle。

我遇到的另一个资源是 http://catherineh.github.io/programming/2016/10/04/parallel-integration-for-mere-mortals

但我需要一个函数，我也可以将限制作为输入传递（因此我使用了部分）。

有谁知道如何解决这些问题？我认为一个解决方案是可以处理多个输入的 pool.map 的某个版本会很棒，但是如果我对 partials 的使用有问题，那也很好找出来。

在此先致谢，如果这里有任何可以解决的问题，请告诉我！

Answer 1

这个答案可能并不令人满意，但希望它能让您深入了解问题属于哪个领域。

重申一下，原来的问题是计算四重积分

integrate(
    integrate(
        integrate(
            integrate(
                f(x1, x2, x3, x4),
                [1+abs(x3), -1-abs(x3)]
                ),
            [x2, -x2]
            ),
         [1-abs(x1), -x1**3]
         ),
    [-3, 1])

从数学上讲，可以将其表述为

integrate(f(x1, x2, x3, x4), Omega)

其中 Omega 是由上述积分限制定义的四维域。如果域是一维、二维或三维，那么您的问题的答案将很明确：

将复杂域离散化为线、三角形或四面体（这些分别是维度 1、2、3 中的单纯形）（使用 one of many mesh tools），然后
对每个 lines/triangles/tetrahedra（例如，来自 here）使用数字正交。

不幸的是，我不知道有任何工具可以将四维域离散为 4 单纯形，也不知道 4 单纯形的正交规则（顶点和中点规则可能除外）。但是，两者都可以在一般情况下创建；特别是一堆正交规则应该很容易想出。

为了完整起见，让我提一下，至少有一个 class 的域在任何维度上都存在积分规则：超立方体。

Answer 2

更新：

经过多次测试和重构，解决这个问题的最好方法似乎不是嵌套函数或定义，而是利用 scipy.integrate.quad 函数中的 args 参数来将外部变量传递给内部集成。

非常感谢评论的人！

python 中的多处理嵌套数值积分

Multiprocessing nested numerical integrals in python

python

parallel-processing

multithreading

scipy

numerical-integration