# HG changeset patch
# User Jeff Hammel <k0scist@gmail.com>
# Date 1409844070 25200
# Node ID 097296d6132e09ec3fa4c2c9ce2a71b1ef60a5cf
# Parent  690778ffd3029736e6d88f0046baf4da5d92a0b9
add interpolation

diff -r 690778ffd302 -r 097296d6132e numerics/interpolation.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/numerics/interpolation.py	Thu Sep 04 08:21:10 2014 -0700
@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+
+
+"""
+interpolation
+"""
+
+__all__ = ['neighbors', 'linear_interpolation']
+
+def neighbors(start, finish):
+    """
+    returns the neighbors in finish from start
+    assumes both are sorted
+    """
+    assert finish
+    index = 0
+    retval = []
+    for x in start:
+        # endpoints
+        if x < finish[0]:
+            retval.append((None, 0))
+            continue
+        if x > finish[-1]:
+            retval.append((len(finish)-1, None))
+            continue
+        # traverse
+        try:
+            while True:
+                if x < finish[index] or x > finish[index+1]:
+                    index += 1
+                    continue
+                else:
+                    break
+            retval.append((index, index+1))
+        except IndexError:
+            retval.append((len(finish)-2, len(finish)-1))
+
+    return retval
+
+
+def linear_interpolation(data, points):
+    """
+    linearly interpolate data to points
+
+    data -- iterable of 2-tuples (or equivalent) of `x,y`
+    points -- `x`-values to interpolate to
+    """
+
+    # ensure we are sorted
+    data = sorted(data, key=lambda x: x[0])
+    points = sorted(points)
+
+    # get the neighbors
+    x = [value[0] for value in data]
+    nearest_neighbors = neighbors(points, x)
+
+    # we don't support endpoints yet; this is interpolation, not extrapolation
+    assert not any([(neighbor[0] is None or neighbor[1] is None)
+                    for neighbor in nearest_neighbors])
+
+    retval = []
+    for index, (left, right) in enumerate(nearest_neighbors):
+        # linearly interpolate
+        ratio = (points[index] - data[left][0])/float(data[right][0] - data[left][0])
+        retval.append(ratio*data[right][1] + (1.-ratio)*data[left][1])
+    return retval
+
+