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