2017-02-11 131 views
0

我正在尝试在Python中编写一个函数,它将参数a作为输入列表并返回依赖于这些参数的向量列表u。然后,我希望能够使用标准矢量操作(例如标量产品)对u采取行动。另外,在以下(从this post第二应答启发)的例子中,该函数将返回两个初始矢量vw的组合:向量化函数返回Python中的向量列表

import numpy as np 
v = np.array([1,0,0]) 
w = np.array([0,1,0]) 
a = np.linspace(0,np.pi,4) 
def u(a): 
    u = np.cos(a)*v + np.sin(a)*w 
    return u 
uvec = np.vectorize(u, otypes=[np.ndarray])  
ufin = np.array(uvec(a).tolist()) 
print "ufin =", ufin 
print "ufin.v =", np.dot(ufin,v) 

这将返回:

ufin = [[ 1.00000000e+00 0.00000000e+00 0.00000000e+00] 
[ 5.00000000e-01 8.66025404e-01 0.00000000e+00] 
[ -5.00000000e-01 8.66025404e-01 0.00000000e+00] 
[ -1.00000000e+00 1.22464680e-16 0.00000000e+00]] 
ufin.v = [ 1. 0.5 -0.5 -1. ] 

这是我想获得自ufin表现得像一个向量。你能告诉我是否有其他更简单的方法来实现这一点?我需要编写一个代码,其中需要定义大量的向量和向量操作,并希望它尽可能紧凑。

预先感谢您!


编辑:

我发现基于上次的答案this post另一个(显然更紧凑)的解决方案。这个想法是将参数列表重塑为一个列数组,使得该函数的输出自动(不需要向量化)将一个向量列表作为二维数组返回。这是这样完成的:

import numpy as np 
from numpy.core.umath_tests import inner1d 

v = np.array([1,0,0]) 
w = np.array([0,1,0]) 
a = np.linspace(0,np.pi,4).reshape((4,1)) 
b = np.linspace(0,np.pi/2,4).reshape((4,1)) 
def u(a): 
    u = np.cos(a)*v + np.sin(a)*w 
    return u 
print "u(a) =",u(a) 
print "u(b) =",u(b) 
print "u(a).v =",np.dot(u(a),v) 
print "u(a)^v =",np.cross(u(a),v) 
# print "u(a).u(b) =",np.dot(u(a),u(b)) # does not work 
print "u(a).u(b) =",inner1d(u(a),u(b)) # works 
print "u(a)^u(b) =",np.cross(u(a),u(b)) 

这将返回:

u(a) = [[ 1.00000000e+00 0.00000000e+00 0.00000000e+00] 
[ 5.00000000e-01 8.66025404e-01 0.00000000e+00] 
[ -5.00000000e-01 8.66025404e-01 0.00000000e+00] 
[ -1.00000000e+00 1.22464680e-16 0.00000000e+00]] 
u(b) = [[ 1.00000000e+00 0.00000000e+00 0.00000000e+00] 
[ 8.66025404e-01 5.00000000e-01 0.00000000e+00] 
[ 5.00000000e-01 8.66025404e-01 0.00000000e+00] 
[ 6.12323400e-17 1.00000000e+00 0.00000000e+00]] 
u(a).v = [ 1. 0.5 -0.5 -1. ] 
u(a)^v = [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 -8.66025404e-01] 
[ 0.00000000e+00 0.00000000e+00 -8.66025404e-01] 
[ 0.00000000e+00 0.00000000e+00 -1.22464680e-16]] 
u(a).u(b) = [ 1.00000000e+00 8.66025404e-01 5.00000000e-01 6.12323400e-17] 
u(a)^u(b) = [[ 0.   0.   0.  ] 
[ 0.   0.  -0.5  ] 
[ 0.   0.  -0.8660254] 
[ 0.   0.  -1.  ]] 

它既是涉及初始(uv)和载体(u(a)u(b)的输出列出操作正确的行为)以及涉及两个输出列表的操作。唯一需要注意的是(输出矢量列表之间的操作)是必须使用inner1d函数代替标准np.dot,因为后者被解释为由于两个矩阵的大小不一致而无法完成的矩阵乘积。

回答

0

如果不反对获得的2D阵列,而不是物体的阵列(前者比后者对于大多数用途更方便,反正)则外积针对此问题有用:

>>> np.outer(np.cos(a), v) + np.outer(np.sin(a), w) 
array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00], 
     [ 5.00000000e-01, 8.66025404e-01, 0.00000000e+00], 
     [ -5.00000000e-01, 8.66025404e-01, 0.00000000e+00], 
     [ -1.00000000e+00, 1.22464680e-16, 0.00000000e+00]]) 
>>> _ @ v 
array([ 1. , 0.5, -0.5, -1. ]) 

