CODEWITHSUNDEEP
pythonnumpyscipy
eldad-a
eldad-a

Reputation: 3191

values assigned to a numpy array do not always equal to the assigned values

The following procedure results in values which do not always match the assigned ones:

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

I expected this procedure to always print zero. However, sometimes it does not! When the output is non-zero, it is of several thousands.

I am not sure whether I am doing something wrong or whether there's some bug here.

I reported this as a scipy issue as well.

solution:

Pauli Virtanen: "The bug is in your code: the i,j,k vectors can contain a given coordinate pair twice." ( https://github.com/scipy/scipy/issues/2581 )

jorgeca posted a similar answer below.

Thanks!

Upvotes: 2

Views: 207

Answers (2)

jorgeca
jorgeca

Reputation: 5522

To summarize your problem, adding v to an array of zeros and then substracting v doesn't always yield an array of zeros:

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!!

That can't be right!

In fact, the indices i, j and k are numpy arrays, and so it's using fancy indexing. What should happen when a triplet of indices is repeated? That is, when the i[m] == i[n], j[m] == j[n] and k[m] == k[n] for some m, n. It turns out that by an implementation detail (that can change any time) only the last triplet (in C order) of indices will get assigned, and in the end vector_field[i,j,k,:] doesn't really contain v.

Upvotes: 2

Saullo G. P. Castro
Saullo G. P. Castro

Reputation: 58915

The problem seems to be caused by round error.

I did a comparison running your code 1000 times using dtypes: float16, float32, float64 and float96; in vector_field and counted the number of times it gives a non zero answer:

float16: 1000
float32: 1000
float64: 142
float96: 115

Upvotes: 1

Related Questions