Python 的 1 元素元组和 SciPy 有问题

Trouble with Python's 1-element tuples and SciPy

我一直在尝试(取得一些成功)使用 numpy vectorize 函数编写矢量化集成调用,但时不时地我会遇到 Python 处理元组。

我想编写 integrate.quad 的变体,可以在点网格上集成向量值函数。类似地,我想创建一个 integrate.nquad 的版本,它在 n 维域上积分,并且可以在点网格上计算这些积分(即具有 n 维域的积分,矢量输出,沿着格子计算k维点space).

例如:

import numpy as np
from scipy import integrate

def vecint(F, I, *args):
    componentintegrals = [integrate.nquad(f, I, args) for f in F]
    retint = [CI[0] for CI in componentintegrals]
    if (len(retint)==1):
        retint = retint[0]
    reterr = np.sqrt(sum(CI[1]**2 for CI in componentintegrals))
    return retint, reterr

vecint 将多变量函数列表作为输入,并将其视为积分问题,其中对向量值函数进行积分。这段代码工作得很好,例如:

print(vecint([lambda x,y: np.sin(x), lambda x,y: np.cos(y)], [[0,np.pi],[0,np.pi]] ) )
print(vecint([lambda x: np.sin(x)], [[0,np.pi]]))
print(vecint([lambda x: np.cos(x), lambda x: np.sin(x)], [[0,np.pi]] ))
## and we can pass additional arguments. 
print(vecint( [lambda x,k: np.sin(x)+k], [[0,np.pi]], 1) )

以上所有调用都按预期工作。当我尝试向量化这些函数时,我的麻烦就开始了。向量化 integrate.quad 没问题...

def quad_1vz1(f, I, *args):
return np.vectorize(lambda n: integrate.quad(f, I[0], I[1], (n,)+args)[0])

符合预期。上面的代码允许这样的调用:

quad_1vz1(lambda x,k: np.sin(kx), [0,np.pi])(K)

其中K=np.mgrid[0:1:6j]等。这些是积分

$$\int_0^{\pi}\sin(kx) dx$$

对于 $k$ 的各种值。

当我尝试用上面的 vecint 函数替换 integrate.quad 时出现问题。例如:

## let's vectorize a 1-dimensional integral of a vector-valued function with one parameter. 
def vecint_2vz1(F, I, *args):
    #print(I, args)
    return np.vectorize(lambda n: vecint( F, I, (n,)+args )[0])

def f1(t,k):
    return np.cos(t)+k
def f2(t,k):
    return np.sin(t)+k

K = np.mgrid[0:1:6j]
print( vecint_2vz1( [f1,f2], [[0,np.pi]] )(K) )

上面的结果是 "ValueError: setting an array element with a sequence."

当在这里矢量化 vecint 时,K 的元素作为 1 元素元组的 1 元素元组发送,即额外的参数可能类似于 ((0,))。

我怀疑要避免这种情况,我必须做一些狡猾的 unpacking/repacking 论证....但我对 Python 的想法有点困惑。

似乎 Python 有时会自动将 1 元组转换为包含在其中的值...有时不会。这让我很困惑。我觉得我缺少一些基本的东西。

Python 的执行输出:

ValueError                                Traceback (most recent call last)
<ipython-input-131-8cb6f52e7d30> in <module>()
 17 ## integral of (cos(t)+k, sin(t)+k)dt for various k's.
 18 
---> 19 print( vecint_2vz1( [f1,f2], [[0,np.pi]] )(K) )

/usr/local/lib/python3.4/dist-packages/numpy/lib/function_base.py in __call__(self, *args, **kwargs)
1809             vargs.extend([kwargs[_n] for _n in names])
1810 
-> 1811         return self._vectorize_call(func=func, args=vargs)
1812 
1813     def _get_ufunc_and_otypes(self, func, args):

