How to implement linear interpolation?

import scipy.interpolate
y_interp = scipy.interpolate.interp1d(x, y)
print y_interp(5.0)

scipy.interpolate.interp1d does linear interpolation by and can be customized to handle error conditions.

I thought up a rather elegant solution (IMHO), so I can't resist posting it:

from bisect import bisect_left

class Interpolate(object):
    def __init__(self, x_list, y_list):
        if any(y - x <= 0 for x, y in zip(x_list, x_list[1:])):
            raise ValueError("x_list must be in strictly ascending order!")
        x_list = self.x_list = map(float, x_list)
        y_list = self.y_list = map(float, y_list)
        intervals = zip(x_list, x_list[1:], y_list, y_list[1:])
        self.slopes = [(y2 - y1)/(x2 - x1) for x1, x2, y1, y2 in intervals]

    def __getitem__(self, x):
        i = bisect_left(self.x_list, x) - 1
        return self.y_list[i] + self.slopes[i] * (x - self.x_list[i])

I map to float so that integer division (python <= 2.7) won't kick in and ruin things if x1, x2, y1 and y2 are all integers for some iterval.

In __getitem__ I'm taking advantage of the fact that self.x_list is sorted in ascending order by using bisect_left to (very) quickly find the index of the largest element smaller than x in self.x_list.

Use the class like this:

i = Interpolate([1, 2.5, 3.4, 5.8, 6], [2, 4, 5.8, 4.3, 4])
# Get the interpolated value at x = 4:
y = i[4]

I've not dealt with the border conditions at all here, for simplicity. As it is, i[x] for x < 1 will work as if the line from (2.5, 4) to (1, 2) had been extended to minus infinity, while i[x] for x == 1 or x > 6 will raise an IndexError. Better would be to raise an IndexError in all cases, but this is left as an exercise for the reader. :)

As I understand your question, you want to write some function y = interpolate(x_values, y_values, x), which will give you the y value at some x? The basic idea then follows these steps:

1. Find the indices of the values in x_values which define an interval containing x. For instance, for x=3 with your example lists, the containing interval would be [x1,x2]=[2.5,3.4], and the indices would be i1=1, i2=2
2. Calculate the slope on this interval by (y_values[i2]-y_values[i1])/(x_values[i2]-x_values[i1]) (ie dy/dx).
3. The value at x is now the value at x1 plus the slope multiplied by the distance from x1.

You will additionally need to decide what happens if x is outside the interval of x_values, either it's an error, or you could interpolate "backwards", assuming the slope is the same as the first/last interval.

Did this help, or did you need more specific advice?