COIN-OR::LEMON - Graph Library

Ignore:
Files:
7 added
19 edited

Legend:

Unmodified
Added
Removed
  • CMakeLists.txt

    r1107 r1111  
    115115SET(LEMON_HAVE_LONG_LONG ${HAVE_LONG_LONG})
    116116
    117 INCLUDE(FindPythonInterp)
    118 
    119117ENABLE_TESTING()
    120118
     
    127125ADD_SUBDIRECTORY(lemon)
    128126IF(${CMAKE_SOURCE_DIR} STREQUAL ${PROJECT_SOURCE_DIR})
     127  ADD_SUBDIRECTORY(contrib)
    129128  ADD_SUBDIRECTORY(demo)
    130129  ADD_SUBDIRECTORY(tools)
  • doc/Doxyfile.in

    r1107 r1111  
    9696                         "@abs_top_srcdir@/lemon/concepts" \
    9797                         "@abs_top_srcdir@/demo" \
     98                         "@abs_top_srcdir@/contrib" \
    9899                         "@abs_top_srcdir@/tools" \
    99100                         "@abs_top_srcdir@/test/test_tools.h" \
  • doc/coding_style.dox

    r463 r1023  
    9999\subsection pri-loc-var Private member variables
    100100
    101 Private member variables should start with underscore
     101Private member variables should start with underscore.
    102102
    103103\code
    104 _start_with_underscores
     104_start_with_underscore
    105105\endcode
    106106
  • doc/dirs.dox

    r463 r1031  
    3232documentation.
    3333*/
     34
     35/**
     36\dir contrib
     37\brief Directory for user contributed source codes.
     38
     39You can place your own C++ code using LEMON into this directory, which
     40will compile to an executable along with LEMON when you build the
     41library. This is probably the easiest way of compiling short to medium
     42codes, for this does require neither a LEMON installed system-wide nor
     43adding several paths to the compiler.
     44
     45Please have a look at <tt>contrib/CMakeLists.txt</tt> for
     46instruction on how to add your own files into the build process.  */
    3447
    3548/**
  • doc/groups.dox

    r959 r1023  
    407407   strongly polynomial \ref klein67primal, \ref goldberg89cyclecanceling.
    408408
    409 In general NetworkSimplex is the most efficient implementation,
    410 but in special cases other algorithms could be faster.
     409In general, \ref NetworkSimplex and \ref CostScaling are the most efficient
     410implementations, but the other two algorithms could be faster in special cases.
    411411For example, if the total supply and/or capacities are rather small,
    412 CapacityScaling is usually the fastest algorithm (without effective scaling).
     412\ref CapacityScaling is usually the fastest algorithm (without effective scaling).
    413413*/
    414414
     
    472472  \ref dasdan98minmeancycle.
    473473
    474 In practice, the \ref HowardMmc "Howard" algorithm proved to be by far the
     474In practice, the \ref HowardMmc "Howard" algorithm turned out to be by far the
    475475most efficient one, though the best known theoretical bound on its running
    476476time is exponential.
     
    540540
    541541/**
    542 @defgroup planar Planarity Embedding and Drawing
     542@defgroup planar Planar Embedding and Drawing
    543543@ingroup algs
    544544\brief Algorithms for planarity checking, embedding and drawing
     
    552552
    553553/**
    554 @defgroup approx Approximation Algorithms
     554@defgroup approx_algs Approximation Algorithms
    555555@ingroup algs
    556556\brief Approximation algorithms.
     
    558558This group contains the approximation and heuristic algorithms
    559559implemented in LEMON.
     560
     561<b>Maximum Clique Problem</b>
     562  - \ref GrossoLocatelliPullanMc An efficient heuristic algorithm of
     563    Grosso, Locatelli, and Pullan.
    560564*/
    561565
  • doc/references.bib

    r802 r999  
    298298  address =      {Dublin, Ireland},
    299299  year =         1991,
    300   month =        sep,
    301 }
     300  month =        sep
     301}
     302
     303%%%%% Other algorithms %%%%%
     304
     305@article{grosso08maxclique,
     306  author =       {Andrea Grosso and Marco Locatelli and Wayne Pullan},
     307  title =        {Simple ingredients leading to very efficient
     308                  heuristics for the maximum clique problem},
     309  journal =      {Journal of Heuristics},
     310  year =         2008,
     311  volume =       14,
     312  number =       6,
     313  pages =        {587--612}
     314}
  • lemon/Makefile.am

    r1107 r1111  
    9292        lemon/graph_to_eps.h \
    9393        lemon/grid_graph.h \
     94        lemon/grosso_locatelli_pullan_mc.h \
    9495        lemon/hartmann_orlin_mmc.h \
    9596        lemon/howard_mmc.h \
     
    108109        lemon/math.h \
    109110        lemon/min_cost_arborescence.h \
     111        lemon/max_cardinality_search.h \
     112        lemon/nagamochi_ibaraki.h \
    110113        lemon/nauty_reader.h \
    111114        lemon/network_simplex.h \
  • lemon/capacity_scaling.h

    r956 r1026  
    8787  /// consider to use the named template parameters instead.
    8888  ///
    89   /// \warning Both number types must be signed and all input data must
     89  /// \warning Both \c V and \c C must be signed number types.
     90  /// \warning All input data (capacities, supply values, and costs) must
    9091  /// be integer.
    91   /// \warning This algorithm does not support negative costs for such
    92   /// arcs that have infinite upper bound.
     92  /// \warning This algorithm does not support negative costs for
     93  /// arcs having infinite upper bound.
    9394#ifdef DOXYGEN
    9495  template <typename GR, typename V, typename C, typename TR>
     
    423424    ///
    424425    /// Using this function has the same effect as using \ref supplyMap()
    425     /// with such a map in which \c k is assigned to \c s, \c -k is
     426    /// with a map in which \c k is assigned to \c s, \c -k is
    426427    /// assigned to \c t and all other nodes have zero supply value.
    427428    ///
  • lemon/core.h

    r1107 r1111  
    447447
    448448  }
     449
     450  /// \brief Check whether a graph is undirected.
     451  ///
     452  /// This function returns \c true if the given graph is undirected.
     453#ifdef DOXYGEN
     454  template <typename GR>
     455  bool undirected(const GR& g) { return false; }
     456#else
     457  template <typename GR>
     458  typename enable_if<UndirectedTagIndicator<GR>, bool>::type
     459  undirected(const GR&) {
     460    return true;
     461  }
     462  template <typename GR>
     463  typename disable_if<UndirectedTagIndicator<GR>, bool>::type
     464  undirected(const GR&) {
     465    return false;
     466  }
     467#endif
    449468
    450469  /// \brief Class to copy a digraph.
  • 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        }
  • lemon/cycle_canceling.h

    r956 r1026  
    6666  /// algorithm. By default, it is the same as \c V.
    6767  ///
    68   /// \warning Both number types must be signed and all input data must
     68  /// \warning Both \c V and \c C must be signed number types.
     69  /// \warning All input data (capacities, supply values, and costs) must
    6970  /// be integer.
    70   /// \warning This algorithm does not support negative costs for such
    71   /// arcs that have infinite upper bound.
     71  /// \warning This algorithm does not support negative costs for
     72  /// arcs having infinite upper bound.
    7273  ///
    7374  /// \note For more information about the three available methods,
     
    117118    /// \ref CycleCanceling provides three different cycle-canceling
    118119    /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
    119     /// is used, which proved to be the most efficient and the most robust
    120     /// on various test inputs.
     120    /// is used, which is by far the most efficient and the most robust.
    121121    /// However, the other methods can be selected using the \ref run()
    122122    /// function with the proper parameter.
     
    350350    ///
    351351    /// Using this function has the same effect as using \ref supplyMap()
    352     /// with such a map in which \c k is assigned to \c s, \c -k is
     352    /// with a map in which \c k is assigned to \c s, \c -k is
    353353    /// assigned to \c t and all other nodes have zero supply value.
    354354    ///
  • lemon/euler.h

    r956 r1023  
    3737  ///Euler tour iterator for digraphs.
    3838
    39   /// \ingroup graph_prop
     39  /// \ingroup graph_properties
    4040  ///This iterator provides an Euler tour (Eulerian circuit) of a \e directed
    4141  ///graph (if there exists) and it converts to the \c Arc type of the digraph.
  • lemon/hao_orlin.h

    r956 r1019  
    5454  /// preflow push-relabel algorithm. Our implementation calculates
    5555  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
    56   /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
    57   /// purpose of such algorithm is e.g. testing network reliability.
     56  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. A notable
     57  /// use of this algorithm is testing network reliability.
    5858  ///
    5959  /// For an undirected graph you can run just the first phase of the
     
    913913    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with
    914914    /// \f$ source \in X \f$ and minimal outgoing capacity).
     915    /// It updates the stored cut if (and only if) the newly found one
     916    /// is better.
    915917    ///
    916918    /// \pre \ref init() must be called before using this function.
     
    925927    /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with
    926928    /// \f$ source \notin X \f$ and minimal outgoing capacity).
     929    /// It updates the stored cut if (and only if) the newly found one
     930    /// is better.
    927931    ///
    928932    /// \pre \ref init() must be called before using this function.
     
    934938    /// \brief Run the algorithm.
    935939    ///
    936     /// This function runs the algorithm. It finds nodes \c source and
    937     /// \c target arbitrarily and then calls \ref init(), \ref calculateOut()
     940    /// This function runs the algorithm. It chooses source node,
     941    /// then calls \ref init(), \ref calculateOut()
    938942    /// and \ref calculateIn().
    939943    void run() {
     
    945949    /// \brief Run the algorithm.
    946950    ///
    947     /// This function runs the algorithm. It uses the given \c source node,
    948     /// finds a proper \c target node and then calls the \ref init(),
    949     /// \ref calculateOut() and \ref calculateIn().
     951    /// This function runs the algorithm. It calls \ref init(),
     952    /// \ref calculateOut() and \ref calculateIn() with the given
     953    /// source node.
    950954    void run(const Node& s) {
    951955      init(s);
     
    966970    /// \brief Return the value of the minimum cut.
    967971    ///
    968     /// This function returns the value of the minimum cut.
     972    /// This function returns the value of the best cut found by the
     973    /// previously called \ref run(), \ref calculateOut() or \ref
     974    /// calculateIn().
    969975    ///
    970976    /// \pre \ref run(), \ref calculateOut() or \ref calculateIn()
     
    977983    /// \brief Return a minimum cut.
    978984    ///
    979     /// This function sets \c cutMap to the characteristic vector of a
    980     /// minimum value cut: it will give a non-empty set \f$ X\subsetneq V \f$
    981     /// with minimal outgoing capacity (i.e. \c cutMap will be \c true exactly
     985    /// This function gives the best cut found by the
     986    /// previously called \ref run(), \ref calculateOut() or \ref
     987    /// calculateIn().
     988    ///
     989    /// It sets \c cutMap to the characteristic vector of the found
     990    /// minimum value cut - a non-empty set \f$ X\subsetneq V \f$
     991    /// of minimum outgoing capacity (i.e. \c cutMap will be \c true exactly
    982992    /// for the nodes of \f$ X \f$).
    983993    ///
  • lemon/kruskal.h

    r631 r1025  
    3131///\file
    3232///\brief Kruskal's algorithm to compute a minimum cost spanning tree
    33 ///
    34 ///Kruskal's algorithm to compute a minimum cost spanning tree.
    35 ///
    3633
    3734namespace lemon {
  • lemon/network_simplex.h

    r978 r1026  
    4848  /// flow problem.
    4949  ///
    50   /// In general, %NetworkSimplex is the fastest implementation available
    51   /// in LEMON for this problem.
    52   /// Moreover, it supports both directions of the supply/demand inequality
    53   /// constraints. For more information, see \ref SupplyType.
     50  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
     51  /// implementations available in LEMON for this problem.
     52  /// Furthermore, this class supports both directions of the supply/demand
     53  /// inequality constraints. For more information, see \ref SupplyType.
    5454  ///
    5555  /// Most of the parameters of the problem (except for the digraph)
     
    6464  /// algorithm. By default, it is the same as \c V.
    6565  ///
    66   /// \warning Both number types must be signed and all input data must
     66  /// \warning Both \c V and \c C must be signed number types.
     67  /// \warning All input data (capacities, supply values, and costs) must
    6768  /// be integer.
    6869  ///
     
    126127    /// of the algorithm.
    127128    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
    128     /// proved to be the most efficient and the most robust on various
     129    /// turend out to be the most efficient and the most robust on various
    129130    /// test inputs.
    130131    /// However, another pivot rule can be selected using the \ref run()
     
    167168    typedef std::vector<Value> ValueVector;
    168169    typedef std::vector<Cost> CostVector;
    169     typedef std::vector<char> BoolVector;
    170     // Note: vector<char> is used instead of vector<bool> for efficiency reasons
     170    typedef std::vector<signed char> CharVector;
     171    // Note: vector<signed char> is used instead of vector<ArcState> and
     172    // vector<ArcDirection> for efficiency reasons
    171173
    172174    // State constants for arcs
     
    177179    };
    178180
    179     typedef std::vector<signed char> StateVector;
    180     // Note: vector<signed char> is used instead of vector<ArcState> for
    181     // efficiency reasons
     181    // Direction constants for tree arcs
     182    enum ArcDirection {
     183      DIR_DOWN = -1,
     184      DIR_UP   =  1
     185    };
    182186
    183187  private:
     
    218222    IntVector _succ_num;
    219223    IntVector _last_succ;
     224    CharVector _pred_dir;
     225    CharVector _state;
    220226    IntVector _dirty_revs;
    221     BoolVector _forward;
    222     StateVector _state;
    223227    int _root;
    224228
    225229    // Temporary data used in the current pivot iteration
    226230    int in_arc, join, u_in, v_in, u_out, v_out;
    227     int first, second, right, last;
    228     int stem, par_stem, new_stem;
    229231    Value delta;
    230232
     
    251253      const IntVector  &_target;
    252254      const CostVector &_cost;
    253       const StateVector &_state;
     255      const CharVector &_state;
    254256      const CostVector &_pi;
    255257      int &_in_arc;
     
    303305      const IntVector  &_target;
    304306      const CostVector &_cost;
    305       const StateVector &_state;
     307      const CharVector &_state;
    306308      const CostVector &_pi;
    307309      int &_in_arc;
     
    342344      const IntVector  &_target;
    343345      const CostVector &_cost;
    344       const StateVector &_state;
     346      const CharVector &_state;
    345347      const CostVector &_pi;
    346348      int &_in_arc;
     
    415417      const IntVector  &_target;
    416418      const CostVector &_cost;
    417       const StateVector &_state;
     419      const CharVector &_state;
    418420      const CostVector &_pi;
    419421      int &_in_arc;
     
    518520      const IntVector  &_target;
    519521      const CostVector &_cost;
    520       const StateVector &_state;
     522      const CharVector &_state;
    521523      const CostVector &_pi;
    522524      int &_in_arc;
     
    571573        // Check the current candidate list
    572574        int e;
     575        Cost c;
    573576        for (int i = 0; i != _curr_length; ++i) {
    574577          e = _candidates[i];
    575           _cand_cost[e] = _state[e] *
    576             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    577           if (_cand_cost[e] >= 0) {
     578          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     579          if (c < 0) {
     580            _cand_cost[e] = c;
     581          } else {
    578582            _candidates[i--] = _candidates[--_curr_length];
    579583          }
     
    585589
    586590        for (e = _next_arc; e != _search_arc_num; ++e) {
    587           _cand_cost[e] = _state[e] *
    588             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    589           if (_cand_cost[e] < 0) {
     591          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     592          if (c < 0) {
     593            _cand_cost[e] = c;
    590594            _candidates[_curr_length++] = e;
    591595          }
     
    634638    ///
    635639    /// \param graph The digraph the algorithm runs on.
    636     /// \param arc_mixing Indicate if the arcs have to be stored in a
     640    /// \param arc_mixing Indicate if the arcs will be stored in a
    637641    /// mixed order in the internal data structure.
    638     /// In special cases, it could lead to better overall performance,
    639     /// but it is usually slower. Therefore it is disabled by default.
    640     NetworkSimplex(const GR& graph, bool arc_mixing = false) :
     642    /// In general, it leads to similar performance as using the original
     643    /// arc order, but it makes the algorithm more robust and in special
     644    /// cases, even significantly faster. Therefore, it is enabled by default.
     645    NetworkSimplex(const GR& graph, bool arc_mixing = true) :
    641646      _graph(graph), _node_id(graph), _arc_id(graph),
    642647      _arc_mixing(arc_mixing),
     
    731736    ///
    732737    /// \return <tt>(*this)</tt>
     738    ///
     739    /// \sa supplyType()
    733740    template<typename SupplyMap>
    734741    NetworkSimplex& supplyMap(const SupplyMap& map) {
     
    747754    ///
    748755    /// Using this function has the same effect as using \ref supplyMap()
    749     /// with such a map in which \c k is assigned to \c s, \c -k is
     756    /// with a map in which \c k is assigned to \c s, \c -k is
    750757    /// assigned to \c t and all other nodes have zero supply value.
    751758    ///
     
    914921      _parent.resize(all_node_num);
    915922      _pred.resize(all_node_num);
    916       _forward.resize(all_node_num);
     923      _pred_dir.resize(all_node_num);
    917924      _thread.resize(all_node_num);
    918925      _rev_thread.resize(all_node_num);
     
    928935      if (_arc_mixing) {
    929936        // Store the arcs in a mixed order
    930         int k = std::max(int(std::sqrt(double(_arc_num))), 10);
     937        const int skip = std::max(_arc_num / _node_num, 3);
    931938        int i = 0, j = 0;
    932939        for (ArcIt a(_graph); a != INVALID; ++a) {
     
    934941          _source[i] = _node_id[_graph.source(a)];
    935942          _target[i] = _node_id[_graph.target(a)];
    936           if ((i += k) >= _arc_num) i = ++j;
     943          if ((i += skip) >= _arc_num) i = ++j;
    937944        }
    938945      } else {
     
    11171124          _state[e] = STATE_TREE;
    11181125          if (_supply[u] >= 0) {
    1119             _forward[u] = true;
     1126            _pred_dir[u] = DIR_UP;
    11201127            _pi[u] = 0;
    11211128            _source[e] = u;
     
    11241131            _cost[e] = 0;
    11251132          } else {
    1126             _forward[u] = false;
     1133            _pred_dir[u] = DIR_DOWN;
    11271134            _pi[u] = ART_COST;
    11281135            _source[e] = _root;
     
    11441151          _last_succ[u] = u;
    11451152          if (_supply[u] >= 0) {
    1146             _forward[u] = true;
     1153            _pred_dir[u] = DIR_UP;
    11471154            _pi[u] = 0;
    11481155            _pred[u] = e;
     
    11541161            _state[e] = STATE_TREE;
    11551162          } else {
    1156             _forward[u] = false;
     1163            _pred_dir[u] = DIR_DOWN;
    11571164            _pi[u] = ART_COST;
    11581165            _pred[u] = f;
     
    11851192          _last_succ[u] = u;
    11861193          if (_supply[u] <= 0) {
    1187             _forward[u] = false;
     1194            _pred_dir[u] = DIR_DOWN;
    11881195            _pi[u] = 0;
    11891196            _pred[u] = e;
     
    11951202            _state[e] = STATE_TREE;
    11961203          } else {
    1197             _forward[u] = true;
     1204            _pred_dir[u] = DIR_UP;
    11981205            _pi[u] = -ART_COST;
    11991206            _pred[u] = f;
     
    12381245      // Initialize first and second nodes according to the direction
    12391246      // of the cycle
     1247      int first, second;
    12401248      if (_state[in_arc] == STATE_LOWER) {
    12411249        first  = _source[in_arc];
     
    12471255      delta = _cap[in_arc];
    12481256      int result = 0;
    1249       Value d;
     1257      Value c, d;
    12501258      int e;
    12511259
    1252       // Search the cycle along the path form the first node to the root
     1260      // Search the cycle form the first node to the join node
    12531261      for (int u = first; u != join; u = _parent[u]) {
    12541262        e = _pred[u];
    1255         d = _forward[u] ?
    1256           _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
     1263        d = _flow[e];
     1264        if (_pred_dir[u] == DIR_DOWN) {
     1265          c = _cap[e];
     1266          d = c >= MAX ? INF : c - d;
     1267        }
    12571268        if (d < delta) {
    12581269          delta = d;
     
    12611272        }
    12621273      }
    1263       // Search the cycle along the path form the second node to the root
     1274
     1275      // Search the cycle form the second node to the join node
    12641276      for (int u = second; u != join; u = _parent[u]) {
    12651277        e = _pred[u];
    1266         d = _forward[u] ?
    1267           (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
     1278        d = _flow[e];
     1279        if (_pred_dir[u] == DIR_UP) {
     1280          c = _cap[e];
     1281          d = c >= MAX ? INF : c - d;
     1282        }
    12681283        if (d <= delta) {
    12691284          delta = d;
     
    12901305        _flow[in_arc] += val;
    12911306        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
    1292           _flow[_pred[u]] += _forward[u] ? -val : val;
     1307          _flow[_pred[u]] -= _pred_dir[u] * val;
    12931308        }
    12941309        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
    1295           _flow[_pred[u]] += _forward[u] ? val : -val;
     1310          _flow[_pred[u]] += _pred_dir[u] * val;
    12961311        }
    12971312      }
     
    13081323    // Update the tree structure
    13091324    void updateTreeStructure() {
    1310       int u, w;
    13111325      int old_rev_thread = _rev_thread[u_out];
    13121326      int old_succ_num = _succ_num[u_out];
     
    13141328      v_out = _parent[u_out];
    13151329
    1316       u = _last_succ[u_in];  // the last successor of u_in
    1317       right = _thread[u];    // the node after it
    1318 
    1319       // Handle the case when old_rev_thread equals to v_in
    1320       // (it also means that join and v_out coincide)
    1321       if (old_rev_thread == v_in) {
    1322         last = _thread[_last_succ[u_out]];
     1330      // Check if u_in and u_out coincide
     1331      if (u_in == u_out) {
     1332        // Update _parent, _pred, _pred_dir
     1333        _parent[u_in] = v_in;
     1334        _pred[u_in] = in_arc;
     1335        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
     1336
     1337        // Update _thread and _rev_thread
     1338        if (_thread[v_in] != u_out) {
     1339          int after = _thread[old_last_succ];
     1340          _thread[old_rev_thread] = after;
     1341          _rev_thread[after] = old_rev_thread;
     1342          after = _thread[v_in];
     1343          _thread[v_in] = u_out;
     1344          _rev_thread[u_out] = v_in;
     1345          _thread[old_last_succ] = after;
     1346          _rev_thread[after] = old_last_succ;
     1347        }
    13231348      } else {
    1324         last = _thread[v_in];
    1325       }
    1326 
    1327       // Update _thread and _parent along the stem nodes (i.e. the nodes
    1328       // between u_in and u_out, whose parent have to be changed)
    1329       _thread[v_in] = stem = u_in;
    1330       _dirty_revs.clear();
    1331       _dirty_revs.push_back(v_in);
    1332       par_stem = v_in;
    1333       while (stem != u_out) {
    1334         // Insert the next stem node into the thread list
    1335         new_stem = _parent[stem];
    1336         _thread[u] = new_stem;
    1337         _dirty_revs.push_back(u);
    1338 
    1339         // Remove the subtree of stem from the thread list
    1340         w = _rev_thread[stem];
    1341         _thread[w] = right;
    1342         _rev_thread[right] = w;
    1343 
    1344         // Change the parent node and shift stem nodes
    1345         _parent[stem] = par_stem;
    1346         par_stem = stem;
    1347         stem = new_stem;
    1348 
    1349         // Update u and right
    1350         u = _last_succ[stem] == _last_succ[par_stem] ?
    1351           _rev_thread[par_stem] : _last_succ[stem];
    1352         right = _thread[u];
    1353       }
    1354       _parent[u_out] = par_stem;
    1355       _thread[u] = last;
    1356       _rev_thread[last] = u;
    1357       _last_succ[u_out] = u;
    1358 
    1359       // Remove the subtree of u_out from the thread list except for
    1360       // the case when old_rev_thread equals to v_in
    1361       // (it also means that join and v_out coincide)
    1362       if (old_rev_thread != v_in) {
    1363         _thread[old_rev_thread] = right;
    1364         _rev_thread[right] = old_rev_thread;
    1365       }
    1366 
    1367       // Update _rev_thread using the new _thread values
    1368       for (int i = 0; i != int(_dirty_revs.size()); ++i) {
    1369         u = _dirty_revs[i];
    1370         _rev_thread[_thread[u]] = u;
    1371       }
    1372 
    1373       // Update _pred, _forward, _last_succ and _succ_num for the
    1374       // stem nodes from u_out to u_in
    1375       int tmp_sc = 0, tmp_ls = _last_succ[u_out];
    1376       u = u_out;
    1377       while (u != u_in) {
    1378         w = _parent[u];
    1379         _pred[u] = _pred[w];
    1380         _forward[u] = !_forward[w];
    1381         tmp_sc += _succ_num[u] - _succ_num[w];
    1382         _succ_num[u] = tmp_sc;
    1383         _last_succ[w] = tmp_ls;
    1384         u = w;
    1385       }
    1386       _pred[u_in] = in_arc;
    1387       _forward[u_in] = (u_in == _source[in_arc]);
    1388       _succ_num[u_in] = old_succ_num;
    1389 
    1390       // Set limits for updating _last_succ form v_in and v_out
    1391       // towards the root
    1392       int up_limit_in = -1;
    1393       int up_limit_out = -1;
    1394       if (_last_succ[join] == v_in) {
    1395         up_limit_out = join;
    1396       } else {
    1397         up_limit_in = join;
     1349        // Handle the case when old_rev_thread equals to v_in
     1350        // (it also means that join and v_out coincide)
     1351        int thread_continue = old_rev_thread == v_in ?
     1352          _thread[old_last_succ] : _thread[v_in];
     1353
     1354        // Update _thread and _parent along the stem nodes (i.e. the nodes
     1355        // between u_in and u_out, whose parent have to be changed)
     1356        int stem = u_in;              // the current stem node
     1357        int par_stem = v_in;          // the new parent of stem
     1358        int next_stem;                // the next stem node
     1359        int last = _last_succ[u_in];  // the last successor of stem
     1360        int before, after = _thread[last];
     1361        _thread[v_in] = u_in;
     1362        _dirty_revs.clear();
     1363        _dirty_revs.push_back(v_in);
     1364        while (stem != u_out) {
     1365          // Insert the next stem node into the thread list
     1366          next_stem = _parent[stem];
     1367          _thread[last] = next_stem;
     1368          _dirty_revs.push_back(last);
     1369
     1370          // Remove the subtree of stem from the thread list
     1371          before = _rev_thread[stem];
     1372          _thread[before] = after;
     1373          _rev_thread[after] = before;
     1374
     1375          // Change the parent node and shift stem nodes
     1376          _parent[stem] = par_stem;
     1377          par_stem = stem;
     1378          stem = next_stem;
     1379
     1380          // Update last and after
     1381          last = _last_succ[stem] == _last_succ[par_stem] ?
     1382            _rev_thread[par_stem] : _last_succ[stem];
     1383          after = _thread[last];
     1384        }
     1385        _parent[u_out] = par_stem;
     1386        _thread[last] = thread_continue;
     1387        _rev_thread[thread_continue] = last;
     1388        _last_succ[u_out] = last;
     1389
     1390        // Remove the subtree of u_out from the thread list except for
     1391        // the case when old_rev_thread equals to v_in
     1392        if (old_rev_thread != v_in) {
     1393          _thread[old_rev_thread] = after;
     1394          _rev_thread[after] = old_rev_thread;
     1395        }
     1396
     1397        // Update _rev_thread using the new _thread values
     1398        for (int i = 0; i != int(_dirty_revs.size()); ++i) {
     1399          int u = _dirty_revs[i];
     1400          _rev_thread[_thread[u]] = u;
     1401        }
     1402
     1403        // Update _pred, _pred_dir, _last_succ and _succ_num for the
     1404        // stem nodes from u_out to u_in
     1405        int tmp_sc = 0, tmp_ls = _last_succ[u_out];
     1406        for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
     1407          _pred[u] = _pred[p];
     1408          _pred_dir[u] = -_pred_dir[p];
     1409          tmp_sc += _succ_num[u] - _succ_num[p];
     1410          _succ_num[u] = tmp_sc;
     1411          _last_succ[p] = tmp_ls;
     1412        }
     1413        _pred[u_in] = in_arc;
     1414        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
     1415        _succ_num[u_in] = old_succ_num;
    13981416      }
    13991417
    14001418      // Update _last_succ from v_in towards the root
    1401       for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
    1402            u = _parent[u]) {
    1403         _last_succ[u] = _last_succ[u_out];
    1404       }
     1419      int up_limit_out = _last_succ[join] == v_in ? join : -1;
     1420      int last_succ_out = _last_succ[u_out];
     1421      for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
     1422        _last_succ[u] = last_succ_out;
     1423      }
     1424
    14051425      // Update _last_succ from v_out towards the root
    14061426      if (join != old_rev_thread && v_in != old_rev_thread) {
    1407         for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
     1427        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
    14081428             u = _parent[u]) {
    14091429          _last_succ[u] = old_rev_thread;
    14101430        }
    1411       } else {
    1412         for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
     1431      }
     1432      else if (last_succ_out != old_last_succ) {
     1433        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
    14131434             u = _parent[u]) {
    1414           _last_succ[u] = _last_succ[u_out];
     1435          _last_succ[u] = last_succ_out;
    14151436        }
    14161437      }
    14171438
    14181439      // Update _succ_num from v_in to join
    1419       for (u = v_in; u != join; u = _parent[u]) {
     1440      for (int u = v_in; u != join; u = _parent[u]) {
    14201441        _succ_num[u] += old_succ_num;
    14211442      }
    14221443      // Update _succ_num from v_out to join
    1423       for (u = v_out; u != join; u = _parent[u]) {
     1444      for (int u = v_out; u != join; u = _parent[u]) {
    14241445        _succ_num[u] -= old_succ_num;
    14251446      }
    14261447    }
    14271448
    1428     // Update potentials
     1449    // Update potentials in the subtree that has been moved
    14291450    void updatePotential() {
    1430       Cost sigma = _forward[u_in] ?
    1431         _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
    1432         _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
    1433       // Update potentials in the subtree, which has been moved
     1451      Cost sigma = _pi[v_in] - _pi[u_in] -
     1452                   _pred_dir[u_in] * _cost[in_arc];
    14341453      int end = _thread[_last_succ[u_in]];
    14351454      for (int u = u_in; u != end; u = _thread[u]) {
  • lemon/path.h

    r956 r1024  
    4444  ///
    4545  /// In a sense, the path can be treated as a list of arcs. The
    46   /// lemon path type stores just this list. As a consequence, it
     46  /// LEMON path type stores just this list. As a consequence, it
    4747  /// cannot enumerate the nodes of the path and the source node of
    4848  /// a zero length path is undefined.
     
    136136    void clear() { head.clear(); tail.clear(); }
    137137
    138     /// \brief The nth arc.
     138    /// \brief The n-th arc.
    139139    ///
    140140    /// \pre \c n is in the <tt>[0..length() - 1]</tt> range.
     
    144144    }
    145145
    146     /// \brief Initialize arc iterator to point to the nth arc
     146    /// \brief Initialize arc iterator to point to the n-th arc
    147147    ///
    148148    /// \pre \c n is in the <tt>[0..length() - 1]</tt> range.
     
    232232  ///
    233233  /// In a sense, the path can be treated as a list of arcs. The
    234   /// lemon path type stores just this list. As a consequence it
     234  /// LEMON path type stores just this list. As a consequence it
    235235  /// cannot enumerate the nodes in the path and the zero length paths
    236236  /// cannot store the source.
     
    328328    void clear() { data.clear(); }
    329329
    330     /// \brief The nth arc.
     330    /// \brief The n-th arc.
    331331    ///
    332332    /// \pre \c n is in the <tt>[0..length() - 1]</tt> range.
     
    335335    }
    336336
    337     /// \brief  Initializes arc iterator to point to the nth arc.
     337    /// \brief  Initializes arc iterator to point to the n-th arc.
    338338    ArcIt nthIt(int n) const {
    339339      return ArcIt(*this, n);
     
    396396  ///
    397397  /// In a sense, the path can be treated as a list of arcs. The
    398   /// lemon path type stores just this list. As a consequence it
     398  /// LEMON path type stores just this list. As a consequence it
    399399  /// cannot enumerate the nodes in the path and the zero length paths
    400400  /// cannot store the source.
     
    505505    };
    506506
    507     /// \brief The nth arc.
    508     ///
    509     /// This function looks for the nth arc in O(n) time.
     507    /// \brief The n-th arc.
     508    ///
     509    /// This function looks for the n-th arc in O(n) time.
    510510    /// \pre \c n is in the <tt>[0..length() - 1]</tt> range.
    511511    const Arc& nth(int n) const {
     
    517517    }
    518518
    519     /// \brief Initializes arc iterator to point to the nth arc.
     519    /// \brief Initializes arc iterator to point to the n-th arc.
    520520    ArcIt nthIt(int n) const {
    521521      Node *node = first;
     
    736736  ///
    737737  /// In a sense, the path can be treated as a list of arcs. The
    738   /// lemon path type stores just this list. As a consequence it
     738  /// LEMON path type stores just this list. As a consequence it
    739739  /// cannot enumerate the nodes in the path and the source node of
    740740  /// a zero length path is undefined.
     
    832832    };
    833833
    834     /// \brief The nth arc.
     834    /// \brief The n-th arc.
    835835    ///
    836836    /// \pre \c n is in the <tt>[0..length() - 1]</tt> range.
     
    839839    }
    840840
    841     /// \brief The arc iterator pointing to the nth arc.
     841    /// \brief The arc iterator pointing to the n-th arc.
    842842    ArcIt nthIt(int n) const {
    843843      return ArcIt(*this, n);
     
    10431043  ///
    10441044  /// In a sense, the path can be treated as a list of arcs. The
    1045   /// lemon path type stores only this list. As a consequence, it
     1045  /// LEMON path type stores only this list. As a consequence, it
    10461046  /// cannot enumerate the nodes in the path and the zero length paths
    10471047  /// cannot have a source node.
  • test/CMakeLists.txt

    r1107 r1111  
    3737  maps_test
    3838  matching_test
     39  max_cardinality_search_test
     40  max_clique_test
    3941  min_cost_arborescence_test
    4042  min_cost_flow_test
    4143  min_mean_cycle_test
     44  nagamochi_ibaraki_test
    4245  path_test
    4346  planarity_test
  • test/Makefile.am

    r1107 r1111  
    3535        test/maps_test \
    3636        test/matching_test \
     37        test/max_cardinality_search_test \
     38        test/max_clique_test \
    3739        test/min_cost_arborescence_test \
    3840        test/min_cost_flow_test \
    3941        test/min_mean_cycle_test \
     42        test/nagamochi_ibaraki_test \
    4043        test/path_test \
    4144        test/planarity_test \
     
    7982test_graph_test_SOURCES = test/graph_test.cc
    8083test_graph_utils_test_SOURCES = test/graph_utils_test.cc
     84test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
    8185test_heap_test_SOURCES = test/heap_test.cc
    8286test_kruskal_test_SOURCES = test/kruskal_test.cc
    83 test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
    8487test_lgf_test_SOURCES = test/lgf_test.cc
    8588test_lp_test_SOURCES = test/lp_test.cc
     
    8790test_mip_test_SOURCES = test/mip_test.cc
    8891test_matching_test_SOURCES = test/matching_test.cc
     92test_max_cardinality_search_test_SOURCES = test/max_cardinality_search_test.cc
     93test_max_clique_test_SOURCES = test/max_clique_test.cc
    8994test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
    9095test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
    9196test_min_mean_cycle_test_SOURCES = test/min_mean_cycle_test.cc
     97test_nagamochi_ibaraki_test_SOURCES = test/nagamochi_ibaraki_test.cc
    9298test_path_test_SOURCES = test/path_test.cc
    9399test_planarity_test_SOURCES = test/planarity_test.cc
  • test/graph_copy_test.cc

    r984 r989  
    1919#include <lemon/smart_graph.h>
    2020#include <lemon/list_graph.h>
     21#include <lemon/static_graph.h>
    2122#include <lemon/lgf_reader.h>
    2223#include <lemon/error.h>
     
    2728using namespace lemon;
    2829
     30template <typename GR>
    2931void digraph_copy_test() {
    3032  const int nn = 10;
     
    5254    }
    5355  }
    54 
     56 
    5557  // Test digraph copy
    56   ListDigraph to;
    57   ListDigraph::NodeMap<int> tnm(to);
    58   ListDigraph::ArcMap<int> tam(to);
    59   ListDigraph::Node tn;
    60   ListDigraph::Arc ta;
    61 
    62   SmartDigraph::NodeMap<ListDigraph::Node> nr(from);
    63   SmartDigraph::ArcMap<ListDigraph::Arc> er(from);
    64 
    65   ListDigraph::NodeMap<SmartDigraph::Node> ncr(to);
    66   ListDigraph::ArcMap<SmartDigraph::Arc> ecr(to);
     58  GR to;
     59  typename GR::template NodeMap<int> tnm(to);
     60  typename GR::template ArcMap<int> tam(to);
     61  typename GR::Node tn;
     62  typename GR::Arc ta;
     63
     64  SmartDigraph::NodeMap<typename GR::Node> nr(from);
     65  SmartDigraph::ArcMap<typename GR::Arc> er(from);
     66
     67  typename GR::template NodeMap<SmartDigraph::Node> ncr(to);
     68  typename GR::template ArcMap<SmartDigraph::Arc> ecr(to);
    6769
    6870  digraphCopy(from, to).
     
    8789  }
    8890
    89   for (ListDigraph::NodeIt it(to); it != INVALID; ++it) {
     91  for (typename GR::NodeIt it(to); it != INVALID; ++it) {
    9092    check(nr[ncr[it]] == it, "Wrong copy.");
    9193  }
    9294
    93   for (ListDigraph::ArcIt it(to); it != INVALID; ++it) {
     95  for (typename GR::ArcIt it(to); it != INVALID; ++it) {
    9496    check(er[ecr[it]] == it, "Wrong copy.");
    9597  }
     
    104106}
    105107
     108template <typename GR>
    106109void graph_copy_test() {
    107110  const int nn = 10;
     
    136139
    137140  // Test graph copy
    138   ListGraph to;
    139   ListGraph::NodeMap<int> tnm(to);
    140   ListGraph::ArcMap<int> tam(to);
    141   ListGraph::EdgeMap<int> tem(to);
    142   ListGraph::Node tn;
    143   ListGraph::Arc ta;
    144   ListGraph::Edge te;
    145 
    146   SmartGraph::NodeMap<ListGraph::Node> nr(from);
    147   SmartGraph::ArcMap<ListGraph::Arc> ar(from);
    148   SmartGraph::EdgeMap<ListGraph::Edge> er(from);
    149 
    150   ListGraph::NodeMap<SmartGraph::Node> ncr(to);
    151   ListGraph::ArcMap<SmartGraph::Arc> acr(to);
    152   ListGraph::EdgeMap<SmartGraph::Edge> ecr(to);
     141  GR to;
     142  typename GR::template NodeMap<int> tnm(to);
     143  typename GR::template ArcMap<int> tam(to);
     144  typename GR::template EdgeMap<int> tem(to);
     145  typename GR::Node tn;
     146  typename GR::Arc ta;
     147  typename GR::Edge te;
     148
     149  SmartGraph::NodeMap<typename GR::Node> nr(from);
     150  SmartGraph::ArcMap<typename GR::Arc> ar(from);
     151  SmartGraph::EdgeMap<typename GR::Edge> er(from);
     152
     153  typename GR::template NodeMap<SmartGraph::Node> ncr(to);
     154  typename GR::template ArcMap<SmartGraph::Arc> acr(to);
     155  typename GR::template EdgeMap<SmartGraph::Edge> ecr(to);
    153156
    154157  graphCopy(from, to).
     
    185188  }
    186189
    187   for (ListGraph::NodeIt it(to); it != INVALID; ++it) {
     190  for (typename GR::NodeIt it(to); it != INVALID; ++it) {
    188191    check(nr[ncr[it]] == it, "Wrong copy.");
    189192  }
    190193
    191   for (ListGraph::ArcIt it(to); it != INVALID; ++it) {
     194  for (typename GR::ArcIt it(to); it != INVALID; ++it) {
    192195    check(ar[acr[it]] == it, "Wrong copy.");
    193196  }
    194   for (ListGraph::EdgeIt it(to); it != INVALID; ++it) {
     197  for (typename GR::EdgeIt it(to); it != INVALID; ++it) {
    195198    check(er[ecr[it]] == it, "Wrong copy.");
    196199  }
     
    209212
    210213int main() {
    211   digraph_copy_test();
    212   graph_copy_test();
     214  digraph_copy_test<SmartDigraph>();
     215  digraph_copy_test<ListDigraph>();
     216  digraph_copy_test<StaticDigraph>();
     217  graph_copy_test<SmartGraph>();
     218  graph_copy_test<ListGraph>();
    213219
    214220  return 0;
Note: See TracChangeset for help on using the changeset viewer.