How to make scipy interpolate give an extrapolated result beyond the input range

Question

I m trying to port a program which uses a hand-rolled interpolator  developed by a mathematician colleage  over to use the interpolators provided by scipy  I d like to use or wrap the scipy interpolator so that it has as close as possible behavior to the old interpolator   A key difference between the two functions is that in our original interpolator - if the input value is above or below the input range  our original interpolator will extrapolate the result  If you try this with the scipy interpolator it raises a ValueError  Consider this program as an example   import numpy as np from scipy import interpolate  x   np arange 0 10  y   np exp -x 3 0  f   interpolate interp1d x  y   print f 9  print f 11    Causes ValueError  because it s greater than max x    Is there a sensible way to make it so that instead of crashing  the final line will simply do a linear extrapolate  continuing the gradients defined by the first and last two points to infinity    Note  that in the real software I m not actually using the exp function - that s here for illustration only

User · Answer

The below code gives you the simple extrapolation module  k is the value to which the data set y has to be extrapolated based on the data set x  The numpy module is required    def extrapol k x y           xm np mean x           ym np mean y           sumnr 0          sumdr 0          length len x           for i in range 0 length               sumnr sumnr   x i -xm   y i -ym                sumdr sumdr   x i -xm   x i -xm             m sumnr sumdr          c ym- m xm           return  m k  c

User · Answer

Standard interpolate   linear extrapolate       def interpola v  x  y           if v  lt   x 0               return y 0   y 1 -y 0    x 1 -x 0    v-x 0           elif v  gt   x -1               return y -2   y -1 -y -2    x -1 -x -2    v-x -2           else              f   interp1d x  y  kind  cubic                return f v

User · Answer

As of SciPy version 0 17 0  there is a new option for scipy interpolate interp1d that allows extrapolation  Simply set fill value  extrapolate  in the call  Modifying your code in this way gives   import numpy as np from scipy import interpolate  x   np arange 0 10  y   np exp -x 3 0  f   interpolate interp1d x  y  fill value  extrapolate    print f 9  print f 11    and the output is   0 0497870683679 0 010394302658

User · Answer

I did it by adding a point to my initial arrays  In this way I avoid defining self-made functions  and the linear extrapolation  in the example below  right extrapolation  looks ok   import numpy as np   from scipy import interp as itp    xnew   np linspace 0 1 51    x1 xold -2    x2 xold -1    y1 yold -2    y2 yold -1    right val y1  xnew -1 -x1   y2-y1   x2-x1    x np append xold xnew -1     y np append yold right val    f   itp xnew x y

User · Answer

Here s an alternative method that uses only the numpy package  It takes advantage of numpy s array functions  so may be faster when interpolating extrapolating large arrays   import numpy as np  def extrap x  xp  yp          np interp function with linear extrapolation        y   np interp x  xp  yp      y   np where x lt xp 0   yp 0   x-xp 0    yp 0 -yp 1    xp 0 -xp 1    y      y   np where x gt xp -1   yp -1   x-xp -1    yp -1 -yp -2    xp -1 -xp -2    y      return y  x   np arange 0 10  y   np exp -x 3 0  xtest   np array  8 5 9 5    print np exp -xtest 3 0  print np interp xtest  x  y  print extrap xtest  x  y    Edit  Mark Mikofski s suggested modification of the  extrap  function   def extrap x  xp  yp          np interp function with linear extrapolation        y   np interp x  xp  yp      y x  lt  xp 0     yp 0     x x lt xp 0  -xp 0      yp 0 -yp 1      xp 0 -xp 1       y x  gt  xp -1    yp -1     x x gt xp -1  -xp -1    yp -1 -yp -2    xp -1 -xp -2       return y

User · Answer

I m afraid that there is no easy to do this in Scipy to my knowledge  You can  as I m fairly sure that you are aware  turn off the bounds errors and fill all function values beyond the range with a constant  but that doesn t really help  See this question on the mailing list for some more ideas  Maybe you could use some kind of piecewise function  but that seems like a major pain

User · Answer

