les valeurs assignées à un tableau numpy ne sont pas toujours égales aux valeurs assignées

Question

La procédure suivante permet d'obtenir des valeurs qui ne correspondent pas toujours à celles qui ont été attribuées :

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

Je m'attendais à ce que cette procédure imprime toujours zéro. Cependant, il arrive que ce ne soit pas le cas ! Lorsque la sortie est non nulle, elle est de plusieurs milliers.

Je ne sais pas si je fais quelque chose de mal ou s'il y a un bug.

Je l'ai signalé comme un problème de Scipy aussi.

solution :

Pauli Virtanen : "Le bug est dans votre code : les vecteurs i,j,k peuvent contenir deux fois une paire de coordonnées donnée." ( https://github.com/scipy/scipy/issues/2581 )

jorgeca a posté une réponse similaire ci-dessous.

Gracias.

Answer 1

Pour résumer votre problème, ajouter v à un tableau de zéros, puis en soustrayant v ne donne pas toujours un tableau de zéros :

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

Ce n'est pas possible !

En fait, les indices i , j y k sont des tableaux numpy, et il utilise donc indexation fantaisiste . Que doit-il se passer lorsqu'un triplet d'indices est répété ? C'est-à-dire, lorsque le i[m] == i[n] , j[m] == j[n] y k[m] == k[n] pour certains m, n. Il s'avère que par une détail de la mise en œuvre (qui peut changer à tout moment), seul le dernier triplet (dans l'ordre C) d'indices sera attribué, et au final vector_field[i,j,k,:] ne contient pas vraiment v .

Answer 2

Le problème semble être causé par une erreur de ronde.

J'ai fait une comparaison en utilisant votre code 1000 fois en utilisant dtypes : float16 , float32 , float64 y float96 ; dans vector_field et compté le nombre de fois où il donne une réponse non nulle :

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