#!/usr/bin/env python

import sys
import numpy
from pulp import *

capacities = numpy.loadtxt(sys.argv[1])
demands = numpy.loadtxt(sys.argv[2])
n = len(capacities)
neighbours = [[v for v in xrange(n) if capacities[u, v] > 0] 
        for u in xrange(n)]

model = LpProblem('Integer multi-flow')
flow = LpVariable.dicts('Flow', [(s, t, u, v)
    for s in xrange(n) for t in xrange(n) 
    for u in xrange(n) for v in neighbours[u]], 
    lowBound = 0, cat = LpInteger)
for s in xrange(n):
  for t in xrange(n):
    for u in xrange(n):
      model += lpSum(flow[s, t, v, u] for v in neighbours[u])\
              -lpSum(flow[s, t, u, z] for z in neighbours[u])\
               ==  (-demands[s, t] if u == s \
                    else demands[s, t] if u == t \
                    else 0)
for u in xrange(n):
    for v in neighbours[u]:
        if v < u: # Constrain each edge only once
            model += lpSum(flow[s, t, v, u] + flow[s, t, u, v]
                    for s in xrange(n) for t in xrange(n)) <= capacities[v, u]

model.writeLP('intflow.lp')
pulp.CPLEX_CMD().solve(model)
print "\nRouting is%s feasible" % \
        ('' if LpStatus[model.status] == 'Optimal' else ' not')

import networkx as nx
import matplotlib.pyplot as plt
from time import sleep

net = nx.Graph(capacities)
pos = nx.spring_layout(net)
nx.draw(net, pos = pos)
plt.draw(); sleep(3); plt.clf()

for s in xrange(n):
  for t in xrange(n):
    if demands[s, t] > 0:
      nx.draw(net, pos = pos)
      route = nx.Graph([(u, v) for u in xrange(n) for v in neighbours[u]
        if flow[s, t, u, v].value() > 0])
      nx.draw(route, pos = pos, edge_color = 'green', width = 5)
      plt.draw(); sleep(2); plt.clf()