From ff0dd2f8c5e9ed5bd63fcc47e09d46a7b2812ff1 Mon Sep 17 00:00:00 2001 From: Aki Date: Tue, 30 Aug 2022 00:15:26 +0200 Subject: Wrapped script in a setup.py package --- .gitignore | 2 + plot.py | 131 --------------------------------------------------- setup.py | 8 ++++ szilagyi/nomogram.py | 131 +++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 141 insertions(+), 131 deletions(-) delete mode 100644 plot.py create mode 100755 setup.py create mode 100644 szilagyi/nomogram.py diff --git a/.gitignore b/.gitignore index 670a936..a6d7b32 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ __pycache__/ .venv/ +*.egg-info/ +build*/ diff --git a/plot.py b/plot.py deleted file mode 100644 index eb2665d..0000000 --- a/plot.py +++ /dev/null @@ -1,131 +0,0 @@ -import csv -import math -import os -import re -from collections import deque - -import matplotlib.pyplot as plot - - -def load(directory): - def _read(iterable): - for x, y in iterable: - yield 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 dist(x1, y1, x2, y2): - return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) - - -def calculate_swi(segments, x, y): - low = segments[0] - high = segments[1] - dist_to_low = min(dist(p[0], p[1], x, y) for p in (low[1][low[2]], low[1][low[2]])) - dist_to_high = min(dist(p[0], p[1], x, y) 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] - - -swis = load("dataset") - - -def _swi(x, y): # dt, depth - segments = find_boundary_curves(swis, x, y) - return calculate_swi(segments, x, y) - - -def grid(start, end, steps): - step = (end - start) / steps - i = start - while i <= end: - yield i - i += step - - -C = [] -X = list(grid(0, 40, 1000)) -Y = list(grid(0, 50000, 1000)) -for y_scaled in range(0, 1000): - y = y_scaled * 50 - row = [] - for x_scaled in range(0, 1000): - x = x_scaled / 25 - try: - swi = _swi(x, y) - except (IndexError, RuntimeError): - swi = -10 - row.append(swi) - C.append(row) -plot.pcolormesh(X, Y, C, cmap='viridis', vmin=-10, vmax=10, rasterized=True) -for _, data in swis: - plot.plot([x[0] for x in data], [x[1] for x in data]) -plot.show() diff --git a/setup.py b/setup.py new file mode 100755 index 0000000..d488731 --- /dev/null +++ b/setup.py @@ -0,0 +1,8 @@ +#!/usr/bin/env python + +from setuptools import setup + +setup( + name="szilagyi", + packages=["szilagyi"], +) diff --git a/szilagyi/nomogram.py b/szilagyi/nomogram.py new file mode 100644 index 0000000..eb2665d --- /dev/null +++ b/szilagyi/nomogram.py @@ -0,0 +1,131 @@ +import csv +import math +import os +import re +from collections import deque + +import matplotlib.pyplot as plot + + +def load(directory): + def _read(iterable): + for x, y in iterable: + yield 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 dist(x1, y1, x2, y2): + return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) + + +def calculate_swi(segments, x, y): + low = segments[0] + high = segments[1] + dist_to_low = min(dist(p[0], p[1], x, y) for p in (low[1][low[2]], low[1][low[2]])) + dist_to_high = min(dist(p[0], p[1], x, y) 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] + + +swis = load("dataset") + + +def _swi(x, y): # dt, depth + segments = find_boundary_curves(swis, x, y) + return calculate_swi(segments, x, y) + + +def grid(start, end, steps): + step = (end - start) / steps + i = start + while i <= end: + yield i + i += step + + +C = [] +X = list(grid(0, 40, 1000)) +Y = list(grid(0, 50000, 1000)) +for y_scaled in range(0, 1000): + y = y_scaled * 50 + row = [] + for x_scaled in range(0, 1000): + x = x_scaled / 25 + try: + swi = _swi(x, y) + except (IndexError, RuntimeError): + swi = -10 + row.append(swi) + C.append(row) +plot.pcolormesh(X, Y, C, cmap='viridis', vmin=-10, vmax=10, rasterized=True) +for _, data in swis: + plot.plot([x[0] for x in data], [x[1] for x in data]) +plot.show() -- cgit v1.1