CODEWITHSUNDEEP
pythonnumpyscipy
Dirklinux
Dirklinux

Reputation: 1155

scipy.weave.inline not working as expected with math library

I have a problem using scipy.weave.inline. I want to program a unitstep function centered around lcenter, and with a with of width_nm. I have two versions: The python version, called pm and an optimized function, called pm_weave, but it looks like abs is not working properly. See the code below. If you run it, you'll get a window of size 1 for the weave variety, no matter what the input is, so it looks like abs doesn't work. If you remove the abs, for example, it works exactly like you expect

How can I fix this?

def pm_weave(w,lcenter,width_nm):
  """ Return a unitstep function that is centered around lcenter with height 1.0 and width width_nm """
  lcenter = float(lcenter)
  w = float(w)
  width_nm = float(width_nm)
  code = """
  #include <math.h>
  float wl = 1.88495559215387594307758602997E3/w;

  if(abs(lcenter-wl) < width_nm) {
     return_val = 1.0;
     }
     else {
     return_val = 0.0;
     }
  """
  res = weave.inline(code,['w','lcenter','width_nm'],
                     type_converters = weave.converters.blitz,
                     compiler = "gcc", headers=["<math.h>"]
                    )
  return res



def pm(w,lcenter,width_nm):
  """ 
  Return a unitstep function centered around lcenter [nm] with width width_nm. w
  should be a radial frequency. 
  """
  return abs(600*np.pi/w - lcenter) < width_nm/2. and 1. or 0.



plot(wavelength_list,map(lambda w:pm(toRadialFrequency(w),778,1),wavelength_list),label="Desired behaviour")
plot(wavelength_list,map(lambda w:pm_weave(toRadialFrequency(w),778,1),wavelength_list),'^',label="weave.inline behaviour")
ylim(0,1.5)
show()

Upvotes: 2

Views: 504

Answers (1)

ely
ely

Reputation: 77504

I think you might need to use fabs() instead of abs() in the C code. abs() will truncate the result, while fabs() will work for floating-point arithmetic.

Upvotes: 2

Related Questions