I don t have enough reputation to comment  but in case somebody is looking for an extrapolation wrapper for a linear 2d-interpolation with scipy  I have adapted the answer that was given here for the 1d interpolation  def extrap2d interpolator   xs   interpolator x ys   interpolator y zs   interpolator z zs   np reshape zs   -1  len xs    def pointwise x  y       if x  lt  xs 0  or y  lt  ys 0           x1 index   np argmin np abs xs - x           x2 index   x1 index   1         y1 index   np argmin np abs ys - y           y2 index   y1 index   1         x1   xs x1 index          x2   xs x2 index          y1   ys y1 index          y2   ys y2 index          z11   zs x1 index  y1 index          z12   zs x1 index  y2 index          z21   zs x2 index  y1 index          z22   zs x2 index  y2 index           return  z11    x2 - x     y2 - y            z21    x - x1     y2 - y            z12    x2 - x     y - y1            z22    x - x1     y - y1               x2 - x1     y2 - y1    0 0        elif x  gt  xs -1  or y  gt  ys -1           x1 index   np argmin np abs xs - x           x2 index   x1 index - 1         y1 index   np argmin np abs ys - y           y2 index   y1 index - 1         x1   xs x1 index          x2   xs x2 index          y1   ys y1 index          y2   ys y2 index          z11   zs x1 index  y1 index          z12   zs x1 index  y2 index          z21   zs x2 index  y1 index          z22   zs x2 index  y2 index            return  z11    x2 - x     y2 - y            z21    x - x1     y2 - y            z12    x2 - x     y - y1            z22    x - x1     y - y1                x2 - x1     y2 - y1    0 0      else          return interpolator x  y  def ufunclike xs  ys       if  isinstance xs  int  or isinstance ys  int  or isinstance xs  np int32  or isinstance ys  np int32           res array   pointwise xs  ys      else          res array   np zeros  len xs   len ys            for x c in range len xs                res array x c       np array  pointwise xs x c   ys y c   for y c in range len ys     T      return res array return ufunclike  I haven t commented a lot and I am aware  that the code isn t super clean  If anybody sees any errors  please let me know  In my current use-case it is working without a problem

User · Answer

What about scipy interpolate splrep  with degree 1 and no smoothing     gt  gt  tck   scipy interpolate splrep  1  2  3  4  5    1  4  9  16  25   k 1  s 0   gt  gt  scipy interpolate splev 6  tck  34 0   It seems to do what you want  since 34   25    25 - 16

User · Answer

You can take a look at InterpolatedUnivariateSpline  Here an example using it    import matplotlib pyplot as plt import numpy as np from scipy interpolate import InterpolatedUnivariateSpline    given values xi   np array  0 2  0 5  0 7  0 9   yi   np array  0 3  -0 1  0 2  0 1     positions to inter extrapolate x   np linspace 0  1  50    spline order  1 linear  2 quadratic  3 cubic      order   1   do inter extrapolation s   InterpolatedUnivariateSpline xi  yi  k order  y   s x     example showing the interpolation for linear  quadratic and cubic interpolation plt figure   plt plot xi  yi  for order in range 1  4       s   InterpolatedUnivariateSpline xi  yi  k order      y   s x      plt plot x  y  plt show

User · Answer

1  Constant extrapolation  You can use interp function from scipy  it extrapolates left and right values as constant beyond the range    gt  gt  gt  from scipy import interp  arange  exp  gt  gt  gt  x   arange 0 10   gt  gt  gt  y   exp -x 3 0   gt  gt  gt  interp  9 10   x  y  array   0 04978707   0 04978707     2  Linear  or other custom  extrapolation  You can write a wrapper around an interpolation function which takes care of linear extrapolation  For example   from scipy interpolate import interp1d from scipy import arange  array  exp  def extrap1d interpolator       xs   interpolator x     ys   interpolator y      def pointwise x           if x  lt  xs 0               return ys 0   x-xs 0    ys 1 -ys 0    xs 1 -xs 0           elif x  gt  xs -1               return ys -1   x-xs -1    ys -1 -ys -2    xs -1 -xs -2           else              return interpolator x       def ufunclike xs           return array list map pointwise  array xs          return ufunclike   extrap1d takes an interpolation function and returns a function which can also extrapolate  And you can use it like this   x   arange 0 10  y   exp -x 3 0  f i   interp1d x  y  f x   extrap1d f i   print f x  9 10     Output     0 04978707  0 03009069

User · Answer

It may be faster to use boolean indexing with large datasets  since the algorithm checks if every point is in outside the interval  whereas boolean indexing allows an easier and faster comparison   For example     Necessary modules import numpy as np from scipy interpolate import interp1d    Original data x   np arange 0 10  y   np exp -x 3 0     Interpolator class f   interp1d x  y     Output range  quite large  xo   np arange 0  10  0 001     Boolean indexing approach    Generate an empty output array for  y  values yo   np empty like xo     Values lower than the minimum  x  are extrapolated at the same time low   xo  lt  f x 0  yo low     f y 0     xo low -f x 0    f y 1 -f y 0    f x 1 -f x 0      Values higher than the maximum  x  are extrapolated at same time high   xo  gt  f x -1  yo high    f y -1     xo high -f x -1    f y -1 -f y -2    f x -1 -f x -2      Values inside the interpolation range are interpolated directly inside   np logical and xo  gt   f x 0   xo  lt   f x -1   yo inside    f xo inside     In my case  with a data set of 300000 points  this means an speed up from 25 8 to 0 094 seconds  this is more than 250 times faster

[python] How to make scipy.interpolate give an extrapolated result beyond the input range?

Examples related to python

Examples related to math

Examples related to numpy

Examples related to scipy