2

以下过程导致的值并不总是与分配的值匹配：

from scipy.interpolate import splprep, splev, splrep
import numpy as np

pos2indx = lambda vec: vec.round().astype(np.int64)

t = np.linspace(1,3,150)
x = 150+100*np.sin(t) + 5*np.random.randn(len(t))
y = 150+100*np.cos(t) + 5*np.random.randn(len(t))
z = 150+100*np.cos(t)*np.sin(t) + 5*np.random.randn(len(t))

vector_field = np.zeros((x.max()+10,y.max()+10,z.max()+10,3), dtype=np.float64)

out = splprep([x,y,z],u=t,k=3, full_output=1, quiet=1)
tck, t = out[0]
v = np.transpose(splev(t,tck, der=1))
i,j,k = pos2indx(x),pos2indx(y),pos2indx(z)

vector_field[i,j,k,:] += v
print np.sum(np.abs(vector_field[i,j,k,:]-v))

我希望这个过程总是打印零。但是，有时它不会！当输出为非零时，为数千。

我不确定我是否做错了什么或者这里是否有一些错误。

我也将此报告为 scipy 问题。

解决方案：

Pauli Virtanen：“错误在您的代码中：i,j,k 向量可以包含两次给定的坐标对。” （https://github.com/scipy/scipy/issues/2581）

jorgeca 在下面发布了类似的答案。

谢谢！

4

2 回答 2

2

总结您的问题，添加v到一个零数组然后减去v并不总是产生一个零数组：

vector_field = np.zeros((x.max()+10,y.max()+10,z.max()+10,3), dtype=np.float64)
vector_field[i,j,k,:] += v
print np.sum(np.abs(vector_field[i,j,k,:]-v))  # sometimes not 0!!

这不可能！

事实上， indexi和j是knumpy 数组，所以它使用了花哨的 indexing。当重复三个索引时会发生什么？也就是说，当i[m] == i[n],j[m] == j[n]并且k[m] == k[n]对于某些 m, n。事实证明，通过实现细节（可以随时更改），只有索引的最后一个三元组（按 C 顺序）将被分配，并且最终vector_field[i,j,k,:]并不真正包含v.

于 2013-06-18T16:52:46.020 回答
1

该问题似乎是由圆形错误引起的。

我1000使用 dtypes: 、 和 ; 对运行您的代码时间float16进行float32了float64比较float96。并vector_field计算它给出非零答案的次数：

float16: 1000
float32: 1000
float64: 142
float96: 115
于 2013-06-18T16:50:33.573 回答