0

I'm using Python2.7 to create a simple vector field and then plotting it...

But Jupyter complains about a division by 0 ("RuntimeWarning: divide by zero encountered in divide"), and I can't find it.

import numpy as np

def field_gen(x0, y0, x, y, q_cons = 1):
    dx = x0-x
    dy = y0-y
    dist = np.sqrt(np.square(dx)+np.square(dy))
    kmod = np.where( dist>0.00001, q_cons / dist, 0  ) 
    kdir = np.where( kmod != 0, (np.arctan2(-dy,-dx) * 180 / np.pi), 0)
    res_X = np.where( kmod !=0, kmod * (np.cos(kdir)) , 0 )
    res_Y = np.where( kmod !=0, kmod * (np.sin(kdir)) , 0 )
    return (res_X, res_Y)

n = 10
X, Y = np.mgrid[0:n, 0:n]

x0=2
y0=2

(u,v)= field_gen(x0, y0, X, Y)
#print(u) #debug
#print
#print(v)
plt.figure()
plt.quiver(X, Y, u, v, units='width')

Any hints on that?

mbauman
  • 30,958
  • 4
  • 88
  • 123
  • 1
    While `np.where` can be used to say "where this is true use this; else do that", it still will evaluate the entire of `q_cons / dist` first. `dist` will have a zero entry in it at `[2,2]`. Thus the error. – freethebees Mar 15 '19 at 17:20
  • Thanks @freethebees , I'm still new at Python. Are there any alternatives that doesn't evaluate the statements first? – Marcelo Silva Mar 15 '19 at 17:33
  • No worries! This is a classic Python gotcha. See my answer below. `np.divide` will take `dist > 0.001` as an input and use that to iterate through your array. Iterating manually would be slow in Python, but NumPy has some great compiled stuff that does that kind of thing real quick for you. – freethebees Mar 15 '19 at 17:39
  • Well, it works, thanks again! I was trying an ugly c-like approach... Using "dist = np.where( dist > 0, dist, np.inf )" inbetween. But, using your solution, in the case of [2,2], it gave back a 0, and I'm wondering if that's a default behavior from np.divide... docs says "values of False indicate to leave the value in the output alone", but since it wasn't provided an out array ("If not provided or None, a freshly-allocated array is returned"), it was allocated with zeros ? – Marcelo Silva Mar 15 '19 at 17:50
  • `out` allows you to define an array to put the result in. From what I understand, if that array was one filled with ones, then the cells that are `False` in `where` will be ones. In you don't provide `out`, then it creates a blank array and uses the values from `q_cons`. You'd have to play around a bit to be sure though. Glad I could help. – freethebees Mar 17 '19 at 08:00

1 Answers1

2

Don't be fooled into think np.where does all the work here. Python will still evaluate all the input arguments first, before running calling np.where.

So in your command kmod = np.where( dist>0.00001, q_cons / dist, 0 ), Python will evaluate dist>0.00001 (ok) and q_cons / dist (bad!) before running np.where.

Try np.divide instead. I think you want something like this:

np.divide(q_cons, dist, where=dist>0.00001 )
freethebees
  • 957
  • 10
  • 24