2d House Robber : Max Flow

import random
"""
2 Dimensional House Robber Problem

The 1D House Robber is : Given a street containing N houses labeled 1 to n
each containing valuables totaling v_i to be stolen. Your task is to rob houses
in order to maximize the amount of loot. However there is a security system
in place and if you rob two adjacent houses an alarm will go off.
Find the optimum set of non-adjacent houses to rob.
    e.g. Given a list of houses represented as their total value:
        houses    := [ 10, 5, 5, 20, 5]
        alignment := [  1, 0, 0,  1, 0]
        That is, selecting the first and fourth house maximizes profit.

This problem can be expanded into 2 dimensions, and is considerably more
difficult.

The 2D House Robber is : Given a City laid out on a grid, there are R rows
of houses and C columns of houses. Each house is known to contain some amount
of valuables totaling v_ij. Your task is to rob as many houses as possible
to maximize the amount of loot. However there is a security system in place
and if you rob two adjacent houses an alarm will go off
Find the optimum set of non-adjacent houses to rob.
    e.g.:  Houses:  alignment:
          10 20 10     0  1  0
          20 30 20 =>  1  0  1
          10 20 10     0  1  0
        This alignment results in the maximum of 80.

Solution using Dynamic Programming:
    The naive DP solution is O(R*2^C) but can be optimized to O(R*φ^C)
    The solution is exponential in N

Solution using LP:
    maximize c^Tx (read: c transpose x (dot product))
    define constraints:
         x + y <= 1 for each pair (x,y) in the grid graph
         xij >= 0 for all variables
    e.g. for a 2x2 grid:
            5  2
            8  4
        maximize :
            5*x11 + 2*x12 + 8*x21 + 4*x22
        subject to these constraints:
            x11 + x12 <= 1
            x11 + x21 <= 1
            x22 + x12 <= 1
            x22 + x21 <= 1
            x11,x12,x21,x22 >= 0
    The fact that the A^T matrix is Totally Unimodular means that the solution
    to this linear program has a linear solution, by duality. This proof
    is rather involved.

    I havent done out the big Oh notation for this but it looks like O(N^6).

Solution using Flow:

This is essentially a Maximum Bipartite matching problem.

The grid graph is bipartite, so separate the nodes into two sets, U and V.
Create the following flow network:
 *  create a vertex for the source S and sink T.
 *  for every node u in U with weight w, add an edge from S to u with capacity w
 *  for every node v in V with weight w, add an edge from v to T with capacity w
 *  for every edge (u,v) in the graph, add an edge from u to v with arbitrarily
        large weight, such that any min cut would be better off cutting all the
        other edges instead of any one of these.

The min cut will cut edges adjacent to S and T.
Suppose we have cut this graph.
Then for every edge (u,v) with large capacity, either (S,u) was cut,
or (u,T) was cut.
This cut is isomorphic to a selection of vertices.
Thus for every edge (u,v) either u was selected or v was selected.
Thus our min cut actually corresponds to a min vertex cover.
Once we have the min vertex cover, we can take its complement to get
the max weighted independent set.
"""

def printm(m):
    print("----"*len(m[0]))
    for row in m:
        print(' '.join("%3d"%x for x in row))
    print("----"*len(m[0]))

def flow_solver_update_main(C,tbl,ali,i,j,x,y):
    """ compute flow from s->u->v->t
    node u = (i,j) and node v=(x,y)
    """
    u = tbl[i][j]
    v = tbl[x][y]
    f = min(u,v)
    if f!=0:
        tbl[i][j] -= f
        tbl[x][y] -= f
        if tbl[x][y]==0:
            print("cut %d->t"%((C+1)*x+y+1))
            ali[x][y]=0
        if tbl[i][j]==0:
            print("cut s->%d"%((C+1)*i+j+1))
            ali[i][j]=0

    return f

def flow_solver_update(R,C,tbl,ali,i,j):
    """compute flow from path s to u (i,j) to each adjacent node v"""
    flow = 0
    if i>0:
        flow+=flow_solver_update_main(C,tbl,ali,i,j,i-1,j)
    if j>0:
        flow+=flow_solver_update_main(C,tbl,ali,i,j,i,j-1)
    if j<C:
        flow+=flow_solver_update_main(C,tbl,ali,i,j,i,j+1)
    if i<R:
        flow+=flow_solver_update_main(C,tbl,ali,i,j,i+1,j)
    return flow

def flow_solver(R,C,tbl):
    """return optimal alignment for a given 2d house-robber matrix"""

    printm(tbl)

    # start with an alignment which selects all indices
    # then eliminate any index corresponding to the min cut
    ali=[[1,]*C for i in range(R)]

    #Given:  Partition:  valid paths:
    # 1 2 3      u v u    s->1->2->t
    # 4 5 6  =>  v u v    s->1->4->t
    # 7 8 9      u v u    s->5->8->t (etc)
    # calculate flow passing through all nodes u in U.
    # this is effectivly Dinic's Algorithm for flow.
    # perform a depth first search s->u->v->t
    # all u->v links are implied connections based on locality in the grid.
    f = 0
    for i in range(R):
        for j in range(i%2,C,2):
            f+=flow_solver_update(R-1,C-1,tbl,ali,i,j)

    print("tbl:")
    printm(tbl)
    print("ali:")
    printm(ali)

    print("max flow / min cut :",f) # equal to the min cut

    return ali


def main():

    #
    # 1 2 3
    # 4 5 6
    # 7 8 9

    tbla = [
          [10,20,10], # min cut / max flow is 70
          [20,30,20],
          [10,20,10],
    ]

    tblb = [
          [25, 1,25], # min cut / max flow is 21
          [ 2, 3, 4],
          [ 5,25, 6],
    ]

    tblc = [
          [10,20,10,20,],
          [20,30,20,90,],
          [10,20,10,20,],
    ]

    #C = 5
    #R = 5
    #tblr = [ [random.randint(0,20) for c in range(C)] for r in range(R) ]

    tbl=tbla
    R = len(tbl)
    C = len(tbl[0])
    flow_solver(R,C,tbl)

if __name__ == '__main__':
    main()