计算二维插值积分时出错.比较 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:Users909404Anaconda3libsite-packagesscipyintegratequadpack.py", line 323, in quad
    points)

  File "C:Users909404Anaconda3libsite-packagesscipyintegratequadpack.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 的值并返回 NaN 以防 xl>=xu ?

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.

推荐答案

除非我省略 np.vectorize 装饰器，否则我无法重现您的错误.虽然设置 xl/xu 值一致确实给了我一个 ZeroDivisionError.

I can't reproduce your error unless I omit the np.vectorize decorator. Setting xl/xu values that coincide does give me a ZeroDivisionError though.

无论如何，没有什么能阻止您在更高级别的函数中检查 xu 与 xl 的值.这样你就可以完全跳过无意义数据点的集成并尽早返回 np.nan :

Anyway, there's nothing stopping you from checking the values of xu vs xl in your higher-level function. That way you can skip integration entirely for nonsensical data points and return np.nan early:

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) 以便于比较:

With these definitions I get (following np.set_printoptions(linewidth=200) for easier comparison:

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]])

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

You can see that the values perfectly agree with your linked image.

现在，我有一个坏消息和一个好消息.坏消息是，虽然 np.vectorize 提供了围绕使用数组输入调用标量集成函数的语法糖，但与本机 for 循环相比，它实际上并没有为您提供加速.好消息是你可以用对 np.exp 的调用替换对 mpmath.exp 的调用，你会更快地得到相同的结果:

Now, I've got bad news and good news. The bad news is that while np.vectorize provides syntactical sugar around calling your scalar integration function with array inputs, it won't actually give you speed-up compared to a native for loop. The good news is that you can replace the calls to mpmath.exp with calls to np.exp and you'll end up with the same result much faster:

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 倍.

So the two methods give the same result (exactly!), but the numpy version is almost 15 times faster.

这篇关于计算二维插值积分时出错.比较 numpy 数组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！