@ infix操作符只适用于python 3.5或更高版本。在旧版本中,您可以使用np.dot。无论在这方面的含义是matrix multiplication

+0

谢谢你的有用答案。我使用更简洁的版本编辑了我的初始文章,该版本基于将参数列表重塑为列数组。你能否告诉我你认为最适合的方法? – user3450569

+0

有没有标准的方法(如函数库)用Python执行这些类型的操作?这种操作在处理物理问题时非常常见(例如,绘制轨道角动量作为轨道参数的函数),所以如果出现这种情况,我会感到惊讶。 – user3450569

+0

@ user3450569是的,在这种情况下(您的功能恰好是“向量就绪”的),您的解决方案会更加优雅。你也可以写'a = np.linspace(0,np.pi,4)[:,None]'来获得我个人觉得在眼睛上更容易的列向量;但这是一个品味问题。你可能想看看'np.ogrid'和'np.ix_'这些在类似的情况下很有用。 –

0

所以你告诉vectorze函数u返回一个数组,实际上是object,因为没有dtype=array

In [475]: uvec=np.vectorize(u, otypes=[object]) 
In [476]: uvec(a) 
Out[476]: 
array([array([ 1., 0., 0.]), array([ 0.5  , 0.8660254, 0.  ]), 
     array([-0.5  , 0.8660254, 0.  ]), 
     array([ -1.00000000e+00, 1.22464680e-16, 0.00000000e+00])], dtype=object) 

所以,如果你想“解开”这一点,创造你必须将其转换为阵列的列表,然后用np.array或(np.stack)加入他们沿着一个新的层面的二维数组:

In [477]: uvec(a).tolist() 
Out[477]: 
[array([ 1., 0., 0.]), 
array([ 0.5  , 0.8660254, 0.  ]), 
array([-0.5  , 0.8660254, 0.  ]), 
array([ -1.00000000e+00, 1.22464680e-16, 0.00000000e+00])] 

但是对于像avectorize这样的单个1d阵列的迭代已被过度杀死。一个简单的列表理解作品一样好

In [483]: np.array([u(i) for i in a]) 
Out[483]: 
array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00], 
     [ 5.00000000e-01, 8.66025404e-01, 0.00000000e+00], 
     [ -5.00000000e-01, 8.66025404e-01, 0.00000000e+00], 
     [ -1.00000000e+00, 1.22464680e-16, 0.00000000e+00]] 

vectorize是比较有用时,有几个输入数组,他们可以有几个方面。它负责广播和迭代多个层面。你也可以使用np.frompyfunc(u,1,1),它总是返回一个对象数组,并且速度稍快。

但你也可以计算这个不用循环:

def ufunc(a,v,w): 
    a = a[:,None] 
    return np.cos(a)*v + np.sin(a)*w 

In [498]: ufunc(a,v,w) 
Out[498]: 
array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00], 
     [ 5.00000000e-01, 8.66025404e-01, 0.00000000e+00], 
     [ -5.00000000e-01, 8.66025404e-01, 0.00000000e+00], 
     [ -1.00000000e+00, 1.22464680e-16, 0.00000000e+00]]) 

这原来a为(4,1)阵列,广播对(3,)v产生(4,3)阵列。


另一种方法是定义一个函数,标量(各地)

import math 
def u1(a,v,w): 
    return math.cos(a)*v+math.sin(a)*w 
u1v = np.vectorize(u1, otypes=[float]) 

In [505]: u1v(a[:,None], v, w) 
Out[505]: 
array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00], 
     [ 5.00000000e-01, 8.66025404e-01, 0.00000000e+00], 
     [ -5.00000000e-01, 8.66025404e-01, 0.00000000e+00], 
     [ -1.00000000e+00, 1.22464680e-16, 0.00000000e+00]]) 

这需要的vectorize充分广播能力的优势。我虽然它会比你的代码慢,因为它必须在两个维度上循环,但它实际上测试更快(尽管不如我的ufunc快)。尝试一半矢量化u函数可能无法为您提供任何速度。


你可以有也定义了原有的功能采取vw作为参数,并使用vectorizeexcluded参数来传递它们原样。我不认为它具有性能优势。

In [521]: def u2(a,v,w): 
    ...:  return np.cos(a)*v + np.sin(a)*w 
    ...: 
In [522]: u2vec=np.vectorize(u2, otypes=[object],excluded=[1,2]) 
In [523]: u2vec(a,v,w) 
Out[523]: 
array([array([ 1., 0., 0.]), array([ 0.5  , 0.8660254, 0.  ]), 
     array([-0.5  , 0.8660254, 0.  ]), 
     array([ -1.00000000e+00, 1.22464680e-16, 0.00000000e+00])], dtype=object)