lib/Documentation/min-cost-max-flow.cpp

   1 //Robado de shygypsy.com/tools/mcmf4.cpp
   2
   3 /**
   4  *   ///////////////////////
   5  *   // MIN COST MAX FLOW //
   6  *   ///////////////////////
   7  *
   8  *   Authors: Frank Chu, Igor Naverniouk
   9  **/
  10
  11 /*********************
  12  * Min cost max flow * (Edmonds-Karp relabelling + fast heap Dijkstra)
  13  *********************
  14  * Takes a directed graph where each edge has a capacity ('cap') and a
  15  * cost per unit of flow ('cost') and returns a maximum flow network
  16  * of minimal cost ('fcost') from s to t. USE mcmf3.cpp FOR DENSE GRAPHS!
  17  *
  18  * PARAMETERS:
  19  *      - cap (global): adjacency matrix where cap[u][v] is the capacity
  20  *          of the edge u->v. cap[u][v] is 0 for non-existent edges.
  21  *      - cost (global): a matrix where cost[u][v] is the cost per unit
  22  *          of flow along the edge u->v. If cap[u][v] == 0, cost[u][v] is
  23  *          ignored. ALL COSTS MUST BE NON-NEGATIVE!
  24  *      - n: the number of vertices ([0, n-1] are considered as vertices).
  25  *      - s: source vertex.
  26  *      - t: sink.
  27  * RETURNS:
  28  *      - the flow
  29  *      - the total cost through 'fcost'
  30  *      - fnet contains the flow network. Careful: both fnet[u][v] and
  31  *          fnet[v][u] could be positive. Take the difference.
  32  * COMPLEXITY:
  33  *      - Worst case: O(m*log(m)*flow  <?  n*m*log(m)*fcost)
  34  * FIELD TESTING:
  35  *      - Valladolid 10594: Data Flow
  36  * REFERENCE:
  37  *      Edmonds, J., Karp, R.  "Theoretical Improvements in Algorithmic
  38  *          Efficieincy for Network Flow Problems".
  39  *      This is a slight improvement of Frank Chu's implementation.
  40  **/
  41
  42 #include <iostream>
  43 using namespace std;
  44
  45 // the maximum number of vertices + 1
  46 #define NN 1024
  47
  48 // adjacency matrix (fill this up)
  49 int cap[NN][NN];
  50
  51 // cost per unit of flow matrix (fill this up)
  52 int cost[NN][NN];
  53
  54 // flow network and adjacency list
  55 int fnet[NN][NN], adj[NN][NN], deg[NN];
  56
  57 // Dijkstra's predecessor, depth and priority queue
  58 int par[NN], d[NN], q[NN], inq[NN], qs;
  59
  60 // Labelling function
  61 int pi[NN];
  62
  63 #define CLR(a, x) memset( a, x, sizeof( a ) )
  64 #define Inf (INT_MAX/2)
  65 #define BUBL { \
  66     t = q[i]; q[i] = q[j]; q[j] = t; \
  67     t = inq[q[i]]; inq[q[i]] = inq[q[j]]; inq[q[j]] = t; }
  68
  69 // Dijkstra's using non-negative edge weights (cost + potential)
  70 #define Pot(u,v) (d[u] + pi[u] - pi[v])
  71 bool dijkstra( int n, int s, int t )
  72 {
  73     CLR( d, 0x3F );
  74     CLR( par, -1 );
  75     CLR( inq, -1 );
  76     //for( int i = 0; i < n; i++ ) d[i] = Inf, par[i] = -1;
  77     d[s] = qs = 0;
  78     inq[q[qs++] = s] = 0;
  79     par[s] = n;
  80
  81     while( qs )
  82     {
  83         // get the minimum from q and bubble down
  84         int u = q[0]; inq[u] = -1;
  85         q[0] = q[--qs];
  86         if( qs ) inq[q[0]] = 0;
  87         for( int i = 0, j = 2*i + 1, t; j < qs; i = j, j = 2*i + 1 )
  88         {
  89             if( j + 1 < qs && d[q[j + 1]] < d[q[j]] ) j++;
  90             if( d[q[j]] >= d[q[i]] ) break;
  91             BUBL;
  92         }
  93
  94         // relax edge (u,i) or (i,u) for all i;
  95         for( int k = 0, v = adj[u][k]; k < deg[u]; v = adj[u][++k] )
  96         {
  97             // try undoing edge v->u
  98             if( fnet[v][u] && d[v] > Pot(u,v) - cost[v][u] )
  99                 d[v] = Pot(u,v) - cost[v][par[v] = u];
 100
 101             // try using edge u->v
 102             if( fnet[u][v] < cap[u][v] && d[v] > Pot(u,v) + cost[u][v] )
 103                 d[v] = Pot(u,v) + cost[par[v] = u][v];
 104
 105             if( par[v] == u )
 106             {
 107                 // bubble up or decrease key
 108                 if( inq[v] < 0 ) { inq[q[qs] = v] = qs; qs++; }
 109                 for( int i = inq[v], j = ( i - 1 )/2, t;
 110                      d[q[i]] < d[q[j]]; i = j, j = ( i - 1 )/2 )
 111                      BUBL;
 112             }
 113         }
 114     }
 115
 116     for( int i = 0; i < n; i++ ) if( pi[i] < Inf ) pi[i] += d[i];
 117
 118     return par[t] >= 0;
 119 }
 120 #undef Pot
 121
 122 int mcmf4( int n, int s, int t, int &fcost )
 123 {
 124     // build the adjacency list
 125     CLR( deg, 0 );
 126     for( int i = 0; i < n; i++ )
 127     for( int j = 0; j < n; j++ )
 128         if( cap[i][j] || cap[j][i] ) adj[i][deg[i]++] = j;
 129
 130     CLR( fnet, 0 );
 131     CLR( pi, 0 );
 132     int flow = fcost = 0;
 133
 134     // repeatedly, find a cheapest path from s to t
 135     while( dijkstra( n, s, t ) )
 136     {
 137         // get the bottleneck capacity
 138         int bot = INT_MAX;
 139         for( int v = t, u = par[v]; v != s; u = par[v = u] )
 140             bot <?= fnet[v][u] ? fnet[v][u] : ( cap[u][v] - fnet[u][v] );
 141
 142         // update the flow network
 143         for( int v = t, u = par[v]; v != s; u = par[v = u] )
 144             if( fnet[v][u] ) { fnet[v][u] -= bot; fcost -= bot * cost[v][u]; }
 145             else { fnet[u][v] += bot; fcost += bot * cost[u][v]; }
 146
 147         flow += bot;
 148     }
 149
 150     return flow;
 151 }
 152
 153 //----------------- EXAMPLE USAGE -----------------
 154 #include <iostream>
 155 #include <stdio.h>
 156 using namespace std;
 157
 158 int main()
 159 {
 160   int numV;
 161   //  while ( cin >> numV && numV ) {
 162   cin >> numV;
 163     memset( cap, 0, sizeof( cap ) );
 164
 165     int m, a, b, c, cp;
 166     int s, t;
 167     cin >> m;
 168     cin >> s >> t;
 169
 170     // fill up cap with existing capacities.
 171     // if the edge u->v has capacity 6, set cap[u][v] = 6.
 172     // for each cap[u][v] > 0, set cost[u][v] to  the
 173     // cost per unit of flow along the edge i->v
 174     for (int i=0; i<m; i++) {
 175       cin >> a >> b >> cp >> c;
 176       cost[a][b] = c; // cost[b][a] = c;
 177       cap[a][b] = cp; // cap[b][a] = cp;
 178     }
 179
 180     int fcost;
 181     int flow = mcmf3( numV, s, t, fcost );
 182     cout << "flow: " << flow << endl;
 183     cout << "cost: " << fcost << endl;
 184
 185     return 0;
 186 }