geometry.py

import sys
import random
import pylab
import math
import flow
import viz
import collections

N = 200
D = 20
RANGE = 20
random.seed(1)

points = [ (random.uniform(-RANGE, RANGE), random.uniform(-RANGE, RANGE)) for _ in xrange(N) ]

viz.figure()
viz.plot_points(points)
viz.show(RANGE)

def max_bipartite_independent(top, bot):

    graph = collections.defaultdict(dict)

    top_set = set( [i for i in xrange( len(top)) ])
    bot_set = set( [i for i in xrange( len(bot)) ])

    for i in xrange(len(top)):
        for j in xrange(len(bot)):

            xi, yi = top[i]
            xj, yj = bot[j]

            if distance(xi, yi, xj, yj) > F:

                node_i = 'TOP' + str(i)
                node_j = 'BOT' + str(j)

                src = 'SOURCE-NODE'
                snk = 'SINK-NODE'

                graph[node_i][node_j] = 1
                graph[src][node_i] =  sys.maxint
                graph[node_j][snk] =  sys.maxint

                top_set.discard(i)
                bot_set.discard(j)

    f = flow.max_flow(graph, 'SOURCE-NODE', 'SINK-NODE')
    solution = f + len(top_set) + len(bot_set) + 2
    return solution

def area(a, b, c):
     return (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])

def distance(xi, yi, xj, yj):
    return math.hypot(xi - xj, yi - yj)

best = 0

for i in xrange(N):
    for j in xrange(i + 1, N):
        xi, yi = points[i]
        xj, yj = points[j]

        if distance(xi, yi, xj, yj) > D: continue

        # Furthest pair in our thing
        F = distance(xi, yi, xj, yj)

        viz.figure()
        viz.plot_points(points)
        viz.plot_line(xi, yi, xj, yj)
        viz.plot_circle(xi, yi, D)
        viz.plot_circle(xj, yj, D)
        viz.show(RANGE)

        lens = []

        for k in xrange(N):

            if k == i: continue
            if k == j: continue

            xk, yk = points[k]

            if distance(xi, yi, xk, yk) > F: continue
            if distance(xj, yj, xk, yk) > F: continue

            lens.append(points[k])


        viz.figure()
        if len(lens) > 0: viz.plot_points(lens)
        viz.plot_line(xi, yi, xj, yj)
        viz.plot_circle(xi, yi, F)
        viz.plot_circle(xj, yj, F)
        viz.show(RANGE)


        top, bot = [], []
        M = len(lens)

        for k in xrange(M):
            xk, yk = lens[k]

            if area((xi, yi), (xj, yj), (xk, yk)) >= 0:
                top.append(lens[k])
            else:
                bot.append(lens[k])

        viz.figure()

        for t in top:
            for b in bot:
                if distance(t[0], t[1], b[0], b[1]) > D:
                    viz.plot_line(t[0], t[1], b[0], b[1], color='k')

        if len(top) > 0: viz.plot_points(top, color='bo')
        if len(bot) > 0: viz.plot_points(bot, color='ro')

        viz.plot_line(xi, yi, xj, yj)
        viz.plot_circle(xi, yi, F)
        viz.plot_circle(xj, yj, F)
        viz.show(RANGE)

        solution = max_bipartite_independent(top, bot)
        best = max(best, solution)
        
print 'MAXIMUM: ', best