lemon/cost_scaling.h
changeset 859 ac5f72c48367
parent 831 cc9e0c15d747
parent 839 f3bc4e9b5f3a
child 863 a93f1a27d831
equal deleted inserted replaced
9:9fdcd1925b95 11:2d6029b252e4
   199   private:
   199   private:
   200 
   200 
   201     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   201     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   202 
   202 
   203     typedef std::vector<int> IntVector;
   203     typedef std::vector<int> IntVector;
   204     typedef std::vector<char> BoolVector;
       
   205     typedef std::vector<Value> ValueVector;
   204     typedef std::vector<Value> ValueVector;
   206     typedef std::vector<Cost> CostVector;
   205     typedef std::vector<Cost> CostVector;
   207     typedef std::vector<LargeCost> LargeCostVector;
   206     typedef std::vector<LargeCost> LargeCostVector;
       
   207     typedef std::vector<char> BoolVector;
       
   208     // Note: vector<char> is used instead of vector<bool> for efficiency reasons
   208 
   209 
   209   private:
   210   private:
   210   
   211   
   211     template <typename KT, typename VT>
   212     template <typename KT, typename VT>
   212     class StaticVectorMap {
   213     class StaticVectorMap {
   246     int _root;
   247     int _root;
   247 
   248 
   248     // Parameters of the problem
   249     // Parameters of the problem
   249     bool _have_lower;
   250     bool _have_lower;
   250     Value _sum_supply;
   251     Value _sum_supply;
       
   252     int _sup_node_num;
   251 
   253 
   252     // Data structures for storing the digraph
   254     // Data structures for storing the digraph
   253     IntNodeMap _node_id;
   255     IntNodeMap _node_id;
   254     IntArcMap _arc_idf;
   256     IntArcMap _arc_idf;
   255     IntArcMap _arc_idb;
   257     IntArcMap _arc_idb;
   274 
   276 
   275     // Data for scaling
   277     // Data for scaling
   276     LargeCost _epsilon;
   278     LargeCost _epsilon;
   277     int _alpha;
   279     int _alpha;
   278 
   280 
       
   281     IntVector _buckets;
       
   282     IntVector _bucket_next;
       
   283     IntVector _bucket_prev;
       
   284     IntVector _rank;
       
   285     int _max_rank;
       
   286   
   279     // Data for a StaticDigraph structure
   287     // Data for a StaticDigraph structure
   280     typedef std::pair<int, int> IntPair;
   288     typedef std::pair<int, int> IntPair;
   281     StaticDigraph _sgr;
   289     StaticDigraph _sgr;
   282     std::vector<IntPair> _arc_vec;
   290     std::vector<IntPair> _arc_vec;
   283     std::vector<LargeCost> _cost_vec;
   291     std::vector<LargeCost> _cost_vec;
   826         for (ArcIt a(_graph); a != INVALID; ++a) {
   834         for (ArcIt a(_graph); a != INVALID; ++a) {
   827           cap[a] = _upper[_arc_idf[a]];
   835           cap[a] = _upper[_arc_idf[a]];
   828         }
   836         }
   829       }
   837       }
   830 
   838 
       
   839       _sup_node_num = 0;
       
   840       for (NodeIt n(_graph); n != INVALID; ++n) {
       
   841         if (sup[n] > 0) ++_sup_node_num;
       
   842       }
       
   843 
   831       // Find a feasible flow using Circulation
   844       // Find a feasible flow using Circulation
   832       Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
   845       Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
   833         circ(_graph, low, cap, sup);
   846         circ(_graph, low, cap, sup);
   834       if (!circ.flowMap(flow).run()) return INFEASIBLE;
   847       if (!circ.flowMap(flow).run()) return INFEASIBLE;
   835 
   848 
   860           _res_cap[_arc_idf[a]] = cap[a] - fa;
   873           _res_cap[_arc_idf[a]] = cap[a] - fa;
   861           _res_cap[_arc_idb[a]] = fa;
   874           _res_cap[_arc_idb[a]] = fa;
   862         }
   875         }
   863         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   876         for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
   864           int ra = _reverse[a];
   877           int ra = _reverse[a];
   865           _res_cap[a] = 1;
   878           _res_cap[a] = 0;
   866           _res_cap[ra] = 0;
   879           _res_cap[ra] = 0;
   867           _cost[a] = 0;
   880           _cost[a] = 0;
   868           _cost[ra] = 0;
   881           _cost[ra] = 0;
   869         }
   882         }
   870       }
   883       }
   874 
   887 
   875     // Execute the algorithm and transform the results
   888     // Execute the algorithm and transform the results
   876     void start(Method method) {
   889     void start(Method method) {
   877       // Maximum path length for partial augment
   890       // Maximum path length for partial augment
   878       const int MAX_PATH_LENGTH = 4;
   891       const int MAX_PATH_LENGTH = 4;
   879       
   892 
       
   893       // Initialize data structures for buckets      
       
   894       _max_rank = _alpha * _res_node_num;
       
   895       _buckets.resize(_max_rank);
       
   896       _bucket_next.resize(_res_node_num + 1);
       
   897       _bucket_prev.resize(_res_node_num + 1);
       
   898       _rank.resize(_res_node_num + 1);
       
   899   
   880       // Execute the algorithm
   900       // Execute the algorithm
   881       switch (method) {
   901       switch (method) {
   882         case PUSH:
   902         case PUSH:
   883           startPush();
   903           startPush();
   884           break;
   904           break;
   913         for (int j = 0; j != limit; ++j) {
   933         for (int j = 0; j != limit; ++j) {
   914           if (!_forward[j]) _res_cap[j] += _lower[j];
   934           if (!_forward[j]) _res_cap[j] += _lower[j];
   915         }
   935         }
   916       }
   936       }
   917     }
   937     }
       
   938     
       
   939     // Initialize a cost scaling phase
       
   940     void initPhase() {
       
   941       // Saturate arcs not satisfying the optimality condition
       
   942       for (int u = 0; u != _res_node_num; ++u) {
       
   943         int last_out = _first_out[u+1];
       
   944         LargeCost pi_u = _pi[u];
       
   945         for (int a = _first_out[u]; a != last_out; ++a) {
       
   946           int v = _target[a];
       
   947           if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
       
   948             Value delta = _res_cap[a];
       
   949             _excess[u] -= delta;
       
   950             _excess[v] += delta;
       
   951             _res_cap[a] = 0;
       
   952             _res_cap[_reverse[a]] += delta;
       
   953           }
       
   954         }
       
   955       }
       
   956       
       
   957       // Find active nodes (i.e. nodes with positive excess)
       
   958       for (int u = 0; u != _res_node_num; ++u) {
       
   959         if (_excess[u] > 0) _active_nodes.push_back(u);
       
   960       }
       
   961 
       
   962       // Initialize the next arcs
       
   963       for (int u = 0; u != _res_node_num; ++u) {
       
   964         _next_out[u] = _first_out[u];
       
   965       }
       
   966     }
       
   967     
       
   968     // Early termination heuristic
       
   969     bool earlyTermination() {
       
   970       const double EARLY_TERM_FACTOR = 3.0;
       
   971 
       
   972       // Build a static residual graph
       
   973       _arc_vec.clear();
       
   974       _cost_vec.clear();
       
   975       for (int j = 0; j != _res_arc_num; ++j) {
       
   976         if (_res_cap[j] > 0) {
       
   977           _arc_vec.push_back(IntPair(_source[j], _target[j]));
       
   978           _cost_vec.push_back(_cost[j] + 1);
       
   979         }
       
   980       }
       
   981       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
       
   982 
       
   983       // Run Bellman-Ford algorithm to check if the current flow is optimal
       
   984       BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
       
   985       bf.init(0);
       
   986       bool done = false;
       
   987       int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
       
   988       for (int i = 0; i < K && !done; ++i) {
       
   989         done = bf.processNextWeakRound();
       
   990       }
       
   991       return done;
       
   992     }
       
   993 
       
   994     // Global potential update heuristic
       
   995     void globalUpdate() {
       
   996       int bucket_end = _root + 1;
       
   997     
       
   998       // Initialize buckets
       
   999       for (int r = 0; r != _max_rank; ++r) {
       
  1000         _buckets[r] = bucket_end;
       
  1001       }
       
  1002       Value total_excess = 0;
       
  1003       for (int i = 0; i != _res_node_num; ++i) {
       
  1004         if (_excess[i] < 0) {
       
  1005           _rank[i] = 0;
       
  1006           _bucket_next[i] = _buckets[0];
       
  1007           _bucket_prev[_buckets[0]] = i;
       
  1008           _buckets[0] = i;
       
  1009         } else {
       
  1010           total_excess += _excess[i];
       
  1011           _rank[i] = _max_rank;
       
  1012         }
       
  1013       }
       
  1014       if (total_excess == 0) return;
       
  1015 
       
  1016       // Search the buckets
       
  1017       int r = 0;
       
  1018       for ( ; r != _max_rank; ++r) {
       
  1019         while (_buckets[r] != bucket_end) {
       
  1020           // Remove the first node from the current bucket
       
  1021           int u = _buckets[r];
       
  1022           _buckets[r] = _bucket_next[u];
       
  1023           
       
  1024           // Search the incomming arcs of u
       
  1025           LargeCost pi_u = _pi[u];
       
  1026           int last_out = _first_out[u+1];
       
  1027           for (int a = _first_out[u]; a != last_out; ++a) {
       
  1028             int ra = _reverse[a];
       
  1029             if (_res_cap[ra] > 0) {
       
  1030               int v = _source[ra];
       
  1031               int old_rank_v = _rank[v];
       
  1032               if (r < old_rank_v) {
       
  1033                 // Compute the new rank of v
       
  1034                 LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
       
  1035                 int new_rank_v = old_rank_v;
       
  1036                 if (nrc < LargeCost(_max_rank))
       
  1037                   new_rank_v = r + 1 + int(nrc);
       
  1038                   
       
  1039                 // Change the rank of v
       
  1040                 if (new_rank_v < old_rank_v) {
       
  1041                   _rank[v] = new_rank_v;
       
  1042                   _next_out[v] = _first_out[v];
       
  1043                   
       
  1044                   // Remove v from its old bucket
       
  1045                   if (old_rank_v < _max_rank) {
       
  1046                     if (_buckets[old_rank_v] == v) {
       
  1047                       _buckets[old_rank_v] = _bucket_next[v];
       
  1048                     } else {
       
  1049                       _bucket_next[_bucket_prev[v]] = _bucket_next[v];
       
  1050                       _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
       
  1051                     }
       
  1052                   }
       
  1053                   
       
  1054                   // Insert v to its new bucket
       
  1055                   _bucket_next[v] = _buckets[new_rank_v];
       
  1056                   _bucket_prev[_buckets[new_rank_v]] = v;
       
  1057                   _buckets[new_rank_v] = v;
       
  1058                 }
       
  1059               }
       
  1060             }
       
  1061           }
       
  1062 
       
  1063           // Finish search if there are no more active nodes
       
  1064           if (_excess[u] > 0) {
       
  1065             total_excess -= _excess[u];
       
  1066             if (total_excess <= 0) break;
       
  1067           }
       
  1068         }
       
  1069         if (total_excess <= 0) break;
       
  1070       }
       
  1071       
       
  1072       // Relabel nodes
       
  1073       for (int u = 0; u != _res_node_num; ++u) {
       
  1074         int k = std::min(_rank[u], r);
       
  1075         if (k > 0) {
       
  1076           _pi[u] -= _epsilon * k;
       
  1077           _next_out[u] = _first_out[u];
       
  1078         }
       
  1079       }
       
  1080     }
   918 
  1081 
   919     /// Execute the algorithm performing augment and relabel operations
  1082     /// Execute the algorithm performing augment and relabel operations
   920     void startAugment(int max_length = std::numeric_limits<int>::max()) {
  1083     void startAugment(int max_length = std::numeric_limits<int>::max()) {
   921       // Paramters for heuristics
  1084       // Paramters for heuristics
   922       const int BF_HEURISTIC_EPSILON_BOUND = 1000;
  1085       const int EARLY_TERM_EPSILON_LIMIT = 1000;
   923       const int BF_HEURISTIC_BOUND_FACTOR  = 3;
  1086       const double GLOBAL_UPDATE_FACTOR = 3.0;
   924 
  1087 
       
  1088       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
       
  1089         (_res_node_num + _sup_node_num * _sup_node_num));
       
  1090       int next_update_limit = global_update_freq;
       
  1091       
       
  1092       int relabel_cnt = 0;
       
  1093       
   925       // Perform cost scaling phases
  1094       // Perform cost scaling phases
   926       IntVector pred_arc(_res_node_num);
  1095       std::vector<int> path;
   927       std::vector<int> path_nodes;
       
   928       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
  1096       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
   929                                         1 : _epsilon / _alpha )
  1097                                         1 : _epsilon / _alpha )
   930       {
  1098       {
   931         // "Early Termination" heuristic: use Bellman-Ford algorithm
  1099         // Early termination heuristic
   932         // to check if the current flow is optimal
  1100         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
   933         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
  1101           if (earlyTermination()) break;
   934           _arc_vec.clear();
       
   935           _cost_vec.clear();
       
   936           for (int j = 0; j != _res_arc_num; ++j) {
       
   937             if (_res_cap[j] > 0) {
       
   938               _arc_vec.push_back(IntPair(_source[j], _target[j]));
       
   939               _cost_vec.push_back(_cost[j] + 1);
       
   940             }
       
   941           }
       
   942           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
       
   943 
       
   944           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
       
   945           bf.init(0);
       
   946           bool done = false;
       
   947           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
       
   948           for (int i = 0; i < K && !done; ++i)
       
   949             done = bf.processNextWeakRound();
       
   950           if (done) break;
       
   951         }
       
   952 
       
   953         // Saturate arcs not satisfying the optimality condition
       
   954         for (int a = 0; a != _res_arc_num; ++a) {
       
   955           if (_res_cap[a] > 0 &&
       
   956               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
       
   957             Value delta = _res_cap[a];
       
   958             _excess[_source[a]] -= delta;
       
   959             _excess[_target[a]] += delta;
       
   960             _res_cap[a] = 0;
       
   961             _res_cap[_reverse[a]] += delta;
       
   962           }
       
   963         }
  1102         }
   964         
  1103         
   965         // Find active nodes (i.e. nodes with positive excess)
  1104         // Initialize current phase
   966         for (int u = 0; u != _res_node_num; ++u) {
  1105         initPhase();
   967           if (_excess[u] > 0) _active_nodes.push_back(u);
  1106         
   968         }
       
   969 
       
   970         // Initialize the next arcs
       
   971         for (int u = 0; u != _res_node_num; ++u) {
       
   972           _next_out[u] = _first_out[u];
       
   973         }
       
   974 
       
   975         // Perform partial augment and relabel operations
  1107         // Perform partial augment and relabel operations
   976         while (true) {
  1108         while (true) {
   977           // Select an active node (FIFO selection)
  1109           // Select an active node (FIFO selection)
   978           while (_active_nodes.size() > 0 &&
  1110           while (_active_nodes.size() > 0 &&
   979                  _excess[_active_nodes.front()] <= 0) {
  1111                  _excess[_active_nodes.front()] <= 0) {
   980             _active_nodes.pop_front();
  1112             _active_nodes.pop_front();
   981           }
  1113           }
   982           if (_active_nodes.size() == 0) break;
  1114           if (_active_nodes.size() == 0) break;
   983           int start = _active_nodes.front();
  1115           int start = _active_nodes.front();
   984           path_nodes.clear();
       
   985           path_nodes.push_back(start);
       
   986 
  1116 
   987           // Find an augmenting path from the start node
  1117           // Find an augmenting path from the start node
       
  1118           path.clear();
   988           int tip = start;
  1119           int tip = start;
   989           while (_excess[tip] >= 0 &&
  1120           while (_excess[tip] >= 0 && int(path.size()) < max_length) {
   990                  int(path_nodes.size()) <= max_length) {
       
   991             int u;
  1121             int u;
   992             LargeCost min_red_cost, rc;
  1122             LargeCost min_red_cost, rc, pi_tip = _pi[tip];
   993             int last_out = _sum_supply < 0 ?
  1123             int last_out = _first_out[tip+1];
   994               _first_out[tip+1] : _first_out[tip+1] - 1;
       
   995             for (int a = _next_out[tip]; a != last_out; ++a) {
  1124             for (int a = _next_out[tip]; a != last_out; ++a) {
   996               if (_res_cap[a] > 0 &&
  1125               u = _target[a];
   997                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
  1126               if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
   998                 u = _target[a];
  1127                 path.push_back(a);
   999                 pred_arc[u] = a;
       
  1000                 _next_out[tip] = a;
  1128                 _next_out[tip] = a;
  1001                 tip = u;
  1129                 tip = u;
  1002                 path_nodes.push_back(tip);
       
  1003                 goto next_step;
  1130                 goto next_step;
  1004               }
  1131               }
  1005             }
  1132             }
  1006 
  1133 
  1007             // Relabel tip node
  1134             // Relabel tip node
  1008             min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
  1135             min_red_cost = std::numeric_limits<LargeCost>::max();
       
  1136             if (tip != start) {
       
  1137               int ra = _reverse[path.back()];
       
  1138               min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
       
  1139             }
  1009             for (int a = _first_out[tip]; a != last_out; ++a) {
  1140             for (int a = _first_out[tip]; a != last_out; ++a) {
  1010               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
  1141               rc = _cost[a] + pi_tip - _pi[_target[a]];
  1011               if (_res_cap[a] > 0 && rc < min_red_cost) {
  1142               if (_res_cap[a] > 0 && rc < min_red_cost) {
  1012                 min_red_cost = rc;
  1143                 min_red_cost = rc;
  1013               }
  1144               }
  1014             }
  1145             }
  1015             _pi[tip] -= min_red_cost + _epsilon;
  1146             _pi[tip] -= min_red_cost + _epsilon;
  1016 
       
  1017             // Reset the next arc of tip
       
  1018             _next_out[tip] = _first_out[tip];
  1147             _next_out[tip] = _first_out[tip];
       
  1148             ++relabel_cnt;
  1019 
  1149 
  1020             // Step back
  1150             // Step back
  1021             if (tip != start) {
  1151             if (tip != start) {
  1022               path_nodes.pop_back();
  1152               tip = _source[path.back()];
  1023               tip = path_nodes.back();
  1153               path.pop_back();
  1024             }
  1154             }
  1025 
  1155 
  1026           next_step: ;
  1156           next_step: ;
  1027           }
  1157           }
  1028 
  1158 
  1029           // Augment along the found path (as much flow as possible)
  1159           // Augment along the found path (as much flow as possible)
  1030           Value delta;
  1160           Value delta;
  1031           int u, v = path_nodes.front(), pa;
  1161           int pa, u, v = start;
  1032           for (int i = 1; i < int(path_nodes.size()); ++i) {
  1162           for (int i = 0; i != int(path.size()); ++i) {
       
  1163             pa = path[i];
  1033             u = v;
  1164             u = v;
  1034             v = path_nodes[i];
  1165             v = _target[pa];
  1035             pa = pred_arc[v];
       
  1036             delta = std::min(_res_cap[pa], _excess[u]);
  1166             delta = std::min(_res_cap[pa], _excess[u]);
  1037             _res_cap[pa] -= delta;
  1167             _res_cap[pa] -= delta;
  1038             _res_cap[_reverse[pa]] += delta;
  1168             _res_cap[_reverse[pa]] += delta;
  1039             _excess[u] -= delta;
  1169             _excess[u] -= delta;
  1040             _excess[v] += delta;
  1170             _excess[v] += delta;
  1041             if (_excess[v] > 0 && _excess[v] <= delta)
  1171             if (_excess[v] > 0 && _excess[v] <= delta)
  1042               _active_nodes.push_back(v);
  1172               _active_nodes.push_back(v);
  1043           }
  1173           }
       
  1174 
       
  1175           // Global update heuristic
       
  1176           if (relabel_cnt >= next_update_limit) {
       
  1177             globalUpdate();
       
  1178             next_update_limit += global_update_freq;
       
  1179           }
  1044         }
  1180         }
  1045       }
  1181       }
  1046     }
  1182     }
  1047 
  1183 
  1048     /// Execute the algorithm performing push and relabel operations
  1184     /// Execute the algorithm performing push and relabel operations
  1049     void startPush() {
  1185     void startPush() {
  1050       // Paramters for heuristics
  1186       // Paramters for heuristics
  1051       const int BF_HEURISTIC_EPSILON_BOUND = 1000;
  1187       const int EARLY_TERM_EPSILON_LIMIT = 1000;
  1052       const int BF_HEURISTIC_BOUND_FACTOR  = 3;
  1188       const double GLOBAL_UPDATE_FACTOR = 2.0;
  1053 
  1189 
       
  1190       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
       
  1191         (_res_node_num + _sup_node_num * _sup_node_num));
       
  1192       int next_update_limit = global_update_freq;
       
  1193 
       
  1194       int relabel_cnt = 0;
       
  1195       
  1054       // Perform cost scaling phases
  1196       // Perform cost scaling phases
  1055       BoolVector hyper(_res_node_num, false);
  1197       BoolVector hyper(_res_node_num, false);
       
  1198       LargeCostVector hyper_cost(_res_node_num);
  1056       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
  1199       for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
  1057                                         1 : _epsilon / _alpha )
  1200                                         1 : _epsilon / _alpha )
  1058       {
  1201       {
  1059         // "Early Termination" heuristic: use Bellman-Ford algorithm
  1202         // Early termination heuristic
  1060         // to check if the current flow is optimal
  1203         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
  1061         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
  1204           if (earlyTermination()) break;
  1062           _arc_vec.clear();
  1205         }
  1063           _cost_vec.clear();
  1206         
  1064           for (int j = 0; j != _res_arc_num; ++j) {
  1207         // Initialize current phase
  1065             if (_res_cap[j] > 0) {
  1208         initPhase();
  1066               _arc_vec.push_back(IntPair(_source[j], _target[j]));
       
  1067               _cost_vec.push_back(_cost[j] + 1);
       
  1068             }
       
  1069           }
       
  1070           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
       
  1071 
       
  1072           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
       
  1073           bf.init(0);
       
  1074           bool done = false;
       
  1075           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
       
  1076           for (int i = 0; i < K && !done; ++i)
       
  1077             done = bf.processNextWeakRound();
       
  1078           if (done) break;
       
  1079         }
       
  1080 
       
  1081         // Saturate arcs not satisfying the optimality condition
       
  1082         for (int a = 0; a != _res_arc_num; ++a) {
       
  1083           if (_res_cap[a] > 0 &&
       
  1084               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
       
  1085             Value delta = _res_cap[a];
       
  1086             _excess[_source[a]] -= delta;
       
  1087             _excess[_target[a]] += delta;
       
  1088             _res_cap[a] = 0;
       
  1089             _res_cap[_reverse[a]] += delta;
       
  1090           }
       
  1091         }
       
  1092 
       
  1093         // Find active nodes (i.e. nodes with positive excess)
       
  1094         for (int u = 0; u != _res_node_num; ++u) {
       
  1095           if (_excess[u] > 0) _active_nodes.push_back(u);
       
  1096         }
       
  1097 
       
  1098         // Initialize the next arcs
       
  1099         for (int u = 0; u != _res_node_num; ++u) {
       
  1100           _next_out[u] = _first_out[u];
       
  1101         }
       
  1102 
  1209 
  1103         // Perform push and relabel operations
  1210         // Perform push and relabel operations
  1104         while (_active_nodes.size() > 0) {
  1211         while (_active_nodes.size() > 0) {
  1105           LargeCost min_red_cost, rc;
  1212           LargeCost min_red_cost, rc, pi_n;
  1106           Value delta;
  1213           Value delta;
  1107           int n, t, a, last_out = _res_arc_num;
  1214           int n, t, a, last_out = _res_arc_num;
  1108 
  1215 
       
  1216         next_node:
  1109           // Select an active node (FIFO selection)
  1217           // Select an active node (FIFO selection)
  1110         next_node:
       
  1111           n = _active_nodes.front();
  1218           n = _active_nodes.front();
  1112           last_out = _sum_supply < 0 ?
  1219           last_out = _first_out[n+1];
  1113             _first_out[n+1] : _first_out[n+1] - 1;
  1220           pi_n = _pi[n];
  1114 
  1221           
  1115           // Perform push operations if there are admissible arcs
  1222           // Perform push operations if there are admissible arcs
  1116           if (_excess[n] > 0) {
  1223           if (_excess[n] > 0) {
  1117             for (a = _next_out[n]; a != last_out; ++a) {
  1224             for (a = _next_out[n]; a != last_out; ++a) {
  1118               if (_res_cap[a] > 0 &&
  1225               if (_res_cap[a] > 0 &&
  1119                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
  1226                   _cost[a] + pi_n - _pi[_target[a]] < 0) {
  1120                 delta = std::min(_res_cap[a], _excess[n]);
  1227                 delta = std::min(_res_cap[a], _excess[n]);
  1121                 t = _target[a];
  1228                 t = _target[a];
  1122 
  1229 
  1123                 // Push-look-ahead heuristic
  1230                 // Push-look-ahead heuristic
  1124                 Value ahead = -_excess[t];
  1231                 Value ahead = -_excess[t];
  1125                 int last_out_t = _sum_supply < 0 ?
  1232                 int last_out_t = _first_out[t+1];
  1126                   _first_out[t+1] : _first_out[t+1] - 1;
  1233                 LargeCost pi_t = _pi[t];
  1127                 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
  1234                 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
  1128                   if (_res_cap[ta] > 0 && 
  1235                   if (_res_cap[ta] > 0 && 
  1129                       _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
  1236                       _cost[ta] + pi_t - _pi[_target[ta]] < 0)
  1130                     ahead += _res_cap[ta];
  1237                     ahead += _res_cap[ta];
  1131                   if (ahead >= delta) break;
  1238                   if (ahead >= delta) break;
  1132                 }
  1239                 }
  1133                 if (ahead < 0) ahead = 0;
  1240                 if (ahead < 0) ahead = 0;
  1134 
  1241 
  1135                 // Push flow along the arc
  1242                 // Push flow along the arc
  1136                 if (ahead < delta) {
  1243                 if (ahead < delta && !hyper[t]) {
  1137                   _res_cap[a] -= ahead;
  1244                   _res_cap[a] -= ahead;
  1138                   _res_cap[_reverse[a]] += ahead;
  1245                   _res_cap[_reverse[a]] += ahead;
  1139                   _excess[n] -= ahead;
  1246                   _excess[n] -= ahead;
  1140                   _excess[t] += ahead;
  1247                   _excess[t] += ahead;
  1141                   _active_nodes.push_front(t);
  1248                   _active_nodes.push_front(t);
  1142                   hyper[t] = true;
  1249                   hyper[t] = true;
       
  1250                   hyper_cost[t] = _cost[a] + pi_n - pi_t;
  1143                   _next_out[n] = a;
  1251                   _next_out[n] = a;
  1144                   goto next_node;
  1252                   goto next_node;
  1145                 } else {
  1253                 } else {
  1146                   _res_cap[a] -= delta;
  1254                   _res_cap[a] -= delta;
  1147                   _res_cap[_reverse[a]] += delta;
  1255                   _res_cap[_reverse[a]] += delta;
  1160             _next_out[n] = a;
  1268             _next_out[n] = a;
  1161           }
  1269           }
  1162 
  1270 
  1163           // Relabel the node if it is still active (or hyper)
  1271           // Relabel the node if it is still active (or hyper)
  1164           if (_excess[n] > 0 || hyper[n]) {
  1272           if (_excess[n] > 0 || hyper[n]) {
  1165             min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
  1273              min_red_cost = hyper[n] ? -hyper_cost[n] :
       
  1274                std::numeric_limits<LargeCost>::max();
  1166             for (int a = _first_out[n]; a != last_out; ++a) {
  1275             for (int a = _first_out[n]; a != last_out; ++a) {
  1167               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
  1276               rc = _cost[a] + pi_n - _pi[_target[a]];
  1168               if (_res_cap[a] > 0 && rc < min_red_cost) {
  1277               if (_res_cap[a] > 0 && rc < min_red_cost) {
  1169                 min_red_cost = rc;
  1278                 min_red_cost = rc;
  1170               }
  1279               }
  1171             }
  1280             }
  1172             _pi[n] -= min_red_cost + _epsilon;
  1281             _pi[n] -= min_red_cost + _epsilon;
       
  1282             _next_out[n] = _first_out[n];
  1173             hyper[n] = false;
  1283             hyper[n] = false;
  1174 
  1284             ++relabel_cnt;
  1175             // Reset the next arc
       
  1176             _next_out[n] = _first_out[n];
       
  1177           }
  1285           }
  1178         
  1286         
  1179           // Remove nodes that are not active nor hyper
  1287           // Remove nodes that are not active nor hyper
  1180         remove_nodes:
  1288         remove_nodes:
  1181           while ( _active_nodes.size() > 0 &&
  1289           while ( _active_nodes.size() > 0 &&
  1182                   _excess[_active_nodes.front()] <= 0 &&
  1290                   _excess[_active_nodes.front()] <= 0 &&
  1183                   !hyper[_active_nodes.front()] ) {
  1291                   !hyper[_active_nodes.front()] ) {
  1184             _active_nodes.pop_front();
  1292             _active_nodes.pop_front();
  1185           }
  1293           }
       
  1294           
       
  1295           // Global update heuristic
       
  1296           if (relabel_cnt >= next_update_limit) {
       
  1297             globalUpdate();
       
  1298             for (int u = 0; u != _res_node_num; ++u)
       
  1299               hyper[u] = false;
       
  1300             next_update_limit += global_update_freq;
       
  1301           }
  1186         }
  1302         }
  1187       }
  1303       }
  1188     }
  1304     }
  1189 
  1305 
  1190   }; //class CostScaling
  1306   }; //class CostScaling