title: For loop seems faster than NumPy/SciPy 3D interpolation
body: <p>I feel confused about NumPy/SciPy interpolation methods. I implemented 3D linear interpolation with <code>LinearNDInterpolator</code> and I found it very slow. Then I wrote a brute-force triple for loop approach in pure Python, and surprisingly it give me 1000x speedup. I gave a shot for the Numba package as well, and it did not turn out to be any faster. </p>
  
  <p>By any source that I found on the internet, Python loops should be superslow compared to NumPy/SciPy and Numba. But this is not what I see.</p>
  
  <p>I post the whole source code that I run. I get these times on my machine:</p>
  
  <pre><code>Numpy ready:  3.94499993324  s,  result[0]=  0.480961746817
  Python for loop...
  Python ready:  0.0299999713898  s,  result[0]=  0.480961746817
  Numba for loop...
  Numba  0  ready:  0.223000049591  s,  result[0]=  0.480961746817
  Numba for loop...
  Numba  1  ready:  0.0360000133514  s,  result[0]=  0.480961746817
  </code></pre>
  
  <p>I use Anaconda Python 2.7. What am I missing here?</p>
  
  <pre><code>import numpy
  import scipy.interpolate
  import time
  from numba import jit
  
  
  # x: a (40,) numpy array of ordered ints
  # y: a (30,) numpy array of ordered ints
  # z: a (10,) numpy array of ordered ints
  # values: a (10,30,40) numpy array of floats
  # targetxs: a (NP,) numpy array of random floats
  # targetys: a (NP,) numpy array of random floats
  # targetzs: a (NP,) numpy array of random floats
  
  
  NP=1000
  
  def numpyInterp(x,y,z,values,targetxs,targetys,targetzs):
  
      start=time.time()
      zz, yy, xx = numpy.broadcast_arrays(z,y[:,numpy.newaxis],x[:,numpy.newaxis,numpy.newaxis])
      grid=numpy.reshape(numpy.array([zz,yy,xx]).swapaxes(1,3),(3,-1)).T
      values3D=numpy.reshape(values,-1)
  
      print 'Reshape matrix: ',time.time()-start
      start=time.time()
      f=scipy.interpolate.LinearNDInterpolator(grid,values3D)
      print 'Interpolation: ',time.time()-start
      #start=time.time()
      #result1=[f(targetzs[i],targetys[i],targetxs[i]) for i in range(len(targetzs))]
      #print 'Evaluation (list comprehension): ',time.time()-start
      # I found that map is slightly (not much) faster on my machine than list comprehension
      start=time.time()    
      result=numpy.squeeze(map(f,targetzs,targetys,targetxs))
      print 'Evaluation (map): ',time.time()-start    
      return result
  
  
  def pythonInterp(x,y,z,values,targetxs,targetys,targetzs):
  
      nx=len(x)
      ny=len(y)
      nz=len(z)
  
      ntarget=targetxs.shape[0]
  
      result=numpy.zeros((ntarget,))
  
      for targ in range(ntarget):
          westix=len(x)-2
          eastix=len(x)-1
          for ix in range(1,nx):
              if targetxs[targ] &lt;= x[ix]:
                  westix=ix-1
                  eastix=ix
                  break
  
          southiy=len(y)-2
          northiy=len(y)-1     
          for iy in range(1,ny):
              if targetys[targ] &lt;= y[iy]:
                  southiy=iy-1
                  northiy=iy
                  break
  
          upiz=len(z)-1
          downiz=len(z)-2
          for iz in range(1,nz):
              if targetzs[targ] &lt;= z[iz]:
                  downiz=iz-1
                  upiz=iz
                  break
  
          xratio=(targetxs[targ]-x[westix])/(x[eastix]-x[westix])
          yratio=(targetys[targ]-y[southiy])/(y[northiy]-y[southiy])
  
          lowerresult=values[downiz,southiy,westix]+(values[downiz,southiy,eastix]-values[downiz,southiy,westix])*xratio+(values[downiz,northiy,westix]-values[downiz,southiy,westix])*yratio+(values[downiz,northiy,eastix]-values[downiz,northiy,westix]-values[downiz,southiy,eastix]+values[downiz,southiy,westix])*xratio*yratio
          upperresult=values[upiz,southiy,westix]+(values[upiz,southiy,eastix]-values[upiz,southiy,westix])*xratio+(values[upiz,northiy,westix]-values[upiz,southiy,westix])*yratio+(values[upiz,northiy,eastix]-values[upiz,northiy,westix]-values[upiz,southiy,eastix]+values[upiz,southiy,westix])*xratio*yratio
  
          result[targ]=lowerresult+(upperresult-lowerresult)*(targetzs[targ]-z[downiz])/(z[upiz]-z[downiz])
  
      return result
  
  @jit
  def numbaInterp(x,y,z,values,targetxs,targetys,targetzs):
  
      nx=len(x)
      ny=len(y)
      nz=len(z)
  
      ntarget=targetxs.shape[0]
  
      result=numpy.zeros((ntarget,))
  
      for targ in range(ntarget):
          westix=len(x)-2
          eastix=len(x)-1
          for ix in range(1,nx):
              if targetxs[targ] &lt;= x[ix]:
                  westix=ix-1
                  eastix=ix
                  break
  
          southiy=len(y)-2
          northiy=len(y)-1     
          for iy in range(1,ny):
              if targetys[targ] &lt;= y[iy]:
                  southiy=iy-1
                  northiy=iy
                  break
  
          upiz=len(z)-1
          downiz=len(z)-2
          for iz in range(1,nz):
              if targetzs[targ] &lt;= z[iz]:
                  downiz=iz-1
                  upiz=iz
                  break
  
          xratio=(targetxs[targ]-x[westix])/(x[eastix]-x[westix])
          yratio=(targetys[targ]-y[southiy])/(y[northiy]-y[southiy])
  
          lowerresult=values[downiz,southiy,westix]+(values[downiz,southiy,eastix]-values[downiz,southiy,westix])*xratio+(values[downiz,northiy,westix]-values[downiz,southiy,westix])*yratio+(values[downiz,northiy,eastix]-values[downiz,northiy,westix]-values[downiz,southiy,eastix]+values[downiz,southiy,westix])*xratio*yratio
          upperresult=values[upiz,southiy,westix]+(values[upiz,southiy,eastix]-values[upiz,southiy,westix])*xratio+(values[upiz,northiy,westix]-values[upiz,southiy,westix])*yratio+(values[upiz,northiy,eastix]-values[upiz,northiy,westix]-values[upiz,southiy,eastix]+values[upiz,southiy,westix])*xratio*yratio
  
          result[targ]=lowerresult+(upperresult-lowerresult)*(targetzs[targ]-z[downiz])/(z[upiz]-z[downiz])
  
      return result
  
  
  
  
  
  
  # Declare input data grid coordinates
  z=numpy.arange(10000,100001,10000)      # 10
  y=numpy.arange(30,60) # 30
  x=numpy.arange(0,40)  # 40
  
  
  # Initialize values (pointwise sin)
  zz, yy, xx = numpy.broadcast_arrays(z,y[:,numpy.newaxis],x[:,numpy.newaxis,numpy.newaxis])
  grid=numpy.array([zz,yy,xx]).swapaxes(1,3)[0,:,:,:]
  values=numpy.sin(grid)
  
  # Initialize points for interpolation
  targetxs=numpy.random.random((NP,))*40
  targetys=numpy.random.random((NP,))*30+30
  targetzs=numpy.random.random((NP,))*90000+10000
  
  
  # Running functions
  start=time.time()
  print 'Numpy...'
  a=numpyInterp(x,y,z,values,targetxs,targetys,targetzs)
  print 'Numpy ready: ',time.time()-start,' s,  result[0]= ',a[0]
  
  
  start=time.time()
  print 'Python for loop...'
  a=pythonInterp(x,y,z,values,targetxs,targetys,targetzs)
  print 'Python ready: ',time.time()-start,' s,  result[0]= ',a[0]
  
  for i in range(5):
      start=time.time()
      print 'Numba for loop...'
      a=numbaInterp(x,y,z,values,targetxs,targetys,targetzs)
      print 'Numba ',i,' ready: ',time.time()-start,' s,  result[0]= ',a[0]
  </code></pre>