计算2D插值的积分时出错.比较numpy数组 [英] Error in calculating integral for 2D interpolation. Comparing numpy arrays
问题描述
我的优化任务处理以下积分的计算，并找到xl和xu的最佳值:
My optimization task deals with calculation of the following integral and finding the best values of xl and xu:
迭代花费的时间太长，因此我决定通过计算所有可能值xl和xu的积分来加快它们的速度，然后在优化过程中对计算出的值进行插值.
Iterations take too long, so I decided to speed them up by calculating integral for all possible values xl and xu and then interpolate calculated values during optimization.
我编写了以下函数:
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，其值具有动态增量.
and two identical arrays grid_xl and grid_xu with dinamic increment of values.
运行代码时，我得到以下信息:
When I run the code I get this:
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?
I guess it comes from the fact that xl should be always less than xu.
Is there any way to compare the values of xl and xu and return NaN in case if xl>=xu?
最后，我想要这样的东西:
In the end I want to have something like this:
并具有使用插值的能力.
And to have the ability to use interpolation.
也许我选择了错误的方式?我将不胜感激.
Maybe I have chosen the wrong way? I'd appreciate any help.
推荐答案
使用低级回调函数进行集成
以下答案是对@Andras Deak答案的评论，但是渴望发表评论.
Using a low-level callback function for integration
The following answer is meant as a comment to @Andras Deak answer, but is to long for a comment.
scipy集成函数多次调用k_integrand_np函数，这非常慢.使用纯Python函数的替代方法是使用低级回调功能.可以使用Numba之类的编译器直接用C或Python编写此函数.以下是对该 answer.
The scipy integrate functions call the k_integrand_np function multiple times, which is quite slow. The alternative of using a pure Python function is to use a low-level callback function. This function can be written directly in C or in Python using a compiler like Numba. The following is a slightly modified version of this answer.
示例
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
这篇关于计算2D插值的积分时出错.比较numpy数组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
