计算 2D 插值积分时出现错误。比较 numpy 数组

Question

我的优化任务涉及计算以下积分并找到

xl

和

xu

的最佳值：

迭代花费的时间太长，因此我决定通过计算所有可能值

xl

和

xu

的积分来加快迭代速度，然后在优化过程中对计算值进行插值。

我写了以下函数：

def k_integrand(x, xl, xu):
    return((x**2)*mpmath.exp(x))/((xu - xl)*(mpmath.exp(x)-1)**2)
@np.vectorize   
def K(xl, xu):
    y, err = integrate.quad(k_integrand, xl, xu, args = (xl, xu))
    return y

以及两个相同的数组

grid_xl

和

grid_xu

，其值动态递增。

当我运行代码时，我得到这个：

K(grid_xl, grid_xu)
Traceback (most recent call last):

  File "<ipython-input-2-5b9df02f12b7>", line 1, in <module>
    K(grid_xl, grid_xu)

  File "C:/Users/909404/OneDrive/Работа/ZnS-FeS/Теплоемкость/Python/CD357/4 - Optimization CD357 interpolation.py", line 75, in K
    y, err = integrate.quad(k_integrand, xl, xu, args = (xl, xu))

  File "C:\Users\909404\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py", line 323, in quad
    points)

  File "C:\Users\909404\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py", line 372, 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()

我猜这是因为

xl

应该总是小于

xu

。 有没有办法比较

xl

和

xu

的值并在

xl>=xu

的情况下返回 NaN？

最后我想要这样的东西：

并具备使用插值的能力。

也许我选择了错误的方式？我将不胜感激任何帮助。

Answer 1

除非我省略

np.vectorize

装饰器，否则我无法重现您的错误。设置一致的

xl

/

xu

值确实会给我一个

ZeroDivisionError

。

无论如何，没有什么可以阻止您在更高级别的函数中检查

xu

与

xl

的值。这样您就可以完全跳过无意义数据点的集成并尽早返回

np.nan

:

import numpy as np
import mpmath
import scipy.integrate as integrate

def k_integrand(x, xl, xu):    
    return ((x**2)*mpmath.exp(x))/((xu - xl)*(mpmath.exp(x)-1)**2)

@np.vectorize   
def K(xl, xu):
    if xu <= xl:
        # don't even try to integrate
        return np.nan
    y, err = integrate.quad(k_integrand, xl, xu, args = (xl, xu))
    return y

grid_xl = np.linspace(0.1,1,10)        # shape (10,) ~ (1,10)
grid_xu = np.linspace(0.5,4,8)[:,None] # shape (8,1)

通过这些定义，我得到了（遵循

np.set_printoptions(linewidth=200)

以便于比较：

In [35]: K(grid_xl, grid_xu)
Out[35]: 
array([[0.99145351, 0.98925197, 0.98650808, 0.98322919,        nan,        nan,        nan,        nan,        nan,        nan],
       [0.97006703, 0.96656815, 0.96254363, 0.95800307, 0.95295785, 0.94742104, 0.94140733, 0.93493293, 0.9280154 ,        nan],
       [0.93730403, 0.93263063, 0.92745487, 0.92178832, 0.91564423, 0.90903747, 0.90198439, 0.89450271, 0.88661141, 0.87833062],
       [0.89565597, 0.88996696, 0.88380385, 0.87717991, 0.87010995, 0.8626103 , 0.85469862, 0.84639383, 0.83771595, 0.82868601],
       [0.84794429, 0.8414176 , 0.83444842, 0.82705134, 0.81924245, 0.81103915, 0.8024601 , 0.79352503, 0.7842547 , 0.77467065],
       [0.79692339, 0.78974   , 0.78214742, 0.77416128, 0.76579857, 0.75707746, 0.74801726, 0.73863822, 0.72896144, 0.71900874],
       [0.7449893 , 0.73732055, 0.7292762 , 0.72087263, 0.71212741, 0.70305921, 0.69368768, 0.68403329, 0.67411725, 0.66396132],
       [0.69402415, 0.68602325, 0.67767956, 0.66900991, 0.66003222, 0.65076537, 0.6412291 , 0.63144388, 0.62143077, 0.61121128]])

您可以看到这些值与您链接的图像完全一致。

现在，我有一个坏消息和一个好消息。坏消息是，虽然

np.vectorize

提供了使用数组输入调用标量积分函数的语法糖，但与本机 for 循环相比，它实际上不会为您带来加速。好消息是，您可以将对

mpmath.exp

的调用替换为对

np.exp

的调用，并且最终会更快地得到相同的结果：

def k_integrand_np(x, xl, xu):    
    return ((x**2)*np.exp(x))/((xu - xl)*(np.exp(x)-1)**2)

@np.vectorize   
def K_np(xl, xu):
    if xu <= xl:
        # don't even try to integrate
        return np.nan
    y, err = integrate.quad(k_integrand_np, xl, xu, args = (xl, xu))
    return y

有了这些定义

In [14]: res_mpmath = K(grid_xl, grid_xu)
    ...: res_np = K_np(grid_xl, grid_xu)
    ...: inds = ~np.isnan(res_mpmath)
    ...: 

In [15]: np.array_equal(res_mpmath[inds], res_np[inds])
Out[15]: True

In [16]: %timeit K(grid_xl, grid_xu)
107 ms ± 521 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [17]: %timeit K_np(grid_xl, grid_xu)
7.26 ms ± 157 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

因此这两种方法给出相同的结果（完全一样！），但 numpy 版本几乎快了 15 倍。

Answer 2

使用低级回调函数进行集成

以下答案是对 @Andras Deak 答案的评论，但太渴望评论了。

scipy集成函数多次调用k_integrand_np函数，这非常慢。使用纯 Python 函数的替代方法是使用低级回调函数。该函数可以使用 Numba 等编译器直接用 C 或 Python 编写。以下是此答案的稍微修改版本。

示例

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


##wrapper for a function that takes 3 input values
def jit_integrand_function(integrand_function):
    jitted_function = numba.njit(integrand_function)
    
    @cfunc(float64(intc, CPointer(float64)))
    def wrapped(n, xx):
        return jitted_function(xx[0], xx[1],xx[2])
    return LowLevelCallable(wrapped.ctypes)

#your function to integrate
def k_integrand_np(x, xl, xu):
  return ((x**2)*np.exp(x))/((xu - xl)*(np.exp(x)-1)**2)

#compile integrand
k_integrand_nb=jit_integrand_function(k_integrand_np)

#now we can use the low-level callable
@np.vectorize
def K_nb(xl, xu):
    if xu <= xl:
        # don't even try to integrate
        return np.nan
    y, err = integrate.quad(k_integrand_nb, xl, xu, args = (xl, xu))
    return y

#for comparison
@np.vectorize
def K_np(xl, xu):
    if xu <= xl:
        # don't even try to integrate
        return np.nan
    y, err = integrate.quad(k_integrand_np, xl, xu, args = (xl, xu))
    return y

性能

#create some data
grid_xl = np.linspace(0.1,1,500)      
grid_xu = np.linspace(0.5,4,800)[:,None] 

t1=time.time()
res_nb = K_nb(grid_xl, grid_xu)
print(time.time()-t1)
t1=time.time()
res_np = K_np(grid_xl, grid_xu)
print(time.time()-t1)

inds = ~np.isnan(res_nb)
np.allclose(res_nb[inds], res_np[inds])

K_np: 24.58s
K_nb: 0.97s (25x speedup)
allclose: True