import csv import math import os import re from collections import deque class Vector(complex): def __getitem__(self, index): if index == 0: return self.real elif index == 1: return self.imag else: raise IndexError def load(directory): def _read(iterable): for x, y in iterable: yield Vector(float(x), float(y)) def _load(filename): with open(filename) as fd: reader = csv.reader(fd) return list(_read(reader)) def _files(directory): for file in os.listdir(directory): match = re.match(r"SWI_(-?\d+)\.csv", file) if match: yield int(match.group(1)), os.path.join(directory, file) return [(x, _load(y)) for x, y in sorted(_files(directory), key=lambda x: x[0])] def look_downwards(data, x, start): for i in range(start, 0, -1): if data[i - 1][0] < x: break else: raise IndexError return i - 1 def look_upwards(data, x, start): for i in range(start, len(data)): if data[i + 1][0] > x: break else: raise IndexError return i def find_segment(data, x): width = data[-1][0] - data[0][0] relative = x - data[0][0] candidate = math.floor(relative / width * len(data)) look = look_downwards if data[candidate][0] > x else look_upwards # May raise IndexError candidate = look(data, x, candidate) return candidate, candidate + 1 def find_boundary_curves(swis, x, y): segments = deque() for index, data in swis: i, j = find_segment(data, x) if data[i][1] > y and data[j][1] > y: segments.append((index, data, i, j)) break if data[i][1] < y and data[j][1] < y: if segments: segments.popleft() segments.append((index, data, i, j)) if len(segments) == 3: middle = segments[1][1] run = middle[j][0] - middle[i][0] if run == 0: raise RuntimeError # tidy up dataset slope = (middle[j][1] - middle[i][1]) / run intercept = middle[j][1] - slope * middle[j][0] value = slope * x + intercept if value == y: raise RuntimeError # Exactly on point; SWI == index if value < y: segments.popleft() else: segments.pop() if len(segments) == 1: raise RuntimeError # SWI == -10 return segments def calculate_swi(segments, x, y): vec = Vector(x, y) low = segments[0] high = segments[1] dist_to_low = min(abs(vec - p) for p in (low[1][low[2]], low[1][low[2]])) dist_to_high = min(abs(vec - p) for p in (high[1][high[2]], high[1][high[2]])) return dist_to_low / (dist_to_low + dist_to_high) * (high[0] - low[0]) + low[0]