mcf: multi-commodity flow formulation for the (asymmetric) traveling salesman problem
    Parameters:
        - n: number of nodes
        - c[i,j]: cost for traversing arc (i,j)
    Returns a model, ready to be solved.

def mcf(n,c):
    """mcf: multi-commodity flow formulation for the (asymmetric) traveling salesman problem
    Parameters:
        - n: number of nodes
        - c[i,j]: cost for traversing arc (i,j)
    Returns a model, ready to be solved.
    """
    model = Model("mcf")

    x,f = {},{}
    for i in range(1,n+1):
        for j in range(1,n+1):
            if i != j:
                x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))
            if i != j and j != 1:
                for k in range(2,n+1):
                    if i != k:
                        f[i,j,k] = model.addVar(ub=1, vtype="C", name="f(%s,%s,%s)"%(i,j,k))

    for i in range(1,n+1):
        model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i)
        model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i)

    for k in range(2,n+1):
        model.addCons(quicksum(f[1,i,k] for i in range(2,n+1) if (1,i,k) in f) == 1, "FlowOut(%s)"%k)
        model.addCons(quicksum(f[i,k,k] for i in range(1,n+1) if (i,k,k) in f) == 1, "FlowIn(%s)"%k)

        for i in range(2,n+1):
            if i != k:
                model.addCons(quicksum(f[j,i,k] for j in range(1,n+1) if (j,i,k) in f) == \
                                quicksum(f[i,j,k] for j in range(1,n+1) if (i,j,k) in f),
                                "FlowCons(%s,%s)"%(i,k))

    for (i,j,k) in f:
        model.addCons(f[i,j,k] <= x[i,j], "FlowUB(%s,%s,%s)"%(i,j,k))

    model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize")

    model.data = x,f
    return model