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

【问题标题】：Multiprocessing nested numerical integrals in pythonpython中的多处理嵌套数值积分
【发布时间】：2017-06-03 04:20:59
【问题描述】：

我正在使用 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])

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

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

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

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

提前感谢，如果这里有什么可以清理的，请告诉我！

【问题讨论】：

在我给出的示例代码中，每个内层的限制取决于外层的结果。如果它们是全局函数，则不能将限制作为浮点数处理，并且不能对积分进行数值计算
这将如何解决我提到的限制问题？就像在书面嵌套积分中一样，它们需要按顺序计算
我可能会过得很好，我可能会在这里学到一两件事:) 最终结果应该是-22576720.048151683?
这是我上次运行时的评估结果！
好的，我在这里看到了你的战斗 :) 我一直在努力重新排列你的代码以使其不被嵌套。能够多处理嵌套函数是一个巨大的进步，甚至可能是不可能的。另外，您的时间花在处理 scipy 上，这不应该有 GIL 限制。我不认为这是可以做到的，但我可能是错的。

标签： python multithreading parallel-processing scipy numerical-integration

【解决方案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），然后
在每条线/三角形/四面体上使用数值求积（例如，来自here）。

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

为了完整起见，让我提一下，在任何维度中至少有一类域存在集成规则：超立方体。

【讨论】：

【解决方案2】：

更新：

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

非常感谢评论的人！

【讨论】：