In [1]: 
import networkx as nx
import random as rnd
import matplotlib.pyplot as plt
from tabulate import tabulate
import math
from IPython.display import SVG
import pygraphviz as pgv

Matrix Rounding

Definition Feasible matrix rounding

  • Given a $m \times n$ matrix $M = \{a_{ij} \}$ of rational numbers.
  • $r_i = \sum_{j=1}^n a_{ij}$,for $i \in 1, \ldots,m $
  • $c_j=\sum_{i=1}^m a_{ij}$, for $j \in 1, \ldots,n $,
  • Round each element $a_{ij}$, $r_i$ and $c_j$ up or down to integer so that sum of rounded elements in each row (column) equal row (column) sum.

Example

Original matrix $$\begin{array}{c|ccc} &17.91 & 16.79 & 22.32\\ \hline 23.61 & 9.13 & 9.75 & 4.73\\ 14.74 & 0.42 & 5.23 & 9.09\\ 18.67 & 8.36 & 1.81 & 8.5\\ \end{array}$$

Rounded matrix $$\begin{array}{c|ccc} & 17 & 16 & 22\\ \hline 23 & 9 & 10 & 4\\ 14 & 0 & 5 & 9\\ 18 & 8 & 1 & 9 \\ \end{array}$$

Theorem

For any matrix $M$ there exists a feasible rounding

Max flow formulation

Instance generation

The random package is used to generate a random matrix

In [2]: 
numrows = 4
numcols = 5

matrix = [[round(rnd.random()*10,2) for x in range(numcols)] for x in range(numrows)] 

sumcol = [0 for j in range(numcols)]

for j in range(numcols):
    for i in range(numrows):
        sumcol [j] += matrix [i][j]

sumrow = [0 for i in range(numrows)]

for i in range(numrows):
    for j in range(numcols):
        sumrow [i] += matrix[i][j]

print "Original matrix"
print tabulate(matrix, tablefmt='grid')
        
print
print "Matrix with columns and rows sum"
table = []
for i in range(numrows):
    newrow = matrix[i][:]
    newrow.insert(0,sumrow[i])
    table.append(newrow)

print tabulate(table, headers = sumcol, tablefmt='grid')

Original matrix
+------+------+------+------+------+
| 0.51 | 9.84 | 7.27 | 0.09 | 9.02 |
+------+------+------+------+------+
| 3.5  | 8.02 | 2.99 | 5.43 | 3.31 |
+------+------+------+------+------+
| 0.89 | 9.18 | 0.8  | 0.8  | 6.31 |
+------+------+------+------+------+
| 4.44 | 0.3  | 1.43 | 7.19 | 9.13 |
+------+------+------+------+------+

Matrix with columns and rows sum
+-------+--------+---------+---------+---------+---------+
|       |   9.34 |   27.34 |   12.49 |   13.51 |   27.77 |
+=======+========+=========+=========+=========+=========+
| 26.73 |   0.51 |    9.84 |    7.27 |    0.09 |    9.02 |
+-------+--------+---------+---------+---------+---------+
| 23.25 |   3.5  |    8.02 |    2.99 |    5.43 |    3.31 |
+-------+--------+---------+---------+---------+---------+
| 17.98 |   0.89 |    9.18 |    0.8  |    0.8  |    6.31 |
+-------+--------+---------+---------+---------+---------+
| 22.49 |   4.44 |    0.3  |    1.43 |    7.19 |    9.13 |
+-------+--------+---------+---------+---------+---------+
In [3]: 
G = nx.DiGraph()

G.add_node('s')
G.add_node('t')

nodecol = ['c'+ str(i) for i in range(numcols)]
noderow = ['r'+ str(i) for i in range(numrows)]

G.add_nodes_from(nodecol)
G.add_nodes_from(noderow)

for i in range(numrows):
    G.add_edge ('s',noderow[i],lb=int (math.floor(sumrow[i])), ub= int(math.ceil(sumrow[i])))
    
