"""Finding the minimum enclosing circle using different optimization
algorithms and different functions
"""

from math import sqrt
import math 
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

solvers = [
    'nelder-mead',
    'powell',
    'cg',
    'bfgs',
#    'Newton-CG',
    'L-BFGS-B',
    'TNC',
    'COBYLA',
    'SLSQP',
#    'dogleg',
#    'trust-ncg',
#    'trust-exact',
#    'trust-krylov'
]

######################################################################

points_X = [1.0, 3.0, 2.5, 4.0, 5.0, 6.0, 5.0]
points_Y = [3.0, 1.0, 3.0, 6.0, 7.0, 7.0, 2.0]

points = list(zip(points_X, points_Y))

######################################################################
# Minimum enclosing circle solver


def distance(p, q):
    """ compute distance between d-dimensional points p and q """
    return sqrt(sum([(i-j)**2 for i, j in zip(p, q)]))

def distance_max(c):
    """computes the distance to the point furthest away from c"""
    return max([distance(p, c) for p in points])

def distance_max_square(c):
    """computes the distance to the point furthest away from c"""
    return max([(p[0]-c[0])**2+(p[1]-c[1])**2 for p in points])

def distance_smooth_max(c):
    """computes the distance to the point furthest away from c"""
    return math.log(sum([math.exp((p[0]-c[0])**2+(p[1]-c[1])**2) for p in points]))

def solve_and_plot(function, method, color, linestyle):
    solution = minimize(function, [0.0, 0.0], method=method)
    cx, cy = solution.x
    c = (cx, cy)
    r = distance_max(c)
    print(method, r)

    plt.plot([cx], [cy], color+"o")     # center point
    circle = plt.Circle(c, r, color=color, fill=False, linestyle=linestyle)
    ax.add_artist(circle)

plt.suptitle('Smallest enclosing circle\n'
             'Function: red = max distance, '
             'green = max distance-squared, '
             'blue = "smooth"-max distance-squared')
    
for idx, opt_algorithm in enumerate(solvers):
    plt.subplot(2, 4, idx+1)
    plt.title(opt_algorithm)
    ax = plt.gca()
    ax.set_xlim((-1,9))
    ax.set_ylim((-1,9))
    ax.set_aspect("equal")
    
    plt.plot(points_X, points_Y, "g.")   # data points (green dots)
    
    solve_and_plot(distance_max, opt_algorithm, 'r', "-")
    solve_and_plot(distance_max_square, opt_algorithm, 'g', "--")
    solve_and_plot(distance_smooth_max, opt_algorithm, 'b', ":")

plt.show()