scipy.integrate.quad上的限制阵列 [英] scipy.integrate.quad on arrays of limits

问题描述

从scipy.integrate四所需要的参数FUNC，A，B。其中，func为集成的函数，a和b是下和上积分限，分别。 a和b必须是数字。

我有一个情况我需要评估函数的积分成千上万不同的A，B，总结的结果。这需要很长的时间来通过循环。我想只给四数组a和b，希望四将返回相应的数组，但没有奏效。

下面是说明什么，我试图做的，同时与Python的循环工作，但速度很慢，和我在尝试量化不起作用一个code。如何解决这个问题有什么建议可以在快速（numpy的-IS）的方式来解决呢？

 导入numpy的是NP
从scipy.integrate进口四＃功能我需要整合：
DEF F（X）：
    返回np.exp（-x * X）＃不同限制了大名单：
的a_list = np.linspace（1，100，1E5）
b_list = np.linspace（2，200，1E5）＃慢环：
所有积分total_slow = 0＃总和。
一，拉链B（的a_list，b_list）：
    total_slow + =四路（F，A，B）[0]
        ＃（四返回一个元组，其中第一指标是结果，
        ＃所以[0]）＃矢量化的方法（不工作）：
total_fast = np.sum（四（楼的a_list，b_list）[0]）返回此错误： 线329，在_quad
    如果（B =天道酬勤和= -Inf！）：
ValueError错误：一个数组有超过真值
一个元件是不明确的。使用a.any（）或a.all（）

 

编辑：

实际的功能，我需要整合（ RHO ），还包含另外两个因素。 rho0 与相同长度的数组的a_list 和 b_list 。 ^ h 是一个标量。

 高清RHO（X，rho0，H）：
    返回rho0 * np.exp（ -  X * X /（2 * H * H））
 

EDIT2：

不同的解决方案的剖析。 'space_sylinder`是其中积分发生的功能。沃伦Weckesser的建议是尽可能快地通过一个简单的分析功能，通过阵列和比Python的慢速度环路（说明呼叫的次数，该计划甚至没有完成，它仍然使用657秒）〜500次。

  ncalls tottime percall cumtime percall文件名：LINENO（功能）
＃解析（但错误）逼近积分的解决方案：
      108 1.850 0.017 2.583 0.024 DensityMap.py:473(space_sylinder）
使用scipy.integrate.quad＃慢蟒蛇循环：
       69 19.223 0.279 657.647 9.531 DensityMap.py:474(space_sylinder）
＃矢量化scipy.special.erf（沃伦Weckesser的建议）：
      108 1.786 0.017 2.517 0.023 DensityMap.py:475(space_sylinder）
 

解决方案

的积分EXP（-x * X）是误差函数，所以你可以使用的 scipy.special.erf 来计算积分。鉴于标量 A 和 B ，你的功能从整体 A 到 b 是 0.5 * np.sqrt（np.pi）*（ERF（二） - ERF（A））。

ERF 是<一个href=\"http://docs.scipy.org/doc/numpy/reference/internals.$c$c-explanations.html#universal-functions\"相对=nofollow>ufunc，这意味着它处理数组参数。鉴于的a_list 和 b_list ，你的计算可以写成

 总= 0.5 * np.sqrt（np.pi）*（ERF（b_list） -  ERF（a_list）将）和（）
 

功能 RHO 也可以用 ERF 处理，通过适当的缩放：

  G = np.sqrt（2）* H
总= G * rho0 * 0.5 * np.sqrt（np.pi）*（ERF（b_list / G） -  ERF（的a_list / G））和（）
 

依靠前检查这对你慢的解决方案。对于一些价值观，在 ERF 功能减法会导致precision的显著损失。

quad from scipy.integrate needs the arguments func, a, b. Where func is a function to integrate and a and b is the lower and upper integration limits, respectively. a and b has to be numbers.

I have a situation where I need to evaluate the integral of a function for hundreds of thousands of different a, b and sum the results. This takes a long time to loop through. I tried to just give quad arrays for a and b, hoping quad would return the corresponding array, but that didn't work.

Here is a code illustrating what I'm trying to do, with both the Python loop that works, but is very slow, and my attempt at vectorizing that doesn't work. Any suggestions on how to this problem can be solved in a fast (numpy-is) way?

import numpy as np
from scipy.integrate import quad

# The function I need to integrate:
def f(x):
    return np.exp(-x*x)

# The large lists of different limits:
a_list = np.linspace(1, 100, 1e5)
b_list = np.linspace(2, 200, 1e5)

# Slow loop:
total_slow = 0  # Sum of all the integrals.
for a, b in zip(a_list, b_list):
    total_slow += quad(f, a, b)[0]
        # (quad returns a tuple where the first index is the result,
        # therefore the [0])

# Vectorized approach (which doesn't work):
total_fast = np.sum(quad(f, a_list, b_list)[0])

"""This error is returned:

 line 329, in _quad
    if (b != Inf and a != -Inf):
ValueError: The truth value of an array with more than 
one element is ambiguous. Use a.any() or a.all()
"""

EDIT:

The actual function I need to integrate (rho) also contains two other factors. rho0 is an array with the same length as a_list and b_list. H is a scalar.

def rho(x, rho0, H):
    return rho0 * np.exp(- x*x / (2*H*H))

EDIT2:

Profiling of different solutions. ´space_sylinder` is the function where the integral happens. Warren Weckesser suggestion is as fast as passing the arrays through a simple analytical function and ~500 times faster than the slow Python loop (note by the number of calls that the program didn't even finish and it still used 657 seconds).

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
# Analytic (but wrong) approximation to the solution of the integral:
      108    1.850    0.017    2.583    0.024 DensityMap.py:473(space_sylinder)
# Slow python loop using scipy.integrate.quad:
       69   19.223    0.279  657.647    9.531 DensityMap.py:474(space_sylinder)
# Vectorized scipy.special.erf (Warren Weckesser's suggestion):
      108    1.786    0.017    2.517    0.023 DensityMap.py:475(space_sylinder)

解决方案

The integral of exp(-x*x) is a scaled version of the error function, so you can use scipy.special.erf to compute the integral. Given scalars a and b, the integral of your function from a to b is 0.5*np.sqrt(np.pi)*(erf(b) - erf(a)).

erf is a "ufunc", which means it handles array arguments. Given a_list and b_list, your calculation can be written as

total = 0.5*np.sqrt(np.pi)*(erf(b_list) - erf(a_list)).sum()

The function rho can also be handled with erf, by using the appropriate scaling:

g = np.sqrt(2)*H
total = g*rho0*0.5*np.sqrt(np.pi)*(erf(b_list/g) - erf(a_list/g)).sum()

Check this against your slow solution before relying on it. For some values, the subtraction of the erf functions can result in a significant loss of precision.

这篇关于scipy.integrate.quad上的限制阵列的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！