COIN-OR::LEMON - Graph Library

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/cost_scaling.h

    r1041 r1049  
    9898  /// "preflow push-relabel" algorithm for the maximum flow problem.
    9999  ///
     100  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
     101  /// implementations available in LEMON for this problem.
     102  ///
    100103  /// Most of the parameters of the problem (except for the digraph)
    101104  /// can be given using separate functions, and the algorithm can be
     
    114117  /// consider to use the named template parameters instead.
    115118  ///
    116   /// \warning Both number types must be signed and all input data must
     119  /// \warning Both \c V and \c C must be signed number types.
     120  /// \warning All input data (capacities, supply values, and costs) must
    117121  /// be integer.
    118   /// \warning This algorithm does not support negative costs for such
    119   /// arcs that have infinite upper bound.
     122  /// \warning This algorithm does not support negative costs for
     123  /// arcs having infinite upper bound.
    120124  ///
    121125  /// \note %CostScaling provides three different internal methods,
     
    179183    /// relabel operation.
    180184    /// By default, the so called \ref PARTIAL_AUGMENT
    181     /// "Partial Augment-Relabel" method is used, which proved to be
     185    /// "Partial Augment-Relabel" method is used, which turned out to be
    182186    /// the most efficient and the most robust on various test inputs.
    183187    /// However, the other methods can be selected using the \ref run()
     
    234238    };
    235239
    236     typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
    237240    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
    238241
     
    285288    int _max_rank;
    286289
    287     // Data for a StaticDigraph structure
    288     typedef std::pair<int, int> IntPair;
    289     StaticDigraph _sgr;
    290     std::vector<IntPair> _arc_vec;
    291     std::vector<LargeCost> _cost_vec;
    292     LargeCostArcMap _cost_map;
    293     LargeCostNodeMap _pi_map;
    294 
    295290  public:
    296291
     
    339334    CostScaling(const GR& graph) :
    340335      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
    341       _cost_map(_cost_vec), _pi_map(_pi),
    342336      INF(std::numeric_limits<Value>::has_infinity ?
    343337          std::numeric_limits<Value>::infinity() :
     
    448442    ///
    449443    /// Using this function has the same effect as using \ref supplyMap()
    450     /// with such a map in which \c k is assigned to \c s, \c -k is
     444    /// with a map in which \c k is assigned to \c s, \c -k is
    451445    /// assigned to \c t and all other nodes have zero supply value.
    452446    ///
     
    494488    /// \param method The internal method that will be used in the
    495489    /// algorithm. For more information, see \ref Method.
    496     /// \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.
    497491    ///
    498492    /// \return \c INFEASIBLE if no feasible flow exists,
     
    508502    /// \see ProblemType, Method
    509503    /// \see resetParams(), reset()
    510     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");
    511506      _alpha = factor;
    512507      ProblemType pt = init();
     
    572567    }
    573568
    574     /// \brief Reset all the parameters that have been given before.
    575     ///
    576     /// This function resets all the paramaters that have been given
    577     /// before using functions \ref lowerMap(), \ref upperMap(),
    578     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
    579     ///
    580     /// It is useful for multiple run() calls. If this function is not
    581     /// used, all the parameters given before are kept for the next
    582     /// \ref run() call.
    583     /// However, the underlying digraph must not be modified after this
    584     /// class have been constructed, since it copies and extends the graph.
     569    /// \brief Reset the internal data structures and all the parameters
     570    /// that have been given before.
     571    ///
     572    /// This function resets the internal data structures and all the
     573    /// paramaters that have been given before using functions \ref lowerMap(),
     574    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     575    ///
     576    /// It is useful for multiple \ref run() calls. By default, all the given
     577    /// parameters are kept for the next \ref run() call, unless
     578    /// \ref resetParams() or \ref reset() is used.
     579    /// If the underlying digraph was also modified after the construction
     580    /// of the class or the last \ref reset() call, then the \ref reset()
     581    /// function must be used, otherwise \ref resetParams() is sufficient.
     582    ///
     583    /// See \ref resetParams() for examples.
     584    ///
    585585    /// \return <tt>(*this)</tt>
     586    ///
     587    /// \see resetParams(), run()
    586588    CostScaling& reset() {
    587589      // Resize vectors
     
    608610      _excess.resize(_res_node_num);
    609611      _next_out.resize(_res_node_num);
    610 
    611       _arc_vec.reserve(_res_arc_num);
    612       _cost_vec.reserve(_res_arc_num);
    613612
    614613      // Copy the graph
     
    887886      }
    888887
    889       return OPTIMAL;
    890     }
    891 
    892     // Execute the algorithm and transform the results
    893     void start(Method method) {
    894       // Maximum path length for partial augment
    895       const int MAX_PATH_LENGTH = 4;
    896 
    897888      // Initialize data structures for buckets
    898889      _max_rank = _alpha * _res_node_num;
     
    902893      _rank.resize(_res_node_num + 1);
    903894
    904       // Execute the algorithm
     895      return OPTIMAL;
     896    }
     897
     898    // Execute the algorithm and transform the results
     899    void start(Method method) {
     900      const int MAX_PARTIAL_PATH_LENGTH = 4;
     901
    905902      switch (method) {
    906903        case PUSH:
     
    911908          break;
    912909        case PARTIAL_AUGMENT:
    913           startAugment(MAX_PATH_LENGTH);
     910          startAugment(MAX_PARTIAL_PATH_LENGTH);
    914911          break;
    915912      }
    916913
    917       // Compute node potentials for the original costs
    918       _arc_vec.clear();
    919       _cost_vec.clear();
    920       for (int j = 0; j != _res_arc_num; ++j) {
    921         if (_res_cap[j] > 0) {
    922           _arc_vec.push_back(IntPair(_source[j], _target[j]));
    923           _cost_vec.push_back(_scost[j]);
    924         }
    925       }
    926       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    927 
    928       typename BellmanFord<StaticDigraph, LargeCostArcMap>
    929         ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
    930       bf.distMap(_pi_map);
    931       bf.init(0);
    932       bf.start();
     914      // Compute node potentials (dual solution)
     915      for (int i = 0; i != _res_node_num; ++i) {
     916        _pi[i] = static_cast<Cost>(_pi[i] / (_res_node_num * _alpha));
     917      }
     918      bool optimal = true;
     919      for (int i = 0; optimal && i != _res_node_num; ++i) {
     920        LargeCost pi_i = _pi[i];
     921        int last_out = _first_out[i+1];
     922        for (int j = _first_out[i]; j != last_out; ++j) {
     923          if (_res_cap[j] > 0 && _scost[j] + pi_i - _pi[_target[j]] < 0) {
     924            optimal = false;
     925            break;
     926          }
     927        }
     928      }
     929
     930      if (!optimal) {
     931        // Compute node potentials for the original costs with BellmanFord
     932        // (if it is necessary)
     933        typedef std::pair<int, int> IntPair;
     934        StaticDigraph sgr;
     935        std::vector<IntPair> arc_vec;
     936        std::vector<LargeCost> cost_vec;
     937        LargeCostArcMap cost_map(cost_vec);
     938
     939        arc_vec.clear();
     940        cost_vec.clear();
     941        for (int j = 0; j != _res_arc_num; ++j) {
     942          if (_res_cap[j] > 0) {
     943            int u = _source[j], v = _target[j];
     944            arc_vec.push_back(IntPair(u, v));
     945            cost_vec.push_back(_scost[j] + _pi[u] - _pi[v]);
     946          }
     947        }
     948        sgr.build(_res_node_num, arc_vec.begin(), arc_vec.end());
     949
     950        typename BellmanFord<StaticDigraph, LargeCostArcMap>::Create
     951          bf(sgr, cost_map);
     952        bf.init(0);
     953        bf.start();
     954
     955        for (int i = 0; i != _res_node_num; ++i) {
     956          _pi[i] += bf.dist(sgr.node(i));
     957        }
     958      }
     959
     960      // Shift potentials to meet the requirements of the GEQ type
     961      // optimality conditions
     962      LargeCost max_pot = _pi[_root];
     963      for (int i = 0; i != _res_node_num; ++i) {
     964        if (_pi[i] > max_pot) max_pot = _pi[i];
     965      }
     966      if (max_pot != 0) {
     967        for (int i = 0; i != _res_node_num; ++i) {
     968          _pi[i] -= max_pot;
     969        }
     970      }
    933971
    934972      // Handle non-zero lower bounds
     
    948986        LargeCost pi_u = _pi[u];
    949987        for (int a = _first_out[u]; a != last_out; ++a) {
    950           int v = _target[a];
    951           if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
    952             Value delta = _res_cap[a];
    953             _excess[u] -= delta;
    954             _excess[v] += delta;
    955             _res_cap[a] = 0;
    956             _res_cap[_reverse[a]] += delta;
     988          Value delta = _res_cap[a];
     989          if (delta > 0) {
     990            int v = _target[a];
     991            if (_cost[a] + pi_u - _pi[v] < 0) {
     992              _excess[u] -= delta;
     993              _excess[v] += delta;
     994              _res_cap[a] = 0;
     995              _res_cap[_reverse[a]] += delta;
     996            }
    957997          }
    958998        }
     
    9701010    }
    9711011
    972     // Early termination heuristic
    973     bool earlyTermination() {
    974       const double EARLY_TERM_FACTOR = 3.0;
    975 
    976       // Build a static residual graph
    977       _arc_vec.clear();
    978       _cost_vec.clear();
    979       for (int j = 0; j != _res_arc_num; ++j) {
    980         if (_res_cap[j] > 0) {
    981           _arc_vec.push_back(IntPair(_source[j], _target[j]));
    982           _cost_vec.push_back(_cost[j] + 1);
    983         }
    984       }
    985       _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    986 
    987       // Run Bellman-Ford algorithm to check if the current flow is optimal
    988       BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
    989       bf.init(0);
    990       bool done = false;
    991       int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
    992       for (int i = 0; i < K && !done; ++i) {
    993         done = bf.processNextWeakRound();
    994       }
    995       return done;
     1012    // Price (potential) refinement heuristic
     1013    bool priceRefinement() {
     1014
     1015      // Stack for stroing the topological order
     1016      IntVector stack(_res_node_num);
     1017      int stack_top;
     1018
     1019      // Perform phases
     1020      while (topologicalSort(stack, stack_top)) {
     1021
     1022        // Compute node ranks in the acyclic admissible network and
     1023        // store the nodes in buckets
     1024        for (int i = 0; i != _res_node_num; ++i) {
     1025          _rank[i] = 0;
     1026        }
     1027        const int bucket_end = _root + 1;
     1028        for (int r = 0; r != _max_rank; ++r) {
     1029          _buckets[r] = bucket_end;
     1030        }
     1031        int top_rank = 0;
     1032        for ( ; stack_top >= 0; --stack_top) {
     1033          int u = stack[stack_top], v;
     1034          int rank_u = _rank[u];
     1035
     1036          LargeCost rc, pi_u = _pi[u];
     1037          int last_out = _first_out[u+1];
     1038          for (int a = _first_out[u]; a != last_out; ++a) {
     1039            if (_res_cap[a] > 0) {
     1040              v = _target[a];
     1041              rc = _cost[a] + pi_u - _pi[v];
     1042              if (rc < 0) {
     1043                LargeCost nrc = static_cast<LargeCost>((-rc - 0.5) / _epsilon);
     1044                if (nrc < LargeCost(_max_rank)) {
     1045                  int new_rank_v = rank_u + static_cast<int>(nrc);
     1046                  if (new_rank_v > _rank[v]) {
     1047                    _rank[v] = new_rank_v;
     1048                  }
     1049                }
     1050              }
     1051            }
     1052          }
     1053
     1054          if (rank_u > 0) {
     1055            top_rank = std::max(top_rank, rank_u);
     1056            int bfirst = _buckets[rank_u];
     1057            _bucket_next[u] = bfirst;
     1058            _bucket_prev[bfirst] = u;
     1059            _buckets[rank_u] = u;
     1060          }
     1061        }
     1062
     1063        // Check if the current flow is epsilon-optimal
     1064        if (top_rank == 0) {
     1065          return true;
     1066        }
     1067
     1068        // Process buckets in top-down order
     1069        for (int rank = top_rank; rank > 0; --rank) {
     1070          while (_buckets[rank] != bucket_end) {
     1071            // Remove the first node from the current bucket
     1072            int u = _buckets[rank];
     1073            _buckets[rank] = _bucket_next[u];
     1074
     1075            // Search the outgoing arcs of u
     1076            LargeCost rc, pi_u = _pi[u];
     1077            int last_out = _first_out[u+1];
     1078            int v, old_rank_v, new_rank_v;
     1079            for (int a = _first_out[u]; a != last_out; ++a) {
     1080              if (_res_cap[a] > 0) {
     1081                v = _target[a];
     1082                old_rank_v = _rank[v];
     1083
     1084                if (old_rank_v < rank) {
     1085
     1086                  // Compute the new rank of node v
     1087                  rc = _cost[a] + pi_u - _pi[v];
     1088                  if (rc < 0) {
     1089                    new_rank_v = rank;
     1090                  } else {
     1091                    LargeCost nrc = rc / _epsilon;
     1092                    new_rank_v = 0;
     1093                    if (nrc < LargeCost(_max_rank)) {
     1094                      new_rank_v = rank - 1 - static_cast<int>(nrc);
     1095                    }
     1096                  }
     1097
     1098                  // Change the rank of node v
     1099                  if (new_rank_v > old_rank_v) {
     1100                    _rank[v] = new_rank_v;
     1101
     1102                    // Remove v from its old bucket
     1103                    if (old_rank_v > 0) {
     1104                      if (_buckets[old_rank_v] == v) {
     1105                        _buckets[old_rank_v] = _bucket_next[v];
     1106                      } else {
     1107                        int pv = _bucket_prev[v], nv = _bucket_next[v];
     1108                        _bucket_next[pv] = nv;
     1109                        _bucket_prev[nv] = pv;
     1110                      }
     1111                    }
     1112
     1113                    // Insert v into its new bucket
     1114                    int nv = _buckets[new_rank_v];
     1115                    _bucket_next[v] = nv;
     1116                    _bucket_prev[nv] = v;
     1117                    _buckets[new_rank_v] = v;
     1118                  }
     1119                }
     1120              }
     1121            }
     1122
     1123            // Refine potential of node u
     1124            _pi[u] -= rank * _epsilon;
     1125          }
     1126        }
     1127
     1128      }
     1129
     1130      return false;
     1131    }
     1132
     1133    // Find and cancel cycles in the admissible network and
     1134    // determine topological order using DFS
     1135    bool topologicalSort(IntVector &stack, int &stack_top) {
     1136      const int MAX_CYCLE_CANCEL = 1;
     1137
     1138      BoolVector reached(_res_node_num, false);
     1139      BoolVector processed(_res_node_num, false);
     1140      IntVector pred(_res_node_num);
     1141      for (int i = 0; i != _res_node_num; ++i) {
     1142        _next_out[i] = _first_out[i];
     1143      }
     1144      stack_top = -1;
     1145
     1146      int cycle_cnt = 0;
     1147      for (int start = 0; start != _res_node_num; ++start) {
     1148        if (reached[start]) continue;
     1149
     1150        // Start DFS search from this start node
     1151        pred[start] = -1;
     1152        int tip = start, v;
     1153        while (true) {
     1154          // Check the outgoing arcs of the current tip node
     1155          reached[tip] = true;
     1156          LargeCost pi_tip = _pi[tip];
     1157          int a, last_out = _first_out[tip+1];
     1158          for (a = _next_out[tip]; a != last_out; ++a) {
     1159            if (_res_cap[a] > 0) {
     1160              v = _target[a];
     1161              if (_cost[a] + pi_tip - _pi[v] < 0) {
     1162                if (!reached[v]) {
     1163                  // A new node is reached
     1164                  reached[v] = true;
     1165                  pred[v] = tip;
     1166                  _next_out[tip] = a;
     1167                  tip = v;
     1168                  a = _next_out[tip];
     1169                  last_out = _first_out[tip+1];
     1170                  break;
     1171                }
     1172                else if (!processed[v]) {
     1173                  // A cycle is found
     1174                  ++cycle_cnt;
     1175                  _next_out[tip] = a;
     1176
     1177                  // Find the minimum residual capacity along the cycle
     1178                  Value d, delta = _res_cap[a];
     1179                  int u, delta_node = tip;
     1180                  for (u = tip; u != v; ) {
     1181                    u = pred[u];
     1182                    d = _res_cap[_next_out[u]];
     1183                    if (d <= delta) {
     1184                      delta = d;
     1185                      delta_node = u;
     1186                    }
     1187                  }
     1188
     1189                  // Augment along the cycle
     1190                  _res_cap[a] -= delta;
     1191                  _res_cap[_reverse[a]] += delta;
     1192                  for (u = tip; u != v; ) {
     1193                    u = pred[u];
     1194                    int ca = _next_out[u];
     1195                    _res_cap[ca] -= delta;
     1196                    _res_cap[_reverse[ca]] += delta;
     1197                  }
     1198
     1199                  // Check the maximum number of cycle canceling
     1200                  if (cycle_cnt >= MAX_CYCLE_CANCEL) {
     1201                    return false;
     1202                  }
     1203
     1204                  // Roll back search to delta_node
     1205                  if (delta_node != tip) {
     1206                    for (u = tip; u != delta_node; u = pred[u]) {
     1207                      reached[u] = false;
     1208                    }
     1209                    tip = delta_node;
     1210                    a = _next_out[tip] + 1;
     1211                    last_out = _first_out[tip+1];
     1212                    break;
     1213                  }
     1214                }
     1215              }
     1216            }
     1217          }
     1218
     1219          // Step back to the previous node
     1220          if (a == last_out) {
     1221            processed[tip] = true;
     1222            stack[++stack_top] = tip;
     1223            tip = pred[tip];
     1224            if (tip < 0) {
     1225              // Finish DFS from the current start node
     1226              break;
     1227            }
     1228            ++_next_out[tip];
     1229          }
     1230        }
     1231
     1232      }
     1233
     1234      return (cycle_cnt == 0);
    9961235    }
    9971236
    9981237    // Global potential update heuristic
    9991238    void globalUpdate() {
    1000       int bucket_end = _root + 1;
     1239      const int bucket_end = _root + 1;
    10011240
    10021241      // Initialize buckets
     
    10051244      }
    10061245      Value total_excess = 0;
     1246      int b0 = bucket_end;
    10071247      for (int i = 0; i != _res_node_num; ++i) {
    10081248        if (_excess[i] < 0) {
    10091249          _rank[i] = 0;
    1010           _bucket_next[i] = _buckets[0];
    1011           _bucket_prev[_buckets[0]] = i;
    1012           _buckets[0] = i;
     1250          _bucket_next[i] = b0;
     1251          _bucket_prev[b0] = i;
     1252          b0 = i;
    10131253        } else {
    10141254          total_excess += _excess[i];
     
    10171257      }
    10181258      if (total_excess == 0) return;
     1259      _buckets[0] = b0;
    10191260
    10201261      // Search the buckets
     
    10381279                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
    10391280                int new_rank_v = old_rank_v;
    1040                 if (nrc < LargeCost(_max_rank))
    1041                   new_rank_v = r + 1 + int(nrc);
     1281                if (nrc < LargeCost(_max_rank)) {
     1282                  new_rank_v = r + 1 + static_cast<int>(nrc);
     1283                }
    10421284
    10431285                // Change the rank of v
     
    10511293                      _buckets[old_rank_v] = _bucket_next[v];
    10521294                    } else {
    1053                       _bucket_next[_bucket_prev[v]] = _bucket_next[v];
    1054                       _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
     1295                      int pv = _bucket_prev[v], nv = _bucket_next[v];
     1296                      _bucket_next[pv] = nv;
     1297                      _bucket_prev[nv] = pv;
    10551298                    }
    10561299                  }
    10571300
    1058                   // Insert v to its new bucket
    1059                   _bucket_next[v] = _buckets[new_rank_v];
    1060                   _bucket_prev[_buckets[new_rank_v]] = v;
     1301                  // Insert v into its new bucket
     1302                  int nv = _buckets[new_rank_v];
     1303                  _bucket_next[v] = nv;
     1304                  _bucket_prev[nv] = v;
    10611305                  _buckets[new_rank_v] = v;
    10621306                }
     
    10871331    void startAugment(int max_length) {
    10881332      // Paramters for heuristics
    1089       const int EARLY_TERM_EPSILON_LIMIT = 1000;
    1090       const double GLOBAL_UPDATE_FACTOR = 3.0;
    1091 
    1092       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1333      const int PRICE_REFINEMENT_LIMIT = 2;
     1334      const double GLOBAL_UPDATE_FACTOR = 1.0;
     1335      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
    10931336        (_res_node_num + _sup_node_num * _sup_node_num));
    1094       int next_update_limit = global_update_freq;
    1095 
     1337      int next_global_update_limit = global_update_skip;
     1338
     1339      // Perform cost scaling phases
     1340      IntVector path;
     1341      BoolVector path_arc(_res_arc_num, false);
    10961342      int relabel_cnt = 0;
    1097 
    1098       // Perform cost scaling phases
    1099       std::vector<int> path;
     1343      int eps_phase_cnt = 0;
    11001344      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    11011345                                        1 : _epsilon / _alpha )
    11021346      {
    1103         // Early termination heuristic
    1104         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
    1105           if (earlyTermination()) break;
     1347        ++eps_phase_cnt;
     1348
     1349        // Price refinement heuristic
     1350        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
     1351          if (priceRefinement()) continue;
    11061352        }
    11071353
     
    11201366
    11211367          // Find an augmenting path from the start node
    1122           path.clear();
    11231368          int tip = start;
    1124           while (_excess[tip] >= 0 && int(path.size()) < max_length) {
     1369          while (int(path.size()) < max_length && _excess[tip] >= 0) {
    11251370            int u;
    1126             LargeCost min_red_cost, rc, pi_tip = _pi[tip];
     1371            LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
     1372            LargeCost pi_tip = _pi[tip];
    11271373            int last_out = _first_out[tip+1];
    11281374            for (int a = _next_out[tip]; a != last_out; ++a) {
    1129               u = _target[a];
    1130               if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
    1131                 path.push_back(a);
    1132                 _next_out[tip] = a;
    1133                 tip = u;
    1134                 goto next_step;
     1375              if (_res_cap[a] > 0) {
     1376                u = _target[a];
     1377                rc = _cost[a] + pi_tip - _pi[u];
     1378                if (rc < 0) {
     1379                  path.push_back(a);
     1380                  _next_out[tip] = a;
     1381                  if (path_arc[a]) {
     1382                    goto augment;   // a cycle is found, stop path search
     1383                  }
     1384                  tip = u;
     1385                  path_arc[a] = true;
     1386                  goto next_step;
     1387                }
     1388                else if (rc < min_red_cost) {
     1389                  min_red_cost = rc;
     1390                }
    11351391              }
    11361392            }
    11371393
    11381394            // Relabel tip node
    1139             min_red_cost = std::numeric_limits<LargeCost>::max();
    11401395            if (tip != start) {
    11411396              int ra = _reverse[path.back()];
    1142               min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
     1397              min_red_cost =
     1398                std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
    11431399            }
     1400            last_out = _next_out[tip];
    11441401            for (int a = _first_out[tip]; a != last_out; ++a) {
    1145               rc = _cost[a] + pi_tip - _pi[_target[a]];
    1146               if (_res_cap[a] > 0 && rc < min_red_cost) {
    1147                 min_red_cost = rc;
     1402              if (_res_cap[a] > 0) {
     1403                rc = _cost[a] + pi_tip - _pi[_target[a]];
     1404                if (rc < min_red_cost) {
     1405                  min_red_cost = rc;
     1406                }
    11481407              }
    11491408            }
     
    11541413            // Step back
    11551414            if (tip != start) {
    1156               tip = _source[path.back()];
     1415              int pa = path.back();
     1416              path_arc[pa] = false;
     1417              tip = _source[pa];
    11571418              path.pop_back();
    11581419            }
     
    11621423
    11631424          // Augment along the found path (as much flow as possible)
     1425        augment:
    11641426          Value delta;
    11651427          int pa, u, v = start;
     
    11681430            u = v;
    11691431            v = _target[pa];
     1432            path_arc[pa] = false;
    11701433            delta = std::min(_res_cap[pa], _excess[u]);
    11711434            _res_cap[pa] -= delta;
     
    11731436            _excess[u] -= delta;
    11741437            _excess[v] += delta;
    1175             if (_excess[v] > 0 && _excess[v] <= delta)
     1438            if (_excess[v] > 0 && _excess[v] <= delta) {
    11761439              _active_nodes.push_back(v);
    1177           }
     1440            }
     1441          }
     1442          path.clear();
    11781443
    11791444          // Global update heuristic
    1180           if (relabel_cnt >= next_update_limit) {
     1445          if (relabel_cnt >= next_global_update_limit) {
    11811446            globalUpdate();
    1182             next_update_limit += global_update_freq;
    1183           }
    1184         }
    1185       }
     1447            next_global_update_limit += global_update_skip;
     1448          }
     1449        }
     1450
     1451      }
     1452
    11861453    }
    11871454
     
    11891456    void startPush() {
    11901457      // Paramters for heuristics
    1191       const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1458      const int PRICE_REFINEMENT_LIMIT = 2;
    11921459      const double GLOBAL_UPDATE_FACTOR = 2.0;
    11931460
    1194       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1461      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
    11951462        (_res_node_num + _sup_node_num * _sup_node_num));
    1196       int next_update_limit = global_update_freq;
    1197 
    1198       int relabel_cnt = 0;
     1463      int next_global_update_limit = global_update_skip;
    11991464
    12001465      // Perform cost scaling phases
    12011466      BoolVector hyper(_res_node_num, false);
    12021467      LargeCostVector hyper_cost(_res_node_num);
     1468      int relabel_cnt = 0;
     1469      int eps_phase_cnt = 0;
    12031470      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    12041471                                        1 : _epsilon / _alpha )
    12051472      {
    1206         // Early termination heuristic
    1207         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
    1208           if (earlyTermination()) break;
     1473        ++eps_phase_cnt;
     1474
     1475        // Price refinement heuristic
     1476        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
     1477          if (priceRefinement()) continue;
    12091478        }
    12101479
     
    12781547               std::numeric_limits<LargeCost>::max();
    12791548            for (int a = _first_out[n]; a != last_out; ++a) {
    1280               rc = _cost[a] + pi_n - _pi[_target[a]];
    1281               if (_res_cap[a] > 0 && rc < min_red_cost) {
    1282                 min_red_cost = rc;
     1549              if (_res_cap[a] > 0) {
     1550                rc = _cost[a] + pi_n - _pi[_target[a]];
     1551                if (rc < min_red_cost) {
     1552                  min_red_cost = rc;
     1553                }
    12831554              }
    12841555            }
     
    12981569
    12991570          // Global update heuristic
    1300           if (relabel_cnt >= next_update_limit) {
     1571          if (relabel_cnt >= next_global_update_limit) {
    13011572            globalUpdate();
    13021573            for (int u = 0; u != _res_node_num; ++u)
    13031574              hyper[u] = false;
    1304             next_update_limit += global_update_freq;
     1575            next_global_update_limit += global_update_skip;
    13051576          }
    13061577        }
Note: See TracChangeset for help on using the changeset viewer.