erreur relevée dans le code qui ne doit pas fonctionner selon l'évaluation paresseuse

Question

I ont le code ci-dessous dans le cadre d'une fonction:

px = x2 - x1 
py = y2 - y1 
pz = z2 - z1 

div = px*px + py*py  

u = ((x0 - x1) * px + (y0 - y1) * py)/div 

la ligne retourne u=RuntimeWarning: invalid value encountered in divide lors de son exécution. C'est parce que parfois la ligne div= renvoie zéro.

Cependant, si je réécris la ligne u= comme:

u = np.where(div != 0, ((x0 - x1) * px + (y0 - y1) * py)/div, 0) 

retourne toujours le même avertissement d'exécution.

Ce code génère les nombres désirés mais je pensais que la fonction np.where était paresseuse. Si ce n'est pas le cas, il y a des accélérations potentielles dans d'autres parties du code que j'ai écrites (d'où je pose cette question). Qu'est-ce qui me manque? Est-ce que la fonction np.where calcule à la fois les entrées 'Vrai' et 'Faux' et en sélectionne une en fonction du booléen?

Notez que ceci est la solution que j'ai fini avec:

np.seterr(invalid='ignore') 
u = np.where(div != 0, ((x0 - x1) * px + (y0 - y1) * py)/div, 0) 
np.seterr(invalid='warn') 

bien que cela fonctionne bien aussi:

u = np.zeros_like(div) 
divb = div != 0 
u[divb] = ((x0[divb] - x1[divb]) * px[divb] + 
     (y0[divb] - y1[divb]) * py[divb])/div[divb] 

(qui est un peu ce que je pensais np.where fait ...

Ces deux solutions sont à peu près à la même vitesse mais les deux sont plus lentes que la seule fonction np.where.

Toutes les explications/suggestions sont les bienvenues! Merci.

Answer 1

Il s'agit d'une division par 0 problème que les programmeurs ont travaillé depuis les premiers langages vectorisés (par exemple APL, MATLAB).

Une solution que je l'ai utilisé dans le passé est (sous condition) ajouter 1 au diviseur:

u = ((x0 - x1) * px + (y0 - y1) * py)/(div + (div==0)) 

Il ne fonctionne pas dans tous les cas, mais il pourrait dans ce domaine, car il semble le div sera 0 seulement si les deux px et py sont 0. Dans ce cas, le numérateur est également 0. 0/1 = 0.

ou clipsage div avec une faible valeur (cela peut être la solution la plus rapide):

..../np.maximum(div,1e-16) 

Une recherche rapide sur le SO pour numpy divide by zero relevai d'autres questions. Par exemple:

https://stackoverflow.com/a/26248892/901925 suggère d'utiliser errstate contexte pour désactiver l'avertissement:

with numpy.errstate(divide='ignore'): 
    u = .../div 

divide œuvres pour des cas comme 1/0, alors que invalid est nécessaire pour 0/0 cas. Mais ignore finit par mettre soit inf ou nan dans le tableau de retour.Donc, vous devez toujours tester (div==0) pour obtenir une valeur 0.

Bien que je préfère comme l'apparition de cette forme:

with np.errstate(invalid='ignore'): 
    u = np.where(div==0, 0, .../div) 

commentaire de Warren explique pourquoi where ne fonctionne pas. Les arguments sont évalués avant d'être passés à la fonction. L'évaluation paresseuse nécessite la coopération de l'interpréteur Python. Il fait normalement partie de la syntaxe (par exemple if, |, &).