COIN-OR::LEMON - Graph Library

Ticket #417: 417-ca9d1dc77026--78d9eddfaf1a.patch

File 417-ca9d1dc77026--78d9eddfaf1a.patch, 29.4 KB (added by Peter Kovacs, 9 years ago)
  • lemon/cost_scaling.h

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1300212980 -3600
    # Node ID 4866b640dba97fae2fdfdd81cf4bec497d1cd48d
    # Parent  ca9d1dc770268c0b476d33857ece1c13471ea447
    Fix and improve refine methods in CostScaling (#417)
    The old implementation could run into infinite loop using the
    AUGMENT method.
    
    diff --git a/lemon/cost_scaling.h b/lemon/cost_scaling.h
    a b  
    11001100    }
    11011101
    11021102    /// Execute the algorithm performing augment and relabel operations
    1103     void startAugment(int max_length = std::numeric_limits<int>::max()) {
     1103    void startAugment(int max_length) {
    11041104      // Paramters for heuristics
    11051105      const int EARLY_TERM_EPSILON_LIMIT = 1000;
    1106       const double GLOBAL_UPDATE_FACTOR = 3.0;
    1107 
    1108       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1106      const double GLOBAL_UPDATE_FACTOR = 1.0;
     1107      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
    11091108        (_res_node_num + _sup_node_num * _sup_node_num));
    1110       int next_update_limit = global_update_freq;
    1111 
    1112       int relabel_cnt = 0;
     1109      int next_global_update_limit = global_update_skip;
    11131110
    11141111      // Perform cost scaling phases
    1115       std::vector<int> path;
     1112      IntVector path;
     1113      BoolVector path_arc(_res_arc_num, false);
     1114      int relabel_cnt = 0;
    11161115      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    11171116                                        1 : _epsilon / _alpha )
    11181117      {
     
    11351134          int start = _active_nodes.front();
    11361135
    11371136          // Find an augmenting path from the start node
    1138           path.clear();
    11391137          int tip = start;
    1140           while (_excess[tip] >= 0 && int(path.size()) < max_length) {
     1138          while (int(path.size()) < max_length && _excess[tip] >= 0) {
    11411139            int u;
    1142             LargeCost min_red_cost, rc, pi_tip = _pi[tip];
     1140            LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
     1141            LargeCost pi_tip = _pi[tip];
    11431142            int last_out = _first_out[tip+1];
    11441143            for (int a = _next_out[tip]; a != last_out; ++a) {
    1145               u = _target[a];
    1146               if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
    1147                 path.push_back(a);
    1148                 _next_out[tip] = a;
    1149                 tip = u;
    1150                 goto next_step;
     1144              if (_res_cap[a] > 0) {
     1145                u = _target[a];
     1146                rc = _cost[a] + pi_tip - _pi[u];
     1147                if (rc < 0) {
     1148                  path.push_back(a);
     1149                  _next_out[tip] = a;
     1150                  if (path_arc[a]) {
     1151                    goto augment;   // a cycle is found, stop path search
     1152                  }
     1153                  tip = u;
     1154                  path_arc[a] = true;
     1155                  goto next_step;
     1156                }
     1157                else if (rc < min_red_cost) {
     1158                  min_red_cost = rc;
     1159                }
    11511160              }
    11521161            }
    11531162
    11541163            // Relabel tip node
    1155             min_red_cost = std::numeric_limits<LargeCost>::max();
    11561164            if (tip != start) {
    11571165              int ra = _reverse[path.back()];
    1158               min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
     1166              min_red_cost =
     1167                std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
    11591168            }
     1169            last_out = _next_out[tip];
    11601170            for (int a = _first_out[tip]; a != last_out; ++a) {
    1161               rc = _cost[a] + pi_tip - _pi[_target[a]];
    1162               if (_res_cap[a] > 0 && rc < min_red_cost) {
    1163                 min_red_cost = rc;
     1171              if (_res_cap[a] > 0) {
     1172                rc = _cost[a] + pi_tip - _pi[_target[a]];
     1173                if (rc < min_red_cost) {
     1174                  min_red_cost = rc;
     1175                }
    11641176              }
    11651177            }
    11661178            _pi[tip] -= min_red_cost + _epsilon;
     
    11691181
    11701182            // Step back
    11711183            if (tip != start) {
    1172               tip = _source[path.back()];
     1184              int pa = path.back();
     1185              path_arc[pa] = false;
     1186              tip = _source[pa];
    11731187              path.pop_back();
    11741188            }
    11751189
     
    11771191          }
    11781192
    11791193          // Augment along the found path (as much flow as possible)
     1194        augment:
    11801195          Value delta;
    11811196          int pa, u, v = start;
    11821197          for (int i = 0; i != int(path.size()); ++i) {
    11831198            pa = path[i];
    11841199            u = v;
    11851200            v = _target[pa];
     1201            path_arc[pa] = false;
    11861202            delta = std::min(_res_cap[pa], _excess[u]);
    11871203            _res_cap[pa] -= delta;
    11881204            _res_cap[_reverse[pa]] += delta;
    11891205            _excess[u] -= delta;
    11901206            _excess[v] += delta;
    1191             if (_excess[v] > 0 && _excess[v] <= delta)
     1207            if (_excess[v] > 0 && _excess[v] <= delta) {
    11921208              _active_nodes.push_back(v);
     1209            }
    11931210          }
     1211          path.clear();
    11941212
    11951213          // Global update heuristic
    1196           if (relabel_cnt >= next_update_limit) {
     1214          if (relabel_cnt >= next_global_update_limit) {
    11971215            globalUpdate();
    1198             next_update_limit += global_update_freq;
     1216            next_global_update_limit += global_update_skip;
    11991217          }
    12001218        }
     1219
    12011220      }
     1221
    12021222    }
    12031223
    12041224    /// Execute the algorithm performing push and relabel operations
     
    12071227      const int EARLY_TERM_EPSILON_LIMIT = 1000;
    12081228      const double GLOBAL_UPDATE_FACTOR = 2.0;
    12091229
    1210       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1230      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
    12111231        (_res_node_num + _sup_node_num * _sup_node_num));
    1212       int next_update_limit = global_update_freq;
    1213 
    1214       int relabel_cnt = 0;
     1232      int next_global_update_limit = global_update_skip;
    12151233
    12161234      // Perform cost scaling phases
    12171235      BoolVector hyper(_res_node_num, false);
    12181236      LargeCostVector hyper_cost(_res_node_num);
     1237      int relabel_cnt = 0;
    12191238      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    12201239                                        1 : _epsilon / _alpha )
    12211240      {
     
    12931312             min_red_cost = hyper[n] ? -hyper_cost[n] :
    12941313               std::numeric_limits<LargeCost>::max();
    12951314            for (int a = _first_out[n]; a != last_out; ++a) {
    1296               rc = _cost[a] + pi_n - _pi[_target[a]];
    1297               if (_res_cap[a] > 0 && rc < min_red_cost) {
    1298                 min_red_cost = rc;
     1315              if (_res_cap[a] > 0) {
     1316                rc = _cost[a] + pi_n - _pi[_target[a]];
     1317                if (rc < min_red_cost) {
     1318                  min_red_cost = rc;
     1319                }
    12991320              }
    13001321            }
    13011322            _pi[n] -= min_red_cost + _epsilon;
     
    13131334          }
    13141335
    13151336          // Global update heuristic
    1316           if (relabel_cnt >= next_update_limit) {
     1337          if (relabel_cnt >= next_global_update_limit) {
    13171338            globalUpdate();
    13181339            for (int u = 0; u != _res_node_num; ++u)
    13191340              hyper[u] = false;
    1320             next_update_limit += global_update_freq;
     1341            next_global_update_limit += global_update_skip;
    13211342          }
    13221343        }
    13231344      }
  • lemon/cost_scaling.h

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1300213941 -3600
    # Node ID f88191cb740f5469104cdde0e26d11d6a4edd3c8
    # Parent  4866b640dba97fae2fdfdd81cf4bec497d1cd48d
    Implement the scaling Price Refinement heuristic in CostScaling (#417)
    instead of Early Termination.
    
    These two heuristics are similar, but the newer one is faster
    and not only makes it possible to skip some epsilon phases, but
    it can improve the performance of the other phases, as well.
    
    diff --git a/lemon/cost_scaling.h b/lemon/cost_scaling.h
    a b  
    980980      }
    981981    }
    982982
    983     // Early termination heuristic
    984     bool earlyTermination() {
    985       const double EARLY_TERM_FACTOR = 3.0;
     983    // Price (potential) refinement heuristic
     984    bool priceRefinement() {
    986985
    987       // Build a static residual graph
    988       _arc_vec.clear();
    989       _cost_vec.clear();
    990       for (int j = 0; j != _res_arc_num; ++j) {
    991         if (_res_cap[j] > 0) {
    992           _arc_vec.push_back(IntPair(_source[j], _target[j]));
    993           _cost_vec.push_back(_cost[j] + 1);
     986      // Stack for stroing the topological order
     987      IntVector stack(_res_node_num);
     988      int stack_top;
     989
     990      // Perform phases
     991      while (topologicalSort(stack, stack_top)) {
     992
     993        // Compute node ranks in the acyclic admissible network and
     994        // store the nodes in buckets
     995        for (int i = 0; i != _res_node_num; ++i) {
     996          _rank[i] = 0;
    994997        }
     998        const int bucket_end = _root + 1;
     999        for (int r = 0; r != _max_rank; ++r) {
     1000          _buckets[r] = bucket_end;
     1001        }
     1002        int top_rank = 0;
     1003        for ( ; stack_top >= 0; --stack_top) {
     1004          int u = stack[stack_top], v;
     1005          int rank_u = _rank[u];
     1006
     1007          LargeCost rc, pi_u = _pi[u];
     1008          int last_out = _first_out[u+1];
     1009          for (int a = _first_out[u]; a != last_out; ++a) {
     1010            if (_res_cap[a] > 0) {
     1011              v = _target[a];
     1012              rc = _cost[a] + pi_u - _pi[v];
     1013              if (rc < 0) {
     1014                LargeCost nrc = static_cast<LargeCost>((-rc - 0.5) / _epsilon);
     1015                if (nrc < LargeCost(_max_rank)) {
     1016                  int new_rank_v = rank_u + static_cast<int>(nrc);
     1017                  if (new_rank_v > _rank[v]) {
     1018                    _rank[v] = new_rank_v;
     1019                  }
     1020                }
     1021              }
     1022            }
     1023          }
     1024
     1025          if (rank_u > 0) {
     1026            top_rank = std::max(top_rank, rank_u);
     1027            int bfirst = _buckets[rank_u];
     1028            _bucket_next[u] = bfirst;
     1029            _bucket_prev[bfirst] = u;
     1030            _buckets[rank_u] = u;
     1031          }
     1032        }
     1033
     1034        // Check if the current flow is epsilon-optimal
     1035        if (top_rank == 0) {
     1036          return true;
     1037        }
     1038
     1039        // Process buckets in top-down order
     1040        for (int rank = top_rank; rank > 0; --rank) {
     1041          while (_buckets[rank] != bucket_end) {
     1042            // Remove the first node from the current bucket
     1043            int u = _buckets[rank];
     1044            _buckets[rank] = _bucket_next[u];
     1045
     1046            // Search the outgoing arcs of u
     1047            LargeCost rc, pi_u = _pi[u];
     1048            int last_out = _first_out[u+1];
     1049            int v, old_rank_v, new_rank_v;
     1050            for (int a = _first_out[u]; a != last_out; ++a) {
     1051              if (_res_cap[a] > 0) {
     1052                v = _target[a];
     1053                old_rank_v = _rank[v];
     1054
     1055                if (old_rank_v < rank) {
     1056
     1057                  // Compute the new rank of node v
     1058                  rc = _cost[a] + pi_u - _pi[v];
     1059                  if (rc < 0) {
     1060                    new_rank_v = rank;
     1061                  } else {
     1062                    LargeCost nrc = rc / _epsilon;
     1063                    new_rank_v = 0;
     1064                    if (nrc < LargeCost(_max_rank)) {
     1065                      new_rank_v = rank - 1 - static_cast<int>(nrc);
     1066                    }
     1067                  }
     1068
     1069                  // Change the rank of node v
     1070                  if (new_rank_v > old_rank_v) {
     1071                    _rank[v] = new_rank_v;
     1072
     1073                    // Remove v from its old bucket
     1074                    if (old_rank_v > 0) {
     1075                      if (_buckets[old_rank_v] == v) {
     1076                        _buckets[old_rank_v] = _bucket_next[v];
     1077                      } else {
     1078                        int pv = _bucket_prev[v], nv = _bucket_next[v];
     1079                        _bucket_next[pv] = nv;
     1080                        _bucket_prev[nv] = pv;
     1081                      }
     1082                    }
     1083
     1084                    // Insert v into its new bucket
     1085                    int nv = _buckets[new_rank_v];
     1086                    _bucket_next[v] = nv;
     1087                    _bucket_prev[nv] = v;
     1088                    _buckets[new_rank_v] = v;
     1089                  }
     1090                }
     1091              }
     1092            }
     1093
     1094            // Refine potential of node u
     1095            _pi[u] -= rank * _epsilon;
     1096          }
     1097        }
     1098
    9951099      }
    996       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    9971100
    998       // Run Bellman-Ford algorithm to check if the current flow is optimal
    999       BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
    1000       bf.init(0);
    1001       bool done = false;
    1002       int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
    1003       for (int i = 0; i < K && !done; ++i) {
    1004         done = bf.processNextWeakRound();
     1101      return false;
     1102    }
     1103
     1104    // Find and cancel cycles in the admissible network and
     1105    // determine topological order using DFS
     1106    bool topologicalSort(IntVector &stack, int &stack_top) {
     1107      const int MAX_CYCLE_CANCEL = 1;
     1108
     1109      BoolVector reached(_res_node_num, false);
     1110      BoolVector processed(_res_node_num, false);
     1111      IntVector pred(_res_node_num);
     1112      for (int i = 0; i != _res_node_num; ++i) {
     1113        _next_out[i] = _first_out[i];
    10051114      }
    1006       return done;
     1115      stack_top = -1;
     1116
     1117      int cycle_cnt = 0;
     1118      for (int start = 0; start != _res_node_num; ++start) {
     1119        if (reached[start]) continue;
     1120
     1121        // Start DFS search from this start node
     1122        pred[start] = -1;
     1123        int tip = start, v;
     1124        while (true) {
     1125          // Check the outgoing arcs of the current tip node
     1126          reached[tip] = true;
     1127          LargeCost pi_tip = _pi[tip];
     1128          int a, last_out = _first_out[tip+1];
     1129          for (a = _next_out[tip]; a != last_out; ++a) {
     1130            if (_res_cap[a] > 0) {
     1131              v = _target[a];
     1132              if (_cost[a] + pi_tip - _pi[v] < 0) {
     1133                if (!reached[v]) {
     1134                  // A new node is reached
     1135                  reached[v] = true;
     1136                  pred[v] = tip;
     1137                  _next_out[tip] = a;
     1138                  tip = v;
     1139                  a = _next_out[tip];
     1140                  last_out = _first_out[tip+1];
     1141                  break;
     1142                }
     1143                else if (!processed[v]) {
     1144                  // A cycle is found
     1145                  ++cycle_cnt;
     1146                  _next_out[tip] = a;
     1147
     1148                  // Find the minimum residual capacity along the cycle
     1149                  Value d, delta = _res_cap[a];
     1150                  int u, delta_node = tip;
     1151                  for (u = tip; u != v; ) {
     1152                    u = pred[u];
     1153                    d = _res_cap[_next_out[u]];
     1154                    if (d <= delta) {
     1155                      delta = d;
     1156                      delta_node = u;
     1157                    }
     1158                  }
     1159
     1160                  // Augment along the cycle
     1161                  _res_cap[a] -= delta;
     1162                  _res_cap[_reverse[a]] += delta;
     1163                  for (u = tip; u != v; ) {
     1164                    u = pred[u];
     1165                    int ca = _next_out[u];
     1166                    _res_cap[ca] -= delta;
     1167                    _res_cap[_reverse[ca]] += delta;
     1168                  }
     1169
     1170                  // Check the maximum number of cycle canceling
     1171                  if (cycle_cnt >= MAX_CYCLE_CANCEL) {
     1172                    return false;
     1173                  }
     1174
     1175                  // Roll back search to delta_node
     1176                  if (delta_node != tip) {
     1177                    for (u = tip; u != delta_node; u = pred[u]) {
     1178                      reached[u] = false;
     1179                    }
     1180                    tip = delta_node;
     1181                    a = _next_out[tip] + 1;
     1182                    last_out = _first_out[tip+1];
     1183                    break;
     1184                  }
     1185                }
     1186              }
     1187            }
     1188          }
     1189
     1190          // Step back to the previous node
     1191          if (a == last_out) {
     1192            processed[tip] = true;
     1193            stack[++stack_top] = tip;
     1194            tip = pred[tip];
     1195            if (tip < 0) {
     1196              // Finish DFS from the current start node
     1197              break;
     1198            }
     1199            ++_next_out[tip];
     1200          }
     1201        }
     1202
     1203      }
     1204
     1205      return (cycle_cnt == 0);
    10071206    }
    10081207
    10091208    // Global potential update heuristic
     
    11021301    /// Execute the algorithm performing augment and relabel operations
    11031302    void startAugment(int max_length) {
    11041303      // Paramters for heuristics
    1105       const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1304      const int PRICE_REFINEMENT_LIMIT = 2;
    11061305      const double GLOBAL_UPDATE_FACTOR = 1.0;
    11071306      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
    11081307        (_res_node_num + _sup_node_num * _sup_node_num));
     
    11121311      IntVector path;
    11131312      BoolVector path_arc(_res_arc_num, false);
    11141313      int relabel_cnt = 0;
     1314      int eps_phase_cnt = 0;
    11151315      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    11161316                                        1 : _epsilon / _alpha )
    11171317      {
    1118         // Early termination heuristic
    1119         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
    1120           if (earlyTermination()) break;
     1318        ++eps_phase_cnt;
     1319
     1320        // Price refinement heuristic
     1321        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
     1322          if (priceRefinement()) continue;
    11211323        }
    11221324
    11231325        // Initialize current phase
     
    12241426    /// Execute the algorithm performing push and relabel operations
    12251427    void startPush() {
    12261428      // Paramters for heuristics
    1227       const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1429      const int PRICE_REFINEMENT_LIMIT = 2;
    12281430      const double GLOBAL_UPDATE_FACTOR = 2.0;
    12291431
    12301432      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
     
    12351437      BoolVector hyper(_res_node_num, false);
    12361438      LargeCostVector hyper_cost(_res_node_num);
    12371439      int relabel_cnt = 0;
     1440      int eps_phase_cnt = 0;
    12381441      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    12391442                                        1 : _epsilon / _alpha )
    12401443      {
    1241         // Early termination heuristic
    1242         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
    1243           if (earlyTermination()) break;
     1444        ++eps_phase_cnt;
     1445
     1446        // Price refinement heuristic
     1447        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
     1448          if (priceRefinement()) continue;
    12441449        }
    12451450
    12461451        // Initialize current phase
  • lemon/cost_scaling.h

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1300215151 -3600
    # Node ID 482ff194df9394723c9db6b69892cc34efa69f54
    # Parent  f88191cb740f5469104cdde0e26d11d6a4edd3c8
    Faster computation of the dual solution in CostScaling (#417)
    
    diff --git a/lemon/cost_scaling.h b/lemon/cost_scaling.h
    a b  
    237237      std::vector<Value>& _v;
    238238    };
    239239
    240     typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
    241240    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
    242241
    243242  private:
     
    288287    IntVector _rank;
    289288    int _max_rank;
    290289
    291     // Data for a StaticDigraph structure
    292     typedef std::pair<int, int> IntPair;
    293     StaticDigraph _sgr;
    294     std::vector<IntPair> _arc_vec;
    295     std::vector<LargeCost> _cost_vec;
    296     LargeCostArcMap _cost_map;
    297     LargeCostNodeMap _pi_map;
    298 
    299290  public:
    300291
    301292    /// \brief Constant for infinite upper bounds (capacities).
     
    342333    /// \param graph The digraph the algorithm runs on.
    343334    CostScaling(const GR& graph) :
    344335      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
    345       _cost_map(_cost_vec), _pi_map(_pi),
    346336      INF(std::numeric_limits<Value>::has_infinity ?
    347337          std::numeric_limits<Value>::infinity() :
    348338          std::numeric_limits<Value>::max())
     
    619609      _excess.resize(_res_node_num);
    620610      _next_out.resize(_res_node_num);
    621611
    622       _arc_vec.reserve(_res_arc_num);
    623       _cost_vec.reserve(_res_arc_num);
    624 
    625612      // Copy the graph
    626613      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
    627614      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     
    923910          break;
    924911      }
    925912
    926       // Compute node potentials for the original costs
    927       _arc_vec.clear();
    928       _cost_vec.clear();
    929       for (int j = 0; j != _res_arc_num; ++j) {
    930         if (_res_cap[j] > 0) {
    931           _arc_vec.push_back(IntPair(_source[j], _target[j]));
    932           _cost_vec.push_back(_scost[j]);
     913      // Compute node potentials (dual solution)
     914      for (int i = 0; i != _res_node_num; ++i) {
     915        _pi[i] = static_cast<Cost>(_pi[i] / (_res_node_num * _alpha));
     916      }
     917      bool optimal = true;
     918      for (int i = 0; optimal && i != _res_node_num; ++i) {
     919        LargeCost pi_i = _pi[i];
     920        int last_out = _first_out[i+1];
     921        for (int j = _first_out[i]; j != last_out; ++j) {
     922          if (_res_cap[j] > 0 && _scost[j] + pi_i - _pi[_target[j]] < 0) {
     923            optimal = false;
     924            break;
     925          }
    933926        }
    934927      }
    935       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    936928
    937       typename BellmanFord<StaticDigraph, LargeCostArcMap>
    938         ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
    939       bf.distMap(_pi_map);
    940       bf.init(0);
    941       bf.start();
     929      if (!optimal) {
     930        // Compute node potentials for the original costs with BellmanFord
     931        // (if it is necessary)
     932        typedef std::pair<int, int> IntPair;
     933        StaticDigraph sgr;
     934        std::vector<IntPair> arc_vec;
     935        std::vector<LargeCost> cost_vec;
     936        LargeCostArcMap cost_map(cost_vec);
     937
     938        arc_vec.clear();
     939        cost_vec.clear();
     940        for (int j = 0; j != _res_arc_num; ++j) {
     941          if (_res_cap[j] > 0) {
     942            int u = _source[j], v = _target[j];
     943            arc_vec.push_back(IntPair(u, v));
     944            cost_vec.push_back(_scost[j] + _pi[u] - _pi[v]);
     945          }
     946        }
     947        sgr.build(_res_node_num, arc_vec.begin(), arc_vec.end());
     948
     949        typename BellmanFord<StaticDigraph, LargeCostArcMap>::Create
     950          bf(sgr, cost_map);
     951        bf.init(0);
     952        bf.start();
     953
     954        for (int i = 0; i != _res_node_num; ++i) {
     955          _pi[i] += bf.dist(sgr.node(i));
     956        }
     957      }
     958
     959      // Shift potentials to meet the requirements of the GEQ type
     960      // optimality conditions
     961      LargeCost max_pot = _pi[_root];
     962      for (int i = 0; i != _res_node_num; ++i) {
     963        if (_pi[i] > max_pot) max_pot = _pi[i];
     964      }
     965      if (max_pot != 0) {
     966        for (int i = 0; i != _res_node_num; ++i) {
     967          _pi[i] -= max_pot;
     968        }
     969      }
    942970
    943971      // Handle non-zero lower bounds
    944972      if (_have_lower) {
  • lemon/cost_scaling.h

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1300215251 -3600
    # Node ID 78d9eddfaf1add430990966beaea64a47eb2d7bf
    # Parent  482ff194df9394723c9db6b69892cc34efa69f54
    Change the default scaling factor in CostScaling (#417)
    
    diff --git a/lemon/cost_scaling.h b/lemon/cost_scaling.h
    a b  
    487487    ///
    488488    /// \param method The internal method that will be used in the
    489489    /// algorithm. For more information, see \ref Method.
    490     /// \param factor The cost scaling factor. It must be larger than one.
     490    /// \param factor The cost scaling factor. It must be at least two.
    491491    ///
    492492    /// \return \c INFEASIBLE if no feasible flow exists,
    493493    /// \n \c OPTIMAL if the problem has optimal solution
     
    501501    ///
    502502    /// \see ProblemType, Method
    503503    /// \see resetParams(), reset()
    504     ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
     504    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 16) {
     505      LEMON_ASSERT(factor >= 2, "The scaling factor must be at least 2");
    505506      _alpha = factor;
    506507      ProblemType pt = init();
    507508      if (pt != OPTIMAL) return pt;
  • lemon/cost_scaling.h

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1300208397 -3600
    # Node ID ca9d1dc770268c0b476d33857ece1c13471ea447
    # Parent  8e39ccaabf484a72dc36608cb5e954a62a0b984d
    Minor improvements in CostScaling (#417)
    
    diff --git a/lemon/cost_scaling.h b/lemon/cost_scaling.h
    a b  
    575575      return *this;
    576576    }
    577577
    578     /// \brief Reset all the parameters that have been given before.
     578    /// \brief Reset the internal data structures and all the parameters
     579    /// that have been given before.
    579580    ///
    580     /// This function resets all the paramaters that have been given
    581     /// before using functions \ref lowerMap(), \ref upperMap(),
    582     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
     581    /// This function resets the internal data structures and all the
     582    /// paramaters that have been given before using functions \ref lowerMap(),
     583    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
    583584    ///
    584     /// It is useful for multiple run() calls. If this function is not
    585     /// used, all the parameters given before are kept for the next
    586     /// \ref run() call.
    587     /// However, the underlying digraph must not be modified after this
    588     /// class have been constructed, since it copies and extends the graph.
     585    /// It is useful for multiple \ref run() calls. By default, all the given
     586    /// parameters are kept for the next \ref run() call, unless
     587    /// \ref resetParams() or \ref reset() is used.
     588    /// If the underlying digraph was also modified after the construction
     589    /// of the class or the last \ref reset() call, then the \ref reset()
     590    /// function must be used, otherwise \ref resetParams() is sufficient.
     591    ///
     592    /// See \ref resetParams() for examples.
     593    ///
    589594    /// \return <tt>(*this)</tt>
     595    ///
     596    /// \see resetParams(), run()
    590597    CostScaling& reset() {
    591598      // Resize vectors
    592599      _node_num = countNodes(_graph);
     
    890897        }
    891898      }
    892899
    893       return OPTIMAL;
    894     }
    895 
    896     // Execute the algorithm and transform the results
    897     void start(Method method) {
    898       // Maximum path length for partial augment
    899       const int MAX_PATH_LENGTH = 4;
    900 
    901900      // Initialize data structures for buckets
    902901      _max_rank = _alpha * _res_node_num;
    903902      _buckets.resize(_max_rank);
     
    905904      _bucket_prev.resize(_res_node_num + 1);
    906905      _rank.resize(_res_node_num + 1);
    907906
    908       // Execute the algorithm
     907      return OPTIMAL;
     908    }
     909
     910    // Execute the algorithm and transform the results
     911    void start(Method method) {
     912      const int MAX_PARTIAL_PATH_LENGTH = 4;
     913
    909914      switch (method) {
    910915        case PUSH:
    911916          startPush();
    912917          break;
    913918        case AUGMENT:
    914           startAugment();
     919          startAugment(_res_node_num - 1);
    915920          break;
    916921        case PARTIAL_AUGMENT:
    917           startAugment(MAX_PATH_LENGTH);
     922          startAugment(MAX_PARTIAL_PATH_LENGTH);
    918923          break;
    919924      }
    920925
     
    951956        int last_out = _first_out[u+1];
    952957        LargeCost pi_u = _pi[u];
    953958        for (int a = _first_out[u]; a != last_out; ++a) {
    954           int v = _target[a];
    955           if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
    956             Value delta = _res_cap[a];
    957             _excess[u] -= delta;
    958             _excess[v] += delta;
    959             _res_cap[a] = 0;
    960             _res_cap[_reverse[a]] += delta;
     959          Value delta = _res_cap[a];
     960          if (delta > 0) {
     961            int v = _target[a];
     962            if (_cost[a] + pi_u - _pi[v] < 0) {
     963              _excess[u] -= delta;
     964              _excess[v] += delta;
     965              _res_cap[a] = 0;
     966              _res_cap[_reverse[a]] += delta;
     967            }
    961968          }
    962969        }
    963970      }
     
    10011008
    10021009    // Global potential update heuristic
    10031010    void globalUpdate() {
    1004       int bucket_end = _root + 1;
     1011      const int bucket_end = _root + 1;
    10051012
    10061013      // Initialize buckets
    10071014      for (int r = 0; r != _max_rank; ++r) {
    10081015        _buckets[r] = bucket_end;
    10091016      }
    10101017      Value total_excess = 0;
     1018      int b0 = bucket_end;
    10111019      for (int i = 0; i != _res_node_num; ++i) {
    10121020        if (_excess[i] < 0) {
    10131021          _rank[i] = 0;
    1014           _bucket_next[i] = _buckets[0];
    1015           _bucket_prev[_buckets[0]] = i;
    1016           _buckets[0] = i;
     1022          _bucket_next[i] = b0;
     1023          _bucket_prev[b0] = i;
     1024          b0 = i;
    10171025        } else {
    10181026          total_excess += _excess[i];
    10191027          _rank[i] = _max_rank;
    10201028        }
    10211029      }
    10221030      if (total_excess == 0) return;
     1031      _buckets[0] = b0;
    10231032
    10241033      // Search the buckets
    10251034      int r = 0;
     
    10411050                // Compute the new rank of v
    10421051                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
    10431052                int new_rank_v = old_rank_v;
    1044                 if (nrc < LargeCost(_max_rank))
    1045                   new_rank_v = r + 1 + int(nrc);
     1053                if (nrc < LargeCost(_max_rank)) {
     1054                  new_rank_v = r + 1 + static_cast<int>(nrc);
     1055                }
    10461056
    10471057                // Change the rank of v
    10481058                if (new_rank_v < old_rank_v) {
     
    10541064                    if (_buckets[old_rank_v] == v) {
    10551065                      _buckets[old_rank_v] = _bucket_next[v];
    10561066                    } else {
    1057                       _bucket_next[_bucket_prev[v]] = _bucket_next[v];
    1058                       _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
     1067                      int pv = _bucket_prev[v], nv = _bucket_next[v];
     1068                      _bucket_next[pv] = nv;
     1069                      _bucket_prev[nv] = pv;
    10591070                    }
    10601071                  }
    10611072
    1062                   // Insert v to its new bucket
    1063                   _bucket_next[v] = _buckets[new_rank_v];
    1064                   _bucket_prev[_buckets[new_rank_v]] = v;
     1073                  // Insert v into its new bucket
     1074                  int nv = _buckets[new_rank_v];
     1075                  _bucket_next[v] = nv;
     1076                  _bucket_prev[nv] = v;
    10651077                  _buckets[new_rank_v] = v;
    10661078                }
    10671079              }