solve_tsp -- solve the traveling salesman problem
       - start with assignment model
       - check flow from a source to every other node;
          - if no flow, a sub-cycle has been found --> add cut
          - otherwise, the solution is optimal
    Parameters:
        - V: set/list of nodes in the graph
        - c[i,j]: cost for traversing edge (i,j)
    Returns the optimum objective value and the list of edges used.

def solve_tsp(V,c):
    """solve_tsp -- solve the traveling salesman problem
       - start with assignment model
       - check flow from a source to every other node;
          - if no flow, a sub-cycle has been found --> add cut
          - otherwise, the solution is optimal
    Parameters:
        - V: set/list of nodes in the graph
        - c[i,j]: cost for traversing edge (i,j)
    Returns the optimum objective value and the list of edges used.
    """

    def addcut(X):
        for sink in V[1:]:
            mflow = maxflow(V,X,V[0],sink)
            mflow.optimize()
            f,cons = mflow.data
            if mflow.ObjVal < 2-EPS:  # no flow to sink, can add cut
                break
        else:
            return False

        #add a cut/constraint
        CutA = set([V[0]])
        for i in cons:
            if cons[i].Pi <= -1+EPS:
                CutA.add(i)
        CutB = set(V) - CutA
        main.addCons(
            quicksum(x[i,j] for i in CutA for j in CutB if j>i) + \
            quicksum(x[j,i] for i in CutA for j in CutB if j<i) >= 2)
        print("mflow:",mflow.getObjVal(),"cut:",CutA,"+",CutB,">= 2")
        print("mflow:",mflow.getObjVal(),"cut:",[(i,j) for i in CutA for j in CutB if j>i],"+",[(j,i) for i in CutA for j in CutB if j<i],">= 2")
        return True

    def isMIP(x):
        for var in x:
            if var.vtype == "CONTINUOUS":
                return False
        return True

    # main part of the solution process:
    main = Model("tsp")
    x = {}
    for i in V:
        for j in V:
            if j > i:
                x[i,j] = main.addVar(ub=1, vtype="C", name="x(%s,%s)"%(i,j))

    for i in V:
        main.addCons(quicksum(x[j,i] for j in V if j < i) + \
                       quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i)

    main.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j > i), "minimize")

    while True:
        main.optimize()
        z = main.getObjVal()
        X = {}
        for (i,j) in x:
            if main.getVal(x[i,j]) > EPS:
                X[i,j] = main.getVal(x[i,j])

        if addcut(X) == False:  # i.e., components are connected
            if isMIP():      # integer variables, components connected: solution found
                break
            for (i,j) in x:     # all components connected, switch to integer model
                main.chgVarType(x[i,j], "BINARY")

    # process solution
    edges = []
    for (i,j) in x:
        if main.getVal(x[i,j]) > EPS:
            edges.append((i,j))
    return main.getObjVal(),edges