for j in range(numcols):
        G.add_edge (nodecol[j],'t',lb=int (math.floor(sumcol[j])), ub= int(math.ceil(sumcol[j])))

for i in range(numrows):
    for j in range(numcols):
        G.add_edge (noderow[i],nodecol[j],lb=int (math.floor(matrix[i][j])), ub= int(math.ceil(matrix[i][j])))
In [4]: 
offset = 0.45
count = 0

lenghtcol = 135.0 * float(len(nodecol))
lenghtrow = 135.0 * float(len(noderow))
lenght = max (lenghtcol, lenghtrow)

offsetcol = lenght / len(nodecol)
offsetrow = lenght / len(noderow)

for i in noderow:
    G.node[i]['pos'] = "%f,%f"%(200, offsetrow * count)
    count += 1

count = 0

for j in nodecol:
    G.node[j]['pos'] = "%f,%f"%(500, offsetcol * count)
    count += 1
    
G.node['s']['pos'] = "%f,%f"%(0.0, offset * count * 300 / 2.0)
G.node['t']['pos'] = "%f,%f"%(700, offset * count * 300 /2.0)
In [5]: 
for i in G.edges():
    G[i[0]][i[1]]['xlabel'] = '%d,%d'%(G[i[0]][i[1]]['lb'],G[i[0]][i[1]]['ub'])

d = nx.to_agraph(G)

d.node_attr.update (fontsize='10', width=0.4, shape='circle')

d.edge_attr.update(fontsize='10', arrowhead='vee', penwidth=0.5)
#d.node_attr['shape']='circle


d.draw ('img.svg', prog='neato', args='-n2')

SVG(filename='img.svg')
#Image ('img.png')
Out[5]:
%3 r0 r0 c3 c3 r0->c3 0,1 c2 c2 r0->c2 7,8 c1 c1 r0->c1 9,10 c0 c0 r0->c0 0,1 c4 c4 r0->c4 9,10 t t c3->t 13,14 c2->t 12,13 c1->t 27,28 c0->t 9,10 c4->t 27,28 r1 r1 r1->c3 5,6 r1->c2 2,3 r1->c1 8,9 r1->c0 3,4 r1->c4 3,4 r2 r2 r2->c3 0,1 r2->c2 0,1 r2->c1 9,10 r2->c0 0,1 r2->c4 6,7 r3 r3 r3->c3 7,8 r3->c2 1,2 r3->c1 0,1 r3->c0 4,5 r3->c4 9,10 s s s->r0 26,27 s->r1 23,24 s->r2 17,18 s->r3 22,23

Finding a feasible rounding

Step 1. Generate a copy $H$ of the graph $G$

Remember The copy() method makes a complete copy of the graph including all of the node or edge attributes.

In [6]: 
H = G.copy()

Step 2. Capacity scaling: capacity of each arc of $H$ is scaled to $u_{ij} - l_{ij}$

In [7]: 
for j in H.edges():
    H[j[0]][j[1]]['capacity'] = \
    H[j[0]][j[1]]['ub'] -  H[j[0]][j[1]]['lb']

Step 3. Add an arc between $t$ and $s$ with infinite capacity (i.e., a capacity large enough)

Remark In this case you can assign to the arc $(t,s)$ a capacity equal to the maximum between $u(\delta^+(s))$ and $u(\delta^-(t))$

In [8]: 
row_cap = 0
for i in noderow:
    row_cap +=  H['s'][i]['ub'] 

col_cap = 0
for i in nodecol:
    col_cap +=  H[i]['t']['ub']

H.add_edge ('t','s',capacity=max(row_cap,col_cap))

Step 4. Two extra nodes $s_1$ and $t_1$ are added to $H$

In [9]: 
H.add_node('s1')
H.add_node('t1')

Step 5. For each node, the flow unbalance is evaluated

Remark Flow unbalance is stored in a dictionary and is evaluated on graph $G$. By constructions, the sum of flow unbalance is equal to 0

In [10]: 
unbalance = {}

