Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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.
share|improve this question
add comment

2 Answers

up vote 3 down vote accepted

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)

share|improve this answer
1  
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. –  Aron Ahmadia Feb 6 '13 at 10:46
    
"Why do you want to generate DivisionByZero exceptions?" Because I didn't know any better. Thanks for the answer. Let me try it out. –  Willie Wong Feb 6 '13 at 11:14
1  
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. –  Willie Wong Feb 6 '13 at 11:37
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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