solve_kcenter -- locate k facilities minimizing distance of most distant customer.
    Parameters:
        I - set of customers
        J - set of potential facilities
        c[i,j] - cost of servicing customer i from facility j
        k - number of facilities to be used
        delta - tolerance for terminating bisection
    Returns:
        - list of facilities to be used
        - edges linking them to customers

def solve_kcenter(I,J,c,k,delta):
    """solve_kcenter -- locate k facilities minimizing distance of most distant customer.
    Parameters:
        I - set of customers
        J - set of potential facilities
        c[i,j] - cost of servicing customer i from facility j
        k - number of facilities to be used
        delta - tolerance for terminating bisection
    Returns:
        - list of facilities to be used
        - edges linking them to customers
    """

    model = kcover(I,J,c,k)
    x,y,z = model.data

    facilities,edges = [],[]
    LB = 0
    UB = max(c[i,j] for (i,j) in c)
    model.setObjlimit(0.1)
    while UB-LB > delta:
        theta = (UB+LB) / 2.
        # print "\n\ncurrent theta:", theta
        for j in J:
            for i in I:
                if c[i,j]>theta:
                    model.chgVarUb(x[i,j], 0.0)
                else:
                    model.chgVarUb(x[i,j], 1.0)

        # model.Params.OutputFlag = 0 # silent mode
        model.setObjlimit(.1)

        model.optimize()

        if model.getStatus == "optimal":
            # infeasibility = sum([z[i].X for i in I])
            # print "infeasibility=",infeasibility
            UB = theta
            facilities = [j for j in y if model.getVal(y[j]) > .5]
            edges = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5]
            # print "updated solution:"
            # print "facilities",facilities
            # print "edges",edges
        else:   # infeasibility > 0:
            LB = theta

    return facilities,edges