|
4
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3
|
|
|
4 """
|
|
|
5 interpolation
|
|
|
6 """
|
|
|
7
|
|
36
|
8 import argparse
|
|
|
9
|
|
10
|
10 __all__ = ['neighbors', 'linear_interpolation', 'InterpolateParser']
|
|
4
|
11
|
|
|
12 def neighbors(start, finish):
|
|
|
13 """
|
|
|
14 returns the neighbors in finish from start
|
|
|
15 assumes both are sorted
|
|
|
16 """
|
|
|
17 assert finish
|
|
|
18 index = 0
|
|
|
19 retval = []
|
|
|
20 for x in start:
|
|
|
21 # endpoints
|
|
|
22 if x < finish[0]:
|
|
|
23 retval.append((None, 0))
|
|
|
24 continue
|
|
|
25 if x > finish[-1]:
|
|
|
26 retval.append((len(finish)-1, None))
|
|
|
27 continue
|
|
|
28 # traverse
|
|
|
29 try:
|
|
|
30 while True:
|
|
|
31 if x < finish[index] or x > finish[index+1]:
|
|
|
32 index += 1
|
|
|
33 continue
|
|
|
34 else:
|
|
|
35 break
|
|
|
36 retval.append((index, index+1))
|
|
|
37 except IndexError:
|
|
|
38 retval.append((len(finish)-2, len(finish)-1))
|
|
|
39
|
|
|
40 return retval
|
|
|
41
|
|
|
42
|
|
|
43 def linear_interpolation(data, points):
|
|
|
44 """
|
|
|
45 linearly interpolate data to points
|
|
|
46
|
|
|
47 data -- iterable of 2-tuples (or equivalent) of `x,y`
|
|
|
48 points -- `x`-values to interpolate to
|
|
|
49 """
|
|
|
50
|
|
|
51 # ensure we are sorted
|
|
|
52 data = sorted(data, key=lambda x: x[0])
|
|
|
53 points = sorted(points)
|
|
|
54
|
|
|
55 # get the neighbors
|
|
|
56 x = [value[0] for value in data]
|
|
|
57 nearest_neighbors = neighbors(points, x)
|
|
|
58
|
|
|
59 # we don't support endpoints yet; this is interpolation, not extrapolation
|
|
|
60 assert not any([(neighbor[0] is None or neighbor[1] is None)
|
|
|
61 for neighbor in nearest_neighbors])
|
|
|
62
|
|
|
63 retval = []
|
|
|
64 for index, (left, right) in enumerate(nearest_neighbors):
|
|
|
65 # linearly interpolate
|
|
|
66 ratio = (points[index] - data[left][0])/float(data[right][0] - data[left][0])
|
|
|
67 retval.append(ratio*data[right][1] + (1.-ratio)*data[left][1])
|
|
|
68 return retval
|
|
|
69
|
|
10
|
70 class InterpolateParser(argparse.ArgumentParser):
|
|
|
71 """CLI option parser"""
|
|
4
|
72
|
|
10
|
73 def __init__(self, **kwargs):
|
|
|
74 kwargs.setdefault('description', __doc__)
|
|
|
75 argparse.ArgumentParser.__init__(self, **kwargs)
|
|
|
76 self.add_argument('input', nargs='?',
|
|
|
77 type=argparse.FileType('r'), default=sys.stdin,
|
|
|
78 help='input file, or read from stdin if ommitted')
|
|
|
79 self.add_argument('-o', '--output', dest='output',
|
|
|
80 type=argparse.FileType('w'), default=sys.stdout,
|
|
|
81 help="output file, or stdout if ommitted")
|
|
|
82 self.add_argument('--points', '--print-points', dest='print_points',
|
|
|
83 action='store_true', default=False,
|
|
|
84 help="print the points to interpolate to and exit")
|
|
|
85 self.options = None
|
|
|
86
|
|
|
87 def parse_args(self, *args, **kw):
|
|
|
88 options = argparse.ArgumentParser.parse_args(self, *args, **kw)
|
|
|
89 self.validate(options)
|
|
|
90 self.options = options
|
|
|
91 return options
|
|
|
92
|
|
|
93 def validate(self, options):
|
|
|
94 """validate options"""
|
|
|
95
|
|
|
96 def main(args=sys.argv[1:]):
|
|
|
97 """CLI"""
|
|
|
98
|
|
|
99 # parse command line options
|
|
|
100 parser = InterpolateParser()
|
|
|
101 options = parser.parse_args(args)
|
|
|
102
|
|
|
103 # read the CSV
|
|
|
104 reader = csv.reader(options.input)
|
|
|
105 data = [[float(col) for col in row] for row in reader]
|
|
|
106 ncols = set([len(row) for row in data])
|
|
|
107 assert len(ncols) == 1
|
|
|
108 ncols = ncols.pop()
|
|
|
109 assert ncols > 1
|
|
|
110
|
|
|
111 # get `x` values
|
|
|
112 data = sorted(data, key=lambda x: x[0])
|
|
|
113 x = [row[0] for row in data]
|
|
|
114 xmin = int(x[0]) + 1
|
|
|
115 xmax = int(x[-1])
|
|
|
116 points = range(xmin, xmax+1)
|
|
|
117 if options.print_points:
|
|
|
118 print ('\n'.join([str(point) for point in points]))
|
|
|
119 return
|
|
|
120
|
|
|
121 # make into x,y series
|
|
|
122 series = [[(row[0], row[col]) for row in data]
|
|
|
123 for col in range(1,ncols)]
|
|
|
124
|
|
|
125 # interpolate
|
|
|
126 interpolated = [linear_interpolation(s, points) for s in series]
|
|
|
127
|
|
|
128 # output interpolated data
|
|
|
129 writer = csv.writer(options.output)
|
|
|
130 for row in zip(points, *interpolated):
|
|
|
131 writer.writerow(row)
|
|
|
132
|
|
|
133 if __name__ == '__main__':
|
|
|
134 main()
|
