COIN-OR::LEMON - Graph Library

Changeset 1047:ddd3c0d3d9bf in lemon for lemon/cost_scaling.h


Ignore:
Timestamp:
03/15/11 19:32:21 (9 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Message:

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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/cost_scaling.h

    r1046 r1047  
    981981    }
    982982
    983     // Early termination heuristic
    984     bool earlyTermination() {
    985       const double EARLY_TERM_FACTOR = 3.0;
    986 
    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);
    994         }
    995       }
    996       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    997 
    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();
    1005       }
    1006       return done;
     983    // Price (potential) refinement heuristic
     984    bool priceRefinement() {
     985
     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;
     997        }
     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
     1099      }
     1100
     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];
     1114      }
     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
     
    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 *
     
    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
     
    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
     
    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
Note: See TracChangeset for help on using the changeset viewer.