/usr/local/lib/python3.4/dist-packages/numpy/lib/function_base.py in _vectorize_call(self, func, args)
1882             if ufunc.nout == 1:
1883                 _res = array(outputs,
-> 1884                              copy=False, subok=True, dtype=otypes[0])
1885             else:
1886                 _res = tuple([array(_x, copy=False, subok=True, dtype=_t)

ValueError: setting an array element with a sequence.

此外,如果我在 vecint 中放置一个 litle print(args) 行,它会打印出诸如 ((a,),) 之类的项目,其中 a 是 mgrid K 的元素。

很难跟踪您的代码中发生了什么。我们正在尝试跟踪 quad 采取了什么,您的 vecint 做了什么,lambdasvectorize.

为什么 quad_1vz1 运行 而 quad_2vz1 不呢? quad 正在解包元组吗?我不知道。

我试图通过以下方式简化事情:

def foo2(*args):
    def foo(*args):
        print(args)
        return 1
    return np.vectorize(lambda n: foo((n,)+args))

产生:

In [148]: foo2()(np.arange(3))
((0,),)
((0,),)
((1,),)
((2,),)
Out[148]: array([1, 1, 1])
In [153]: foo2(3)(np.arange(2))
((0, 3),)
((0, 3),)
((1, 3),)
Out[153]: array([1, 1])

请注意,由于我没有在 vectorize 中指定输出类型,因此 运行 只需一步即可弄清楚输出是什么样的。因此我看到 ((0,),) 两次。如果初始类型与后来的类型不同(例如整数 v 浮点数),则会出现其他 SO 问题。

这2个元组层次是我自己(n,)+args的产物,也是*argsvectorize的产物。我必须做更多的实验来理清第二个两个的责任。

如果您只想迭代一个变量,

vectorize 是一个糟糕的工具。它不会增加速度,而且很难正确应用。如果您有多个变量要一起广播,它会更有用。

vectorize 有 2 个参数:

def foo2(*args):
    def foo(*args):
        print(args)
        return sum(*args)
    return np.vectorize(lambda x,y: foo((x,y)+args))

In [164]: foo2(10)(np.arange(2),np.arange(3,5)[:,None])
((0, 3, 10),)
((0, 3, 10),)
((1, 3, 10),)
((0, 4, 10),)
((1, 4, 10),)
Out[164]: 
array([[13, 14],
       [14, 15]])
In [166]: 10+np.arange(2)+np.arange(3,5)[:,None]
Out[166]: 
array([[13, 14],
       [14, 15]])

注意我把内打印改成了print(*args),我看到了(0, 3, 10)。使用 *argsargs 时需要小心 - 一个是元组,另一个不是(或者相反?)。


如果我的函数 returns 是一个数组,我得到你的错误:

def foo2(*args):
    def foo(*args):
        print(*args)
        return np.array(*args)
    return np.vectorize(lambda x,y: foo((x,y)+args))
   .....: 
In [197]: foo2(10)(np.arange(2),np.arange(3,5)[:,None])
...
ValueError: setting an array element with a sequence.

但一步就可以了

In [203]: foo2(3)(1,2)
(1, 2, 3)
(1, 2, 3)
Out[203]: array([1, 2, 3])

我可以指定 object otypes:

def foo2(*args):
    def foo(*args):
        print(*args)
        return np.array(*args)
    return np.vectorize(lambda x,y: foo((x,y)+args),otypes=[object])
   .....: 
In [206]: foo2(3)(1,2)
(1, 2, 3)
Out[206]: array([1, 2, 3], dtype=object)
In [207]: foo2(3)([1,1],[2,3])
(1, 2, 3)
(1, 3, 3)
Out[207]: array([array([1, 2, 3]), array([1, 3, 3])], dtype=object)
In [208]: foo2(10)(np.arange(2),np.arange(3,5)[:,None])
(0, 3, 10)
(1, 3, 10)
(0, 4, 10)
(1, 4, 10)
Out[208]: 
array([[array([ 0,  3, 10]), array([ 1,  3, 10])],
       [array([ 0,  4, 10]), array([ 1,  4, 10])]], dtype=object)

这些可以堆叠成一个数组,但在这个二维结果中我也必须展平(前提是内部数组的大小都相同)。

In [225]: np.vstack(np.ravel(A))
Out[225]: 
array([[ 0,  3, 10],
       [ 1,  3, 10],
       [ 0,  4, 10],
       [ 1,  4, 10]])

我依稀记得关于 vectorize 对象 otypes 的问题,但我不记得这是否存在问题(在某些版本中)或者它只是一个解决方案。