Fast Interpolation of 1D Vectors

QUESTION: Is there a faster way to interpolate 1D vectors than INTERPOL?

ANSWER: This answer is provided by Paul Ricchiazzi of the Institute for Computational Earth System Science at the University of California, Santa Barbara in a IDL newsgroup post on 21 February 1997. Here is his answer.

Editor's Note: Changes to INTERPOL in 1998 and 1999 have increased its speed by at least 100-fold, so that the FINDEX routine described below is no longer needed to achieve fast interpolation speeds.

"I find that using INTERPOL to do 1D interpolations on large irregular grids can be extremely time consuming. The problem with INTERPOL is that it uses a linear search to find where a given field point fits into the irregular grid. Below you'll find my solution to the problem. The procedure FINDEX uses a binary search to obtain a "floating point index" which can be used with INTERPOLATE. I have found that FINDEX, used with INTERPOLATE, can be up to 70 times faster than using INTERPOL by itself. I am donating this procedure to the IDL community in hopes of saving untold millions of machine cycles that would otherwise have been wasted in futile linear searches."

Here is the documentation header and FINDEX program Paul submitted. To download the lastest version of the file, which resides on the Earth Space Research Group's anonymous FTP server, click here.

   ;+
   ; ROUTINE:  findex
   ;
   ; PURPOSE:  Compute "floating point index" into a table using binary
   ;           search.  The resulting output may be used with INTERPOLATE.
   ;
   ; USEAGE:   result = findex(u,v)
   ;
   ; INPUT:    
   ;   u       a monitically increasing or decreasing 1-D grid
   ;   v       a scalor, or array of values
   ;
   ; OUTPUT:
   ;   result  Floating point index. Integer part of RESULT(i) gives
   ;           the index into to U such that V(i) is between
   ;           U(RESULT(i)) and U(RESULT(i)+1).  The fractional part
   ;           is the weighting factor
   ;
   ;                          V(i)-U(RESULT(i))
   ;                          ---------------------
   ;                     U(RESULT(i)+1)-U(RESULT(i))
   ;
   ;
   ; DISCUSSION: 
   ;           This routine is used to expedite one dimensional
   ;           interpolation on irregular 1-d grids.  Using this routine
   ;           with INTERPOLATE is much faster then IDL's INTERPOL
   ;           procedure because it uses a binary instead of linear
   ;           search algorithm.  The speedup is even more dramatic when
   ;           the same independent variable (V) and grid (U) are used
   ;           for several dependent variable interpolations.
   ;
   ;  
   ; EXAMPLE:  
   ;
   ;; In this example I found the FINDEX + INTERPOLATE combination
   ;; to be about 60 times faster then INTERPOL.
   ;
   ;  u=randomu(iseed,200000) & u=u(sort(u))
   ;  v=randomu(iseed,10)     & v=v(sort(v))
   ;  y=randomu(iseed,200000) & y=y(sort(y))
   ;
   ;  t=systime(1) & y1=interpolate(y,findex(u,v)) & print,systime(1)-t
   ;  t=systime(1) & y2=interpol(y,u,v)            & print,systime(1)-t
   ;  print,f='(3(a,10f7.4/))','findex:   ',y1,'interpol: ',y2,'diff:     ',y1-y2
   ;
   ; AUTHOR:   Paul Ricchiazzi                        21 Feb 97
   ;           Institute for Computational Earth System Science
   ;           University of California, Santa Barbara
   ;           paul@icess.ucsb.edu
   ;
   ; REVISIONS:
   ;
   ;-
   ;
   function findex,u,v
   nu=n_elements(u)
   nv=n_elements(v)

   us=u-shift(u,+1)
   us=us(1:*)
   umx=max(us,min=umn)
   if umx gt 0 and umn lt 0 then message,'u must be monotonic'
   if umx gt 0 then inc=1 else inc=0

   maxcomp=fix(alog(float(nu))/alog(2.)+.5) 

   ; maxcomp = maximum number of binary search iteratios

   jlim=lonarr(2,nv)
   jlim(0,*)=0          ; array of lower limits
   jlim(1,*)=nu-1       ; array of upper limits

   iter=0
   repeat begin
     jj=(jlim(0,*)+jlim(1,*))/2
     ii=where(v ge u(jj),n) & if n gt 0 then jlim(1-inc,ii)=jj(ii)
     ii=where(v lt u(jj),n) & if n gt 0 then jlim(inc,ii)=jj(ii)
     jdif=max(jlim(1,*)-jlim(0,*))
     if iter gt maxcomp then begin
       print,maxcomp,iter, jdif
       message,'binary search failed'
     endif
     iter=iter+1
   endrep until jdif eq 1 

   w=v-v
   w(*)=(v-u(jlim(0,*)))/(u(jlim(0,*)+1)-u(jlim(0,*)))+jlim(0,*)

   return,w
   end