SA with scipy

even_points.py

""" Placing points in a polygon evenly with simulated annealing.
"""

import math
import random

import descartes
from matplotlib import pyplot
import scipy.optimize
import shapely.geometry


blue = '#6699cc'
red = '#de0000'
points = 16
vertices = [
    [0, 0],
    [90, 10],
    [120, 30],
    [110, 60],
    [80, 90],
    [40, 110],
    [10, 50],
    [15, 20]
]
fixed_initial_points = [
    (20, 10),
    (20, 50),
    (40, 80),
    (60, 60),
    (70, 80),
    (70, 70),
    (70, 80),
    (100, 40),
    (110, 35),
    (70, 20)
]


def random_point_in(polygon):
  """Finds a random (x, y) point in a shapely Polygon instance."""
  while True:
    x = random.uniform(polygon.bounds[0], polygon.bounds[2])
    y = random.uniform(polygon.bounds[1], polygon.bounds[3])
    if polygon.contains(shapely.geometry.Point(x, y)):
      return (x, y)


def distance(a, b):
  """Finds distance between two (x, y) points."""
  return math.sqrt((a[0] - b[0])**2 + (a[1] - b[1])**2)


def fitness(points, *params):
  polygon = params[0]
  number_of_neighbors = int(len(points) / 3) + 1
  multiplier = 1
  # Penalize points that are outside the polygon.
  for coordinates in points:
    if not polygon.contains(shapely.geometry.Point(coordinates)):
      multiplier *= 50
  # Penalize points that do not fill the polygon.
  min_x = min([p[0] for p in points])
  max_x = max([p[0] for p in points])
  min_y = min([p[1] for p in points])
  max_y = max([p[1] for p in points])
  bounding_box_area = (max_x - min_x) * (max_y - min_y)
  multiplier *= (polygon.area / bounding_box_area)
  # MMSD: Test random points in the polygon.
  test_points = [random_point_in(polygon) for _ in range(10)]
  distances = []
  for test in test_points:
    neighbors = sorted(points, key=lambda p2: distance(test, p2))
    nearest_neighbors = neighbors[1:number_of_neighbors+1]
    d = [distance(test, p2) for p2 in nearest_neighbors]
    distances.append(sum(d) / len(d))
  average = sum(distances) / len(distances)
  # Variance check.
  nearest_distances = []
  for test in test_points:
    nearest = sorted(points, key=lambda p2: distance(test, p2))[0]
    nearest_distances.append(distance(test, nearest))
  nearest_average = sum(nearest_distances) / len(nearest_distances)
  variance = (sum([(p - nearest_average)**2 for p in nearest_distances]) /
              (len(nearest_distances) - 1))
  return multiplier * average * variance


if __name__ == '__main__':
  # Setup shapes.
  polygon = shapely.geometry.Polygon(vertices)
  patch = descartes.PolygonPatch(polygon, fc=blue, ec=blue, alpha=0.5)
  initial_points = [random_point_in(polygon) for _ in range(points)]
  #initial_points = fixed_initial_points
  initial_temp = fitness(initial_points, polygon)

  # Anneal.
  result = scipy.optimize.anneal(fitness, initial_points, args=(polygon,),
                                 learn_rate=0.2, schedule='boltzmann',
                                 T0=initial_temp, full_output=True, disp=True,
                                 )
  print result
  points = [list(p) for p in result[0]]

  # Setup the plot.
  figure = pyplot.figure(dpi=90)
  figure.set_size_inches(9, 9)
  axes = figure.add_subplot(1, 1, 1)

  # Plot polygon and points.
  axes.add_patch(patch)
  for point in initial_points:
    axes.plot(point[0], point[1], 'b.')
  for point in points:
    axes.plot(point[0], point[1], 'g*')

  # Label, limit and save.
  axes.set_xlim(polygon.bounds[0], polygon.bounds[2])
  axes.set_ylim(polygon.bounds[1], polygon.bounds[3])
  axes.set_title('blue dots -> green stars via SA')
  pyplot.axis('equal')
  figure.savefig('out.png')

out.png