from math import sqrt
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

points = [(1.0, 3.0), (3.0, 1.0), (2.5, 3.0), (4.0, 6.0),
          (5.0, 7.0), (6.0, 7.0), (5.0, 2.0)]

# Minimum enclosing circle solver

trace = []

def distance(p, q):
    return sqrt((p[0]-q[0])**2 + (p[1]-q[1])**2)
    
def distance_max(q):
    '''computes the distance to the point furthest away from c'''
    dist = max(distance(p, q) for p in points)
    trace.append((*q, dist))
    return dist

solution = minimize(distance_max, [0.0, 0.0], method='nelder-mead')

center = solution.x
radius = solution.fun

#  Minimum enclosing circle - 3D surface plot

points_x, points_y = zip(*points)
trace_x, trace_y, trace_z = zip(*trace)

xs = points_x + trace_x
ys = points_y + trace_y

x_min = min(xs)
x_max = max(xs)
y_min = min(ys)
y_max = max(ys)

# enforce apsect ratio
x_max = max(x_max, x_min + y_max - y_min)
y_max = max(y_max, y_min + x_max - x_min)

grid_x = np.linspace(x_min, x_max, 100)
grid_y = np.linspace(y_min, y_max, 100)

X, Y = np.meshgrid(grid_x, grid_y)

Z = np.zeros(X.shape)
for px, py in points:
    Z = np.maximum(Z, (X - px)**2 + (Y - py)**2)
Z = np.sqrt(Z)

# plot_surface requires X, Y, Z are 2D numpy.arrays 
ax = plt.subplot(1, 2, 1, projection='3d')
ax.plot_surface(X, Y, Z, cmap='plasma', alpha=0.7)
ax.plot(trace_x, trace_y, trace_z, '.-', c='darkblue')
ax.scatter(*center, radius, 'o', c='red')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('max distance')
ax.set_title('plot_surface')

# Minimum enclosing circle - contour plot

plt.subplot(1, 2, 2)
plt.title('pyplot.contour')
plt.plot(trace_x, trace_y, '.-', color='darkblue')
plt.plot(points_x, points_y, 'o', color='darkgreen')
plt.plot(*center, 'o', c='red')
qcs = plt.contour(X, Y, Z, levels=30, cmap='plasma')
plt.clabel(qcs, inline=1, fontsize=8, fmt='%.1f')

plt.suptitle('Maximum distance to an input point')
plt.tight_layout()
plt.show()
