COIN-OR::LEMON - Graph Library

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/cost_scaling.h

    r1041 r1271  
    33 * This file is a part of LEMON, a generic C++ optimization library.
    44 *
    5  * Copyright (C) 2003-2010
     5 * Copyright (C) 2003-2013
    66 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    77 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     
    9292  /// \ref CostScaling implements a cost scaling algorithm that performs
    9393  /// push/augment and relabel operations for finding a \ref min_cost_flow
    94   /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
    95   /// \ref goldberg97efficient, \ref bunnagel98efficient.
     94  /// "minimum cost flow" \cite amo93networkflows,
     95  /// \cite goldberg90approximation,
     96  /// \cite goldberg97efficient, \cite bunnagel98efficient.
    9697  /// It is a highly efficient primal-dual solution method, which
    9798  /// can be viewed as the generalization of the \ref Preflow
    9899  /// "preflow push-relabel" algorithm for the maximum flow problem.
     100  /// It is a polynomial algorithm, its running time complexity is
     101  /// \f$O(n^2m\log(nK))\f$, where <i>K</i> denotes the maximum arc cost.
     102  ///
     103  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
     104  /// implementations available in LEMON for solving this problem.
     105  /// (For more information, see \ref min_cost_flow_algs "the module page".)
    99106  ///
    100107  /// Most of the parameters of the problem (except for the digraph)
     
    114121  /// consider to use the named template parameters instead.
    115122  ///
    116   /// \warning Both number types must be signed and all input data must
     123  /// \warning Both \c V and \c C must be signed number types.
     124  /// \warning All input data (capacities, supply values, and costs) must
    117125  /// be integer.
    118   /// \warning This algorithm does not support negative costs for such
    119   /// arcs that have infinite upper bound.
     126  /// \warning This algorithm does not support negative costs for
     127  /// arcs having infinite upper bound.
    120128  ///
    121129  /// \note %CostScaling provides three different internal methods,
     
    146154    typedef typename TR::LargeCost LargeCost;
    147155
    148     /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
     156    /// \brief The \ref lemon::CostScalingDefaultTraits "traits class"
     157    /// of the algorithm
    149158    typedef TR Traits;
    150159
     
    179188    /// relabel operation.
    180189    /// By default, the so called \ref PARTIAL_AUGMENT
    181     /// "Partial Augment-Relabel" method is used, which proved to be
     190    /// "Partial Augment-Relabel" method is used, which turned out to be
    182191    /// the most efficient and the most robust on various test inputs.
    183192    /// However, the other methods can be selected using the \ref run()
     
    206215    typedef std::vector<LargeCost> LargeCostVector;
    207216    typedef std::vector<char> BoolVector;
    208     // Note: vector<char> is used instead of vector<bool> for efficiency reasons
     217    // Note: vector<char> is used instead of vector<bool>
     218    // for efficiency reasons
    209219
    210220  private:
     
    234244    };
    235245
    236     typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
    237246    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
    238247
     
    285294    int _max_rank;
    286295
    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 
    295296  public:
    296297
     
    339340    CostScaling(const GR& graph) :
    340341      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
    341       _cost_map(_cost_vec), _pi_map(_pi),
    342342      INF(std::numeric_limits<Value>::has_infinity ?
    343343          std::numeric_limits<Value>::infinity() :
     
    448448    ///
    449449    /// 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
     450    /// with a map in which \c k is assigned to \c s, \c -k is
    451451    /// assigned to \c t and all other nodes have zero supply value.
    452452    ///
     
    494494    /// \param method The internal method that will be used in the
    495495    /// algorithm. For more information, see \ref Method.
    496     /// \param factor The cost scaling factor. It must be larger than one.
     496    /// \param factor The cost scaling factor. It must be at least two.
    497497    ///
    498498    /// \return \c INFEASIBLE if no feasible flow exists,
     
    508508    /// \see ProblemType, Method
    509509    /// \see resetParams(), reset()
    510     ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
     510    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 16) {
     511      LEMON_ASSERT(factor >= 2, "The scaling factor must be at least 2");
    511512      _alpha = factor;
    512513      ProblemType pt = init();
     
    572573    }
    573574
    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.
     575    /// \brief Reset the internal data structures and all the parameters
     576    /// that have been given before.
     577    ///
     578    /// This function resets the internal data structures and all the
     579    /// paramaters that have been given before using functions \ref lowerMap(),
     580    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     581    ///
     582    /// It is useful for multiple \ref run() calls. By default, all the given
     583    /// parameters are kept for the next \ref run() call, unless
     584    /// \ref resetParams() or \ref reset() is used.
     585    /// If the underlying digraph was also modified after the construction
     586    /// of the class or the last \ref reset() call, then the \ref reset()
     587    /// function must be used, otherwise \ref resetParams() is sufficient.
     588    ///
     589    /// See \ref resetParams() for examples.
     590    ///
    585591    /// \return <tt>(*this)</tt>
     592    ///
     593    /// \see resetParams(), run()
    586594    CostScaling& reset() {
    587595      // Resize vectors
     
    608616      _excess.resize(_res_node_num);
    609617      _next_out.resize(_res_node_num);
    610 
    611       _arc_vec.reserve(_res_arc_num);
    612       _cost_vec.reserve(_res_arc_num);
    613618
    614619      // Copy the graph
     
    668673    ///
    669674    /// This function returns the total cost of the found flow.
    670     /// Its complexity is O(e).
     675    /// Its complexity is O(m).
    671676    ///
    672677    /// \note The return type of the function can be specified as a
     
    706711    }
    707712
    708     /// \brief Return the flow map (the primal solution).
     713    /// \brief Copy the flow values (the primal solution) into the
     714    /// given map.
    709715    ///
    710716    /// This function copies the flow value on each arc into the given
     
    730736    }
    731737
    732     /// \brief Return the potential map (the dual solution).
     738    /// \brief Copy the potential values (the dual solution) into the
     739    /// given map.
    733740    ///
    734741    /// This function copies the potential (dual value) of each node
     
    759766      }
    760767      if (_sum_supply > 0) return INFEASIBLE;
     768
     769      // Check lower and upper bounds
     770      LEMON_DEBUG(checkBoundMaps(),
     771          "Upper bounds must be greater or equal to the lower bounds");
    761772
    762773
     
    887898      }
    888899
    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 
    897900      // Initialize data structures for buckets
    898901      _max_rank = _alpha * _res_node_num;
     
    902905      _rank.resize(_res_node_num + 1);
    903906
    904       // Execute the algorithm
     907      return OPTIMAL;
     908    }
     909
     910    // Check if the upper bound is greater or equal to the lower bound
     911    // on each arc.
     912    bool checkBoundMaps() {
     913      for (int j = 0; j != _res_arc_num; ++j) {
     914        if (_upper[j] < _lower[j]) return false;
     915      }
     916      return true;
     917    }
     918
     919    // Execute the algorithm and transform the results
     920    void start(Method method) {
     921      const int MAX_PARTIAL_PATH_LENGTH = 4;
     922
    905923      switch (method) {
    906924        case PUSH:
     
    911929          break;
    912930        case PARTIAL_AUGMENT:
    913           startAugment(MAX_PATH_LENGTH);
     931          startAugment(MAX_PARTIAL_PATH_LENGTH);
    914932          break;
    915933      }
    916934
    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();
     935      // Compute node potentials (dual solution)
     936      for (int i = 0; i != _res_node_num; ++i) {
     937        _pi[i] = static_cast<Cost>(_pi[i] / (_res_node_num * _alpha));
     938      }
     939      bool optimal = true;
     940      for (int i = 0; optimal && i != _res_node_num; ++i) {
     941        LargeCost pi_i = _pi[i];
     942        int last_out = _first_out[i+1];
     943        for (int j = _first_out[i]; j != last_out; ++j) {
     944          if (_res_cap[j] > 0 && _scost[j] + pi_i - _pi[_target[j]] < 0) {
     945            optimal = false;
     946            break;
     947          }
     948        }
     949      }
     950
     951      if (!optimal) {
     952        // Compute node potentials for the original costs with BellmanFord
     953        // (if it is necessary)
     954        typedef std::pair<int, int> IntPair;
     955        StaticDigraph sgr;
     956        std::vector<IntPair> arc_vec;
     957        std::vector<LargeCost> cost_vec;
     958        LargeCostArcMap cost_map(cost_vec);
     959
     960        arc_vec.clear();
     961        cost_vec.clear();
     962        for (int j = 0; j != _res_arc_num; ++j) {
     963          if (_res_cap[j] > 0) {
     964            int u = _source[j], v = _target[j];
     965            arc_vec.push_back(IntPair(u, v));
     966            cost_vec.push_back(_scost[j] + _pi[u] - _pi[v]);
     967          }
     968        }
     969        sgr.build(_res_node_num, arc_vec.begin(), arc_vec.end());
     970
     971        typename BellmanFord<StaticDigraph, LargeCostArcMap>::Create
     972          bf(sgr, cost_map);
     973        bf.init(0);
     974        bf.start();
     975
     976        for (int i = 0; i != _res_node_num; ++i) {
     977          _pi[i] += bf.dist(sgr.node(i));
     978        }
     979      }
     980
     981      // Shift potentials to meet the requirements of the GEQ type
     982      // optimality conditions
     983      LargeCost max_pot = _pi[_root];
     984      for (int i = 0; i != _res_node_num; ++i) {
     985        if (_pi[i] > max_pot) max_pot = _pi[i];
     986      }
     987      if (max_pot != 0) {
     988        for (int i = 0; i != _res_node_num; ++i) {
     989          _pi[i] -= max_pot;
     990        }
     991      }
    933992
    934993      // Handle non-zero lower bounds
     
    9481007        LargeCost pi_u = _pi[u];
    9491008        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;
     1009          Value delta = _res_cap[a];
     1010          if (delta > 0) {
     1011            int v = _target[a];
     1012            if (_cost[a] + pi_u - _pi[v] < 0) {
     1013              _excess[u] -= delta;
     1014              _excess[v] += delta;
     1015              _res_cap[a] = 0;
     1016              _res_cap[_reverse[a]] += delta;
     1017            }
    9571018          }
    9581019        }
     
    9701031    }
    9711032
    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;
     1033    // Price (potential) refinement heuristic
     1034    bool priceRefinement() {
     1035
     1036      // Stack for stroing the topological order
     1037      IntVector stack(_res_node_num);
     1038      int stack_top;
     1039
     1040      // Perform phases
     1041      while (topologicalSort(stack, stack_top)) {
     1042
     1043        // Compute node ranks in the acyclic admissible network and
     1044        // store the nodes in buckets
     1045        for (int i = 0; i != _res_node_num; ++i) {
     1046          _rank[i] = 0;
     1047        }
     1048        const int bucket_end = _root + 1;
     1049        for (int r = 0; r != _max_rank; ++r) {
     1050          _buckets[r] = bucket_end;
     1051        }
     1052        int top_rank = 0;
     1053        for ( ; stack_top >= 0; --stack_top) {
     1054          int u = stack[stack_top], v;
     1055          int rank_u = _rank[u];
     1056
     1057          LargeCost rc, pi_u = _pi[u];
     1058          int last_out = _first_out[u+1];
     1059          for (int a = _first_out[u]; a != last_out; ++a) {
     1060            if (_res_cap[a] > 0) {
     1061              v = _target[a];
     1062              rc = _cost[a] + pi_u - _pi[v];
     1063              if (rc < 0) {
     1064                LargeCost nrc = static_cast<LargeCost>((-rc - 0.5) / _epsilon);
     1065                if (nrc < LargeCost(_max_rank)) {
     1066                  int new_rank_v = rank_u + static_cast<int>(nrc);
     1067                  if (new_rank_v > _rank[v]) {
     1068                    _rank[v] = new_rank_v;
     1069                  }
     1070                }
     1071              }
     1072            }
     1073          }
     1074
     1075          if (rank_u > 0) {
     1076            top_rank = std::max(top_rank, rank_u);
     1077            int bfirst = _buckets[rank_u];
     1078            _bucket_next[u] = bfirst;
     1079            _bucket_prev[bfirst] = u;
     1080            _buckets[rank_u] = u;
     1081          }
     1082        }
     1083
     1084        // Check if the current flow is epsilon-optimal
     1085        if (top_rank == 0) {
     1086          return true;
     1087        }
     1088
     1089        // Process buckets in top-down order
     1090        for (int rank = top_rank; rank > 0; --rank) {
     1091          while (_buckets[rank] != bucket_end) {
     1092            // Remove the first node from the current bucket
     1093            int u = _buckets[rank];
     1094            _buckets[rank] = _bucket_next[u];
     1095
     1096            // Search the outgoing arcs of u
     1097            LargeCost rc, pi_u = _pi[u];
     1098            int last_out = _first_out[u+1];
     1099            int v, old_rank_v, new_rank_v;
     1100            for (int a = _first_out[u]; a != last_out; ++a) {
     1101              if (_res_cap[a] > 0) {
     1102                v = _target[a];
     1103                old_rank_v = _rank[v];
     1104
     1105                if (old_rank_v < rank) {
     1106
     1107                  // Compute the new rank of node v
     1108                  rc = _cost[a] + pi_u - _pi[v];
     1109                  if (rc < 0) {
     1110                    new_rank_v = rank;
     1111                  } else {
     1112                    LargeCost nrc = rc / _epsilon;
     1113                    new_rank_v = 0;
     1114                    if (nrc < LargeCost(_max_rank)) {
     1115                      new_rank_v = rank - 1 - static_cast<int>(nrc);
     1116                    }
     1117                  }
     1118
     1119                  // Change the rank of node v
     1120                  if (new_rank_v > old_rank_v) {
     1121                    _rank[v] = new_rank_v;
     1122
     1123                    // Remove v from its old bucket
     1124                    if (old_rank_v > 0) {
     1125                      if (_buckets[old_rank_v] == v) {
     1126                        _buckets[old_rank_v] = _bucket_next[v];
     1127                      } else {
     1128                        int pv = _bucket_prev[v], nv = _bucket_next[v];
     1129                        _bucket_next[pv] = nv;
     1130                        _bucket_prev[nv] = pv;
     1131                      }
     1132                    }
     1133
     1134                    // Insert v into its new bucket
     1135                    int nv = _buckets[new_rank_v];
     1136                    _bucket_next[v] = nv;
     1137                    _bucket_prev[nv] = v;
     1138                    _buckets[new_rank_v] = v;
     1139                  }
     1140                }
     1141              }
     1142            }
     1143
     1144            // Refine potential of node u
     1145            _pi[u] -= rank * _epsilon;
     1146          }
     1147        }
     1148
     1149      }
     1150
     1151      return false;
     1152    }
     1153
     1154    // Find and cancel cycles in the admissible network and
     1155    // determine topological order using DFS
     1156    bool topologicalSort(IntVector &stack, int &stack_top) {
     1157      const int MAX_CYCLE_CANCEL = 1;
     1158
     1159      BoolVector reached(_res_node_num, false);
     1160      BoolVector processed(_res_node_num, false);
     1161      IntVector pred(_res_node_num);
     1162      for (int i = 0; i != _res_node_num; ++i) {
     1163        _next_out[i] = _first_out[i];
     1164      }
     1165      stack_top = -1;
     1166
     1167      int cycle_cnt = 0;
     1168      for (int start = 0; start != _res_node_num; ++start) {
     1169        if (reached[start]) continue;
     1170
     1171        // Start DFS search from this start node
     1172        pred[start] = -1;
     1173        int tip = start, v;
     1174        while (true) {
     1175          // Check the outgoing arcs of the current tip node
     1176          reached[tip] = true;
     1177          LargeCost pi_tip = _pi[tip];
     1178          int a, last_out = _first_out[tip+1];
     1179          for (a = _next_out[tip]; a != last_out; ++a) {
     1180            if (_res_cap[a] > 0) {
     1181              v = _target[a];
     1182              if (_cost[a] + pi_tip - _pi[v] < 0) {
     1183                if (!reached[v]) {
     1184                  // A new node is reached
     1185                  reached[v] = true;
     1186                  pred[v] = tip;
     1187                  _next_out[tip] = a;
     1188                  tip = v;
     1189                  a = _next_out[tip];
     1190                  last_out = _first_out[tip+1];
     1191                  break;
     1192                }
     1193                else if (!processed[v]) {
     1194                  // A cycle is found
     1195                  ++cycle_cnt;
     1196                  _next_out[tip] = a;
     1197
     1198                  // Find the minimum residual capacity along the cycle
     1199                  Value d, delta = _res_cap[a];
     1200                  int u, delta_node = tip;
     1201                  for (u = tip; u != v; ) {
     1202                    u = pred[u];
     1203                    d = _res_cap[_next_out[u]];
     1204                    if (d <= delta) {
     1205                      delta = d;
     1206                      delta_node = u;
     1207                    }
     1208                  }
     1209
     1210                  // Augment along the cycle
     1211                  _res_cap[a] -= delta;
     1212                  _res_cap[_reverse[a]] += delta;
     1213                  for (u = tip; u != v; ) {
     1214                    u = pred[u];
     1215                    int ca = _next_out[u];
     1216                    _res_cap[ca] -= delta;
     1217                    _res_cap[_reverse[ca]] += delta;
     1218                  }
     1219
     1220                  // Check the maximum number of cycle canceling
     1221                  if (cycle_cnt >= MAX_CYCLE_CANCEL) {
     1222                    return false;
     1223                  }
     1224
     1225                  // Roll back search to delta_node
     1226                  if (delta_node != tip) {
     1227                    for (u = tip; u != delta_node; u = pred[u]) {
     1228                      reached[u] = false;
     1229                    }
     1230                    tip = delta_node;
     1231                    a = _next_out[tip] + 1;
     1232                    last_out = _first_out[tip+1];
     1233                    break;
     1234                  }
     1235                }
     1236              }
     1237            }
     1238          }
     1239
     1240          // Step back to the previous node
     1241          if (a == last_out) {
     1242            processed[tip] = true;
     1243            stack[++stack_top] = tip;
     1244            tip = pred[tip];
     1245            if (tip < 0) {
     1246              // Finish DFS from the current start node
     1247              break;
     1248            }
     1249            ++_next_out[tip];
     1250          }
     1251        }
     1252
     1253      }
     1254
     1255      return (cycle_cnt == 0);
    9961256    }
    9971257
    9981258    // Global potential update heuristic
    9991259    void globalUpdate() {
    1000       int bucket_end = _root + 1;
     1260      const int bucket_end = _root + 1;
    10011261
    10021262      // Initialize buckets
     
    10051265      }
    10061266      Value total_excess = 0;
     1267      int b0 = bucket_end;
    10071268      for (int i = 0; i != _res_node_num; ++i) {
    10081269        if (_excess[i] < 0) {
    10091270          _rank[i] = 0;
    1010           _bucket_next[i] = _buckets[0];
    1011           _bucket_prev[_buckets[0]] = i;
    1012           _buckets[0] = i;
     1271          _bucket_next[i] = b0;
     1272          _bucket_prev[b0] = i;
     1273          b0 = i;
    10131274        } else {
    10141275          total_excess += _excess[i];
     
    10171278      }
    10181279      if (total_excess == 0) return;
     1280      _buckets[0] = b0;
    10191281
    10201282      // Search the buckets
     
    10261288          _buckets[r] = _bucket_next[u];
    10271289
    1028           // Search the incomming arcs of u
     1290          // Search the incoming arcs of u
    10291291          LargeCost pi_u = _pi[u];
    10301292          int last_out = _first_out[u+1];
     
    10381300                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
    10391301                int new_rank_v = old_rank_v;
    1040                 if (nrc < LargeCost(_max_rank))
    1041                   new_rank_v = r + 1 + int(nrc);
     1302                if (nrc < LargeCost(_max_rank)) {
     1303                  new_rank_v = r + 1 + static_cast<int>(nrc);
     1304                }
    10421305
    10431306                // Change the rank of v
     
    10511314                      _buckets[old_rank_v] = _bucket_next[v];
    10521315                    } else {
    1053                       _bucket_next[_bucket_prev[v]] = _bucket_next[v];
    1054                       _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
     1316                      int pv = _bucket_prev[v], nv = _bucket_next[v];
     1317                      _bucket_next[pv] = nv;
     1318                      _bucket_prev[nv] = pv;
    10551319                    }
    10561320                  }
    10571321
    1058                   // Insert v to its new bucket
    1059                   _bucket_next[v] = _buckets[new_rank_v];
    1060                   _bucket_prev[_buckets[new_rank_v]] = v;
     1322                  // Insert v into its new bucket
     1323                  int nv = _buckets[new_rank_v];
     1324                  _bucket_next[v] = nv;
     1325                  _bucket_prev[nv] = v;
    10611326                  _buckets[new_rank_v] = v;
    10621327                }
     
    10871352    void startAugment(int max_length) {
    10881353      // 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 *
     1354      const int PRICE_REFINEMENT_LIMIT = 2;
     1355      const double GLOBAL_UPDATE_FACTOR = 1.0;
     1356      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
    10931357        (_res_node_num + _sup_node_num * _sup_node_num));
    1094       int next_update_limit = global_update_freq;
    1095 
     1358      int next_global_update_limit = global_update_skip;
     1359
     1360      // Perform cost scaling phases
     1361      IntVector path;
     1362      BoolVector path_arc(_res_arc_num, false);
    10961363      int relabel_cnt = 0;
    1097 
    1098       // Perform cost scaling phases
    1099       std::vector<int> path;
     1364      int eps_phase_cnt = 0;
    11001365      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    11011366                                        1 : _epsilon / _alpha )
    11021367      {
    1103         // Early termination heuristic
    1104         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
    1105           if (earlyTermination()) break;
     1368        ++eps_phase_cnt;
     1369
     1370        // Price refinement heuristic
     1371        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
     1372          if (priceRefinement()) continue;
    11061373        }
    11071374
     
    11201387
    11211388          // Find an augmenting path from the start node
    1122           path.clear();
    11231389          int tip = start;
    1124           while (_excess[tip] >= 0 && int(path.size()) < max_length) {
     1390          while (int(path.size()) < max_length && _excess[tip] >= 0) {
    11251391            int u;
    1126             LargeCost min_red_cost, rc, pi_tip = _pi[tip];
     1392            LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
     1393            LargeCost pi_tip = _pi[tip];
    11271394            int last_out = _first_out[tip+1];
    11281395            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;
     1396              if (_res_cap[a] > 0) {
     1397                u = _target[a];
     1398                rc = _cost[a] + pi_tip - _pi[u];
     1399                if (rc < 0) {
     1400                  path.push_back(a);
     1401                  _next_out[tip] = a;
     1402                  if (path_arc[a]) {
     1403                    goto augment;   // a cycle is found, stop path search
     1404                  }
     1405                  tip = u;
     1406                  path_arc[a] = true;
     1407                  goto next_step;
     1408                }
     1409                else if (rc < min_red_cost) {
     1410                  min_red_cost = rc;
     1411                }
    11351412              }
    11361413            }
    11371414
    11381415            // Relabel tip node
    1139             min_red_cost = std::numeric_limits<LargeCost>::max();
    11401416            if (tip != start) {
    11411417              int ra = _reverse[path.back()];
    1142               min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
     1418              min_red_cost =
     1419                std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
    11431420            }
     1421            last_out = _next_out[tip];
    11441422            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;
     1423              if (_res_cap[a] > 0) {
     1424                rc = _cost[a] + pi_tip - _pi[_target[a]];
     1425                if (rc < min_red_cost) {
     1426                  min_red_cost = rc;
     1427                }
    11481428              }
    11491429            }
     
    11541434            // Step back
    11551435            if (tip != start) {
    1156               tip = _source[path.back()];
     1436              int pa = path.back();
     1437              path_arc[pa] = false;
     1438              tip = _source[pa];
    11571439              path.pop_back();
    11581440            }
     
    11621444
    11631445          // Augment along the found path (as much flow as possible)
     1446        augment:
    11641447          Value delta;
    11651448          int pa, u, v = start;
     
    11681451            u = v;
    11691452            v = _target[pa];
     1453            path_arc[pa] = false;
    11701454            delta = std::min(_res_cap[pa], _excess[u]);
    11711455            _res_cap[pa] -= delta;
     
    11731457            _excess[u] -= delta;
    11741458            _excess[v] += delta;
    1175             if (_excess[v] > 0 && _excess[v] <= delta)
     1459            if (_excess[v] > 0 && _excess[v] <= delta) {
    11761460              _active_nodes.push_back(v);
    1177           }
     1461            }
     1462          }
     1463          path.clear();
    11781464
    11791465          // Global update heuristic
    1180           if (relabel_cnt >= next_update_limit) {
     1466          if (relabel_cnt >= next_global_update_limit) {
    11811467            globalUpdate();
    1182             next_update_limit += global_update_freq;
    1183           }
    1184         }
    1185       }
     1468            next_global_update_limit += global_update_skip;
     1469          }
     1470        }
     1471
     1472      }
     1473
    11861474    }
    11871475
     
    11891477    void startPush() {
    11901478      // Paramters for heuristics
    1191       const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1479      const int PRICE_REFINEMENT_LIMIT = 2;
    11921480      const double GLOBAL_UPDATE_FACTOR = 2.0;
    11931481
    1194       const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1482      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
    11951483        (_res_node_num + _sup_node_num * _sup_node_num));
    1196       int next_update_limit = global_update_freq;
    1197 
    1198       int relabel_cnt = 0;
     1484      int next_global_update_limit = global_update_skip;
    11991485
    12001486      // Perform cost scaling phases
    12011487      BoolVector hyper(_res_node_num, false);
    12021488      LargeCostVector hyper_cost(_res_node_num);
     1489      int relabel_cnt = 0;
     1490      int eps_phase_cnt = 0;
    12031491      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    12041492                                        1 : _epsilon / _alpha )
    12051493      {
    1206         // Early termination heuristic
    1207         if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
    1208           if (earlyTermination()) break;
     1494        ++eps_phase_cnt;
     1495
     1496        // Price refinement heuristic
     1497        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
     1498          if (priceRefinement()) continue;
    12091499        }
    12101500
     
    12781568               std::numeric_limits<LargeCost>::max();
    12791569            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;
     1570              if (_res_cap[a] > 0) {
     1571                rc = _cost[a] + pi_n - _pi[_target[a]];
     1572                if (rc < min_red_cost) {
     1573                  min_red_cost = rc;
     1574                }
    12831575              }
    12841576            }
     
    12981590
    12991591          // Global update heuristic
    1300           if (relabel_cnt >= next_update_limit) {
     1592          if (relabel_cnt >= next_global_update_limit) {
    13011593            globalUpdate();
    13021594            for (int u = 0; u != _res_node_num; ++u)
    13031595              hyper[u] = false;
    1304             next_update_limit += global_update_freq;
     1596            next_global_update_limit += global_update_skip;
    13051597          }
    13061598        }
Note: See TracChangeset for help on using the changeset viewer.