dish-height.py

#!/usr/bin/env python3

import numpy as np
import sys
import matplotlib.pylab as plt

# p(x) = sum([ p * x**i for i, p in enumerate(ps)])
def differentiate(ps):
    return [ p * i for i, p in enumerate(ps) ][1:]

def roots(ps):
    return np.roots(ps[::-1])

def evaluate(ps, x):
    return sum([ p * x**i for i, p in enumerate(ps)])

def distance(dish_a, dish_b):
    radius_0, shape_0 = dish_a
    radius_1, shape_1 = dish_b
    radius = min(radius_0, radius_1)
    diff = shape_0 - shape_1
    dd = differentiate(diff)
    p_roots = roots(dd)

    def is_maximum(root):
        return evaluate(dd, root) > 0.0

    def is_in_range(root):
        return 0.0 < root <= radius

    def is_real(root):
        return abs(root.imag) <= 0.00000001

    def root_filter(root):
        return is_maximum(root) and is_in_range(root) and is_real(root)

    #  print("roots: ", p_roots)
    p_roots = [ root for root in p_roots if root_filter(root) ]
    distances = [ evaluate(diff, root) for root in p_roots ]
    distances += [ evaluate(diff, 0.0), evaluate(diff, radius) ]
    min_distances = max(distances)
    return min_distances

def read_dish(line):
    radius, ps = line.split(":")
    radius = float(radius)
    even_ps = np.array([float(c) for c in ps.split(" ")])
    ps = np.zeros(2*len(even_ps)+1)
    ps[2::2] = even_ps
    return radius, ps

def pad(ps, n):
    padded_ps = np.zeros(n)
    padded_ps[:len(ps)] = ps
    return padded_ps

def plot_dishes(dishes):
    assert len(dishes) > 0
    for radius, shape in dishes: 
        n = 100
        xs = np.linspace(-radius, radius, n)
        ys = [ evaluate(shape, x) for x in xs ]
        plt.plot(xs, ys)
    plt.legend(range(len(dishes)))
    plt.show()

if __name__ == "__main__":
    with open(sys.argv[1]) as f:
        dishes = [ read_dish(line) for line in f.readlines() if not line.startswith("#")  ]

    n = max([len(shape) for r, shape in dishes])
    dishes = [ (r, pad(shape, n)) for r, shape in dishes ]

    y_offsets = [ 0.0 ]
    for i in range(1, len(dishes)):
        distances = [ distance(dishes[j], dishes[i]) for j in range(i) ]
        y_offsets.append(max(distances))
        dishes[i][1][0] += y_offsets[i]

    stack_height = y_offsets[-1]
    print(stack_height)
    plot_dishes(dishes)