"""
center.py:  Using Gurobi for solving the k-center problem:
    bisection on the minimum max distance.

Problem:
minimize the maximum travel distance from n customers to k facilities.
            
Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2010
"""
from gurobipy import *

def kcenter(m, n, c, k, max_c):
    model = Model("k-center")

    z, y, x = {}, {}, {}
    for i in range(n):
        z[i] = model.addVar(obj=1, vtype="B", name="z[%s]"%i)
    for j in range(m):
        y[j] = model.addVar(obj=0, vtype="B", name="y[%s]"%j)
        for i in range(n):
            x[i,j] = model.addVar(obj=0, vtype="B", name="x[%s,%s]"%(i,j))
    model.update()

    for i in range(n):
        coef = [1 for j in range(m+1)]
        var = [x[i,j] for j in range(m)]
        var.append(z[i])
        model.addConstr(LinExpr(coef,var), "=", 1, name="Assign[%s]"%i)
    for j in range(m):
        for i in range(n):
            model.addConstr(x[i,j], "<", y[j], name="Strong[%s,%s]"%(i,j))
    coef = [1 for j in range(m)]
    var = [y[j] for j in range(m)]
    model.addConstr(LinExpr(coef,var), "=", rhs=k, name="k_center")              

    model.update()
    model.__data = x,y,z
    return model


def solve_kcenter(m, n, c, k, max_c, delta):
    model = kcenter(m, n, c, k, max_c)
    x,y,z = model.__data

    nodes = []
    edges = []
    
    LB = 0
    UB = max_c
    while UB-LB > delta:        
        theta = (UB+LB) / 2.
        # print "\n\ncurrent theta:", theta
        for j in range(m):
            for i in range(n):
                if c[i,j]>theta:
                    x[i,j].UB = 0
                else:
                    x[i,j].UB = 1.0
        model.update()
        model.Params.OutputFlag = 0 # silent mode
        model.optimize()

        infeasibility = sum([z[i].X for i in range(m)])
        # print "infeasibility=",infeasibility
        if infeasibility > 0:
            LB = theta
        else:
            UB = theta
            nodes = [j for j in y if y[j].X > .5]
            edges = [(i,j) for (i,j) in x if x[i,j].X > .5]
            # print "updated solution:"
            # print "nodes", nodes
            # print "edges", edges

    return nodes, edges





import math
import random
def distance(x1, y1, x2, y2):
    return math.sqrt((x2-x1)**2 + (y2-y1)**2)


def make_data(n):
    # positions of the points in the plane
    x = [random.random() for i in range(n)]
    y = [random.random() for i in range(n)]

    c = {}
    max_c = 0.
    for i in range(n):
        for j in range(n):
            c[i,j] = distance(x[i],y[i],x[j],y[j])
            max_c = max(c[i,j],max_c)

    return c, x, y, max_c

                
if __name__ == "__main__":
    random.seed(67)
    n = 200
    c, x_pos, y_pos, max_c = make_data(n)
    m = n
    k = 20
    delta = 1.e-4
    nodes, edges = solve_kcenter(m, n, c, k, max_c, delta)
    print "Selected nodes:", nodes
    print "Edges:", edges
    print [((i,j),c[i,j]) for (i,j) in edges]
    print "max c:", max([c[i,j] for (i,j) in edges])

    try: # plot the result using networkx and matplotlib
        import networkx as NX
        import matplotlib.pyplot as P 
        P.ion() # interactive mode on
        G = NX.Graph()

        other = [j for j in range(m)  if j not in nodes]
        G.add_nodes_from(nodes)
        G.add_nodes_from(other)
        for (i,j) in edges:
            G.add_edge(i,j)

        position = {}
        for i in range(n):
            position[i]=(x_pos[i],y_pos[i])

        NX.draw(G, position, node_color='y', nodelist=nodes)
        NX.draw(G, position, node_color='g', nodelist=other)
        raw_input("press [enter] to continue")
    except ImportError:
        print "install 'networkx' and 'matplotlib' for plotting"