for i in G.nodes():
        auxunb = 0
        for j in G.in_edges(i):
            auxunb += G[j[0]][j[1]]['lb']
        
        for j in G.out_edges(i):
            auxunb -= G[j[0]][j[1]]['lb']
        
        unbalance[i] = auxunb

Step 6. Two set of arcs are added to $H$:

  1. Arcs $(s_1, i)$ if flow unbalance of $i>0$
  2. Arcs $(i, t_1)$ if flow unbalance of $i<0$
In [11]: 
for i in unbalance:

    if unbalance[i] > 0:
        H.add_edge ('s1',i,capacity = unbalance[i])
    if unbalance[i] < 0:
        H.add_edge (i,'t1',capacity =- unbalance[i])

Step 7. Evaluate the maximum $s_1-t_1$ flow

A valid circulation always exists

In [12]: 
value, flow = nx.maximum_flow(H,'s1','t1', 'capacity')

for i in H.out_edges('s1'):
    print "Flow and capacity of arc (%s,%s):" % (i[0],i[1]), H[i[0]][i[1]]['capacity'], flow[i[0]][i[1]]
    if H[i[0]][i[1]]['capacity'] != flow[i[0]][i[1]]:
        print "Valid circulation not found"
        break
    

Flow and capacity of arc (s1,t): 88 88
Flow and capacity of arc (s1,r0): 1 1
Flow and capacity of arc (s1,r1): 2 2
Flow and capacity of arc (s1,r2): 2 2
Flow and capacity of arc (s1,r3): 1 1

Step 8. Print the rounded matrix

In [13]: 
rounded_matrix = []

for i in range(numrows):
    rounded_matrix.append([])
    for j in range(numcols):
        rounded_matrix[i].append(int(math.floor(matrix[i][j])) + \
        flow[noderow[i]][nodecol[j]])

sumroundedcol = [0 for j in range(numcols)]

for j in range(numcols):
    for i in range(numrows):
        sumroundedcol [j] += rounded_matrix [i][j]

sumroundedrow = [0 for i in range(numrows)]

for i in range(numrows):
    for j in range(numcols):
        sumroundedrow [i] += rounded_matrix[i][j]

print "Original matrix"
table = []
for i in range(numrows):
    newrow = matrix[i][:]
    newrow.insert(0,sumrow[i])
    table.append(newrow)

print tabulate(table, headers = sumcol, tablefmt="grid")

print
print "Rounded matrix"
table = []
for i in range(numrows):
    newrow = rounded_matrix[i][:]
    newrow.insert(0,sumroundedrow[i])
    table.append(newrow)

print tabulate(table, headers = sumroundedcol, tablefmt="grid")

Original matrix
+-------+--------+---------+---------+---------+---------+
|       |   9.34 |   27.34 |   12.49 |   13.51 |   27.77 |
+=======+========+=========+=========+=========+=========+
| 26.73 |   0.51 |    9.84 |    7.27 |    0.09 |    9.02 |
+-------+--------+---------+---------+---------+---------+
| 23.25 |   3.5  |    8.02 |    2.99 |    5.43 |    3.31 |
+-------+--------+---------+---------+---------+---------+
| 17.98 |   0.89 |    9.18 |    0.8  |    0.8  |    6.31 |
+-------+--------+---------+---------+---------+---------+
| 22.49 |   4.44 |    0.3  |    1.43 |    7.19 |    9.13 |
+-------+--------+---------+---------+---------+---------+

Rounded matrix
+----+-----+------+------+------+------+
|    |   9 |   27 |   12 |   13 |   27 |
+====+=====+======+======+======+======+
| 26 |   1 |    9 |    7 |    0 |    9 |
+----+-----+------+------+------+------+
| 23 |   4 |    8 |    3 |    5 |    3 |
+----+-----+------+------+------+------+
| 17 |   0 |   10 |    1 |    0 |    6 |
+----+-----+------+------+------+------+
| 22 |   4 |    0 |    1 |    8 |    9 |
+----+-----+------+------+------+------+