I want to populate a numpy array with values from the smooth bump function

f(x) = exp ( - 1 / (1 - x^2) )     if |x| < 1,  f(x) =  0 otherwise

Currently I have something that works (as in gives me the right numbers on my platform)

x = linspace(-1.1 , 1.1, 300)   #Sample 300 points between [-1.1,1.1]
bump = exp( 1 - 1 / (1 - clip(square(x), 0,1)) ) 

when the absolute value of an entry in x is at least 1, its square gets clipped to 1, and we have "1/(1-1) = 1/0 = +inf" as "expected" on my platform, which then gets set by "exp(1 - inf) = 0" which is exactly the behaviour I want.

My questions:

  1. I suspect that the above is not the best practice. Am I correct in my suspicions?
  2. Are there better ways of handling this division by zero? At the end of the day the array x may not be just simply a linear list of values. So I want something that can compute f(x) from x efficiently.

Why do you want to generate DivisionByZero exceptions?

I would use masked arrays:

import numpy as np
x= np.linspace(-1.1,1.1,300)
masked_idx = (np.abs(x)>1)
masked_x = np.ma.array(x,mask=idx)

def f(x):
   return np.exp(-1.0/(1.0-x**2))

masked_f = f(masked_x)

plot(masked_x,masked_f)    # in IPython / pyplot

If you want, you can do the masking in your function (by having boundary arguments)

  • 1
    $\begingroup$ shakes fist at GertVdE for beating me to the answer. Gert has it, use array masks, don't rely on dangerous boundaries of the floating-point standard. For example, your code didn't work on my OS X.8 Intel Mac. $\endgroup$ Feb 6 '13 at 10:46
  • $\begingroup$ "Why do you want to generate DivisionByZero exceptions?" Because I didn't know any better. Thanks for the answer. Let me try it out. $\endgroup$ Feb 6 '13 at 11:14
  • 1
    $\begingroup$ Okay, so after I do the arithmetic on the non-zero values, I can just just use the filled() function to set the masked values to 0. That works great. Thanks again. $\endgroup$ Feb 6 '13 at 11:37

Another take:

_f = lambda x: np.exp(-1.0/(1.0-x*x))
f = lambda x: np.piecewise(x, [np.abs(x) < 1, np.abs(x) >= 1], [_f, 0.0])

x = np.linspace(-1.1,1.1,300)
bump = f(x)

my point here is that what you really need for efficiency is a ufunc, possibly implemented in C. In the absence of a true ufunc you can go with any numpy trick.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.