COIN-OR::LEMON - Graph Library

Changes in / [971:22201ee8e437:970:bd523dbc7033] in lemon-main


Ignore:
Files:
7 deleted
19 edited

Legend:

Unmodified
Added
Removed
  • CMakeLists.txt

    r966 r964  
    115115SET(LEMON_HAVE_LONG_LONG ${HAVE_LONG_LONG})
    116116
     117INCLUDE(FindPythonInterp)
     118
    117119ENABLE_TESTING()
    118120
     
    125127ADD_SUBDIRECTORY(lemon)
    126128IF(${CMAKE_SOURCE_DIR} STREQUAL ${PROJECT_SOURCE_DIR})
    127   ADD_SUBDIRECTORY(contrib)
    128129  ADD_SUBDIRECTORY(demo)
    129130  ADD_SUBDIRECTORY(tools)
  • doc/Doxyfile.in

    r966 r964  
    9696                         "@abs_top_srcdir@/lemon/concepts" \
    9797                         "@abs_top_srcdir@/demo" \
    98                          "@abs_top_srcdir@/contrib" \
    9998                         "@abs_top_srcdir@/tools" \
    10099                         "@abs_top_srcdir@/test/test_tools.h" \
  • doc/coding_style.dox

    r919 r440  
    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_underscore
     104_start_with_underscores
    105105\endcode
    106106
  • doc/dirs.dox

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

    r919 r879  
    407407   strongly polynomial \ref klein67primal, \ref goldberg89cyclecanceling.
    408408
    409 In general, \ref NetworkSimplex and \ref CostScaling are the most efficient
    410 implementations, but the other two algorithms could be faster in special cases.
     409In general NetworkSimplex is the most efficient implementation,
     410but in special cases other algorithms could be faster.
    411411For example, if the total supply and/or capacities are rather small,
    412 \ref CapacityScaling is usually the fastest algorithm (without effective scaling).
     412CapacityScaling is usually the fastest algorithm (without effective scaling).
    413413*/
    414414
     
    472472  \ref dasdan98minmeancycle.
    473473
    474 In practice, the \ref HowardMmc "Howard" algorithm turned out to be by far the
     474In practice, the \ref HowardMmc "Howard" algorithm proved 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 Planar Embedding and Drawing
     542@defgroup planar Planarity Embedding and Drawing
    543543@ingroup algs
    544544\brief Algorithms for planarity checking, embedding and drawing
     
    552552
    553553/**
    554 @defgroup approx_algs Approximation Algorithms
     554@defgroup approx 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.
    564560*/
    565561
  • doc/references.bib

    r904 r755  
    298298  address =      {Dublin, Ireland},
    299299  year =         1991,
    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 }
     300  month =        sep,
     301}
  • lemon/Makefile.am

    r966 r964  
    9292        lemon/graph_to_eps.h \
    9393        lemon/grid_graph.h \
    94         lemon/grosso_locatelli_pullan_mc.h \
    9594        lemon/hartmann_orlin_mmc.h \
    9695        lemon/howard_mmc.h \
     
    109108        lemon/math.h \
    110109        lemon/min_cost_arborescence.h \
    111         lemon/max_cardinality_search.h \
    112         lemon/nagamochi_ibaraki.h \
    113110        lemon/nauty_reader.h \
    114111        lemon/network_simplex.h \
  • lemon/capacity_scaling.h

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

    r966 r964  
    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
    468449
    469450  /// \brief Class to copy a digraph.
  • lemon/cost_scaling.h

    r938 r931  
    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   ///
    103100  /// Most of the parameters of the problem (except for the digraph)
    104101  /// can be given using separate functions, and the algorithm can be
     
    117114  /// consider to use the named template parameters instead.
    118115  ///
    119   /// \warning Both \c V and \c C must be signed number types.
    120   /// \warning All input data (capacities, supply values, and costs) must
     116  /// \warning Both number types must be signed and all input data must
    121117  /// be integer.
    122   /// \warning This algorithm does not support negative costs for
    123   /// arcs having infinite upper bound.
     118  /// \warning This algorithm does not support negative costs for such
     119  /// arcs that have infinite upper bound.
    124120  ///
    125121  /// \note %CostScaling provides three different internal methods,
     
    183179    /// relabel operation.
    184180    /// By default, the so called \ref PARTIAL_AUGMENT
    185     /// "Partial Augment-Relabel" method is used, which turned out to be
     181    /// "Partial Augment-Relabel" method is used, which proved to be
    186182    /// the most efficient and the most robust on various test inputs.
    187183    /// However, the other methods can be selected using the \ref run()
     
    238234    };
    239235
     236    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
    240237    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
    241238
     
    288285    int _max_rank;
    289286
     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
    290295  public:
    291296
     
    334339    CostScaling(const GR& graph) :
    335340      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
     341      _cost_map(_cost_vec), _pi_map(_pi),
    336342      INF(std::numeric_limits<Value>::has_infinity ?
    337343          std::numeric_limits<Value>::infinity() :
     
    442448    ///
    443449    /// Using this function has the same effect as using \ref supplyMap()
    444     /// with a map in which \c k is assigned to \c s, \c -k is
     450    /// with such a map in which \c k is assigned to \c s, \c -k is
    445451    /// assigned to \c t and all other nodes have zero supply value.
    446452    ///
     
    488494    /// \param method The internal method that will be used in the
    489495    /// algorithm. For more information, see \ref Method.
    490     /// \param factor The cost scaling factor. It must be at least two.
     496    /// \param factor The cost scaling factor. It must be larger than one.
    491497    ///
    492498    /// \return \c INFEASIBLE if no feasible flow exists,
     
    502508    /// \see ProblemType, Method
    503509    /// \see resetParams(), reset()
    504     ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 16) {
    505       LEMON_ASSERT(factor >= 2, "The scaling factor must be at least 2");
     510    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
    506511      _alpha = factor;
    507512      ProblemType pt = init();
     
    567572    }
    568573
    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     ///
     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.
    585585    /// \return <tt>(*this)</tt>
    586     ///
    587     /// \see resetParams(), run()
    588586    CostScaling& reset() {
    589587      // Resize vectors
     
    610608      _excess.resize(_res_node_num);
    611609      _next_out.resize(_res_node_num);
     610
     611      _arc_vec.reserve(_res_arc_num);
     612      _cost_vec.reserve(_res_arc_num);
    612613
    613614      // Copy the graph
     
    886887      }
    887888
     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
    888897      // Initialize data structures for buckets
    889898      _max_rank = _alpha * _res_node_num;
     
    893902      _rank.resize(_res_node_num + 1);
    894903
    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 
     904      // Execute the algorithm
    902905      switch (method) {
    903906        case PUSH:
     
    908911          break;
    909912        case PARTIAL_AUGMENT:
    910           startAugment(MAX_PARTIAL_PATH_LENGTH);
     913          startAugment(MAX_PATH_LENGTH);
    911914          break;
    912915      }
    913916
    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       }
     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();
    971933
    972934      // Handle non-zero lower bounds
     
    986948        LargeCost pi_u = _pi[u];
    987949        for (int a = _first_out[u]; a != last_out; ++a) {
    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             }
     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;
    997957          }
    998958        }
     
    1010970    }
    1011971
    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);
     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;
    1235996    }
    1236997
    1237998    // Global potential update heuristic
    1238999    void globalUpdate() {
    1239       const int bucket_end = _root + 1;
     1000      int bucket_end = _root + 1;
    12401001
    12411002      // Initialize buckets
     
    12441005      }
    12451006      Value total_excess = 0;
    1246       int b0 = bucket_end;
    12471007      for (int i = 0; i != _res_node_num; ++i) {
    12481008        if (_excess[i] < 0) {
    12491009          _rank[i] = 0;
    1250           _bucket_next[i] = b0;
    1251           _bucket_prev[b0] = i;
    1252           b0 = i;
     1010          _bucket_next[i] = _buckets[0];
     1011          _bucket_prev[_buckets[0]] = i;
     1012          _buckets[0] = i;
    12531013        } else {
    12541014          total_excess += _excess[i];
     
    12571017      }
    12581018      if (total_excess == 0) return;
    1259       _buckets[0] = b0;
    12601019
    12611020      // Search the buckets
     
    12791038                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
    12801039                int new_rank_v = old_rank_v;
    1281                 if (nrc < LargeCost(_max_rank)) {
    1282                   new_rank_v = r + 1 + static_cast<int>(nrc);
    1283                 }
     1040                if (nrc < LargeCost(_max_rank))
     1041                  new_rank_v = r + 1 + int(nrc);
    12841042
    12851043                // Change the rank of v
     
    12931051                      _buckets[old_rank_v] = _bucket_next[v];
    12941052                    } else {
    1295                       int pv = _bucket_prev[v], nv = _bucket_next[v];
    1296                       _bucket_next[pv] = nv;
    1297                       _bucket_prev[nv] = pv;
     1053                      _bucket_next[_bucket_prev[v]] = _bucket_next[v];
     1054                      _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
    12981055                    }
    12991056                  }
    13001057
    1301                   // Insert v into its new bucket
    1302                   int nv = _buckets[new_rank_v];
    1303                   _bucket_next[v] = nv;
    1304                   _bucket_prev[nv] = v;
     1058                  // Insert v to its new bucket
     1059                  _bucket_next[v] = _buckets[new_rank_v];
     1060                  _bucket_prev[_buckets[new_rank_v]] = v;
    13051061                  _buckets[new_rank_v] = v;
    13061062                }
     
    13311087    void startAugment(int max_length) {
    13321088      // Paramters for heuristics
    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 *
     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 *
    13361093        (_res_node_num + _sup_node_num * _sup_node_num));
    1337       int next_global_update_limit = global_update_skip;
     1094      int next_update_limit = global_update_freq;
     1095
     1096      int relabel_cnt = 0;
    13381097
    13391098      // Perform cost scaling phases
    1340       IntVector path;
    1341       BoolVector path_arc(_res_arc_num, false);
    1342       int relabel_cnt = 0;
    1343       int eps_phase_cnt = 0;
     1099      std::vector<int> path;
    13441100      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    13451101                                        1 : _epsilon / _alpha )
    13461102      {
    1347         ++eps_phase_cnt;
    1348 
    1349         // Price refinement heuristic
    1350         if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
    1351           if (priceRefinement()) continue;
     1103        // Early termination heuristic
     1104        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
     1105          if (earlyTermination()) break;
    13521106        }
    13531107
     
    13661120
    13671121          // Find an augmenting path from the start node
     1122          path.clear();
    13681123          int tip = start;
    1369           while (int(path.size()) < max_length && _excess[tip] >= 0) {
     1124          while (_excess[tip] >= 0 && int(path.size()) < max_length) {
    13701125            int u;
    1371             LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
    1372             LargeCost pi_tip = _pi[tip];
     1126            LargeCost min_red_cost, rc, pi_tip = _pi[tip];
    13731127            int last_out = _first_out[tip+1];
    13741128            for (int a = _next_out[tip]; a != last_out; ++a) {
    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                 }
     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;
    13911135              }
    13921136            }
    13931137
    13941138            // Relabel tip node
     1139            min_red_cost = std::numeric_limits<LargeCost>::max();
    13951140            if (tip != start) {
    13961141              int ra = _reverse[path.back()];
    1397               min_red_cost =
    1398                 std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
     1142              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
    13991143            }
    1400             last_out = _next_out[tip];
    14011144            for (int a = _first_out[tip]; a != last_out; ++a) {
    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                 }
     1145              rc = _cost[a] + pi_tip - _pi[_target[a]];
     1146              if (_res_cap[a] > 0 && rc < min_red_cost) {
     1147                min_red_cost = rc;
    14071148              }
    14081149            }
     
    14131154            // Step back
    14141155            if (tip != start) {
    1415               int pa = path.back();
    1416               path_arc[pa] = false;
    1417               tip = _source[pa];
     1156              tip = _source[path.back()];
    14181157              path.pop_back();
    14191158            }
     
    14231162
    14241163          // Augment along the found path (as much flow as possible)
    1425         augment:
    14261164          Value delta;
    14271165          int pa, u, v = start;
     
    14301168            u = v;
    14311169            v = _target[pa];
    1432             path_arc[pa] = false;
    14331170            delta = std::min(_res_cap[pa], _excess[u]);
    14341171            _res_cap[pa] -= delta;
     
    14361173            _excess[u] -= delta;
    14371174            _excess[v] += delta;
    1438             if (_excess[v] > 0 && _excess[v] <= delta) {
     1175            if (_excess[v] > 0 && _excess[v] <= delta)
    14391176              _active_nodes.push_back(v);
    1440             }
    14411177          }
    1442           path.clear();
    14431178
    14441179          // Global update heuristic
    1445           if (relabel_cnt >= next_global_update_limit) {
     1180          if (relabel_cnt >= next_update_limit) {
    14461181            globalUpdate();
    1447             next_global_update_limit += global_update_skip;
     1182            next_update_limit += global_update_freq;
    14481183          }
    14491184        }
    1450 
    1451       }
    1452 
     1185      }
    14531186    }
    14541187
     
    14561189    void startPush() {
    14571190      // Paramters for heuristics
    1458       const int PRICE_REFINEMENT_LIMIT = 2;
     1191      const int EARLY_TERM_EPSILON_LIMIT = 1000;
    14591192      const double GLOBAL_UPDATE_FACTOR = 2.0;
    14601193
    1461       const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
     1194      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
    14621195        (_res_node_num + _sup_node_num * _sup_node_num));
    1463       int next_global_update_limit = global_update_skip;
     1196      int next_update_limit = global_update_freq;
     1197
     1198      int relabel_cnt = 0;
    14641199
    14651200      // Perform cost scaling phases
    14661201      BoolVector hyper(_res_node_num, false);
    14671202      LargeCostVector hyper_cost(_res_node_num);
    1468       int relabel_cnt = 0;
    1469       int eps_phase_cnt = 0;
    14701203      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    14711204                                        1 : _epsilon / _alpha )
    14721205      {
    1473         ++eps_phase_cnt;
    1474 
    1475         // Price refinement heuristic
    1476         if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
    1477           if (priceRefinement()) continue;
     1206        // Early termination heuristic
     1207        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
     1208          if (earlyTermination()) break;
    14781209        }
    14791210
     
    15471278               std::numeric_limits<LargeCost>::max();
    15481279            for (int a = _first_out[n]; a != last_out; ++a) {
    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                 }
     1280              rc = _cost[a] + pi_n - _pi[_target[a]];
     1281              if (_res_cap[a] > 0 && rc < min_red_cost) {
     1282                min_red_cost = rc;
    15541283              }
    15551284            }
     
    15691298
    15701299          // Global update heuristic
    1571           if (relabel_cnt >= next_global_update_limit) {
     1300          if (relabel_cnt >= next_update_limit) {
    15721301            globalUpdate();
    15731302            for (int u = 0; u != _res_node_num; ++u)
    15741303              hyper[u] = false;
    1575             next_global_update_limit += global_update_skip;
     1304            next_update_limit += global_update_freq;
    15761305          }
    15771306        }
  • lemon/cycle_canceling.h

    r922 r877  
    6666  /// algorithm. By default, it is the same as \c V.
    6767  ///
    68   /// \warning Both \c V and \c C must be signed number types.
    69   /// \warning All input data (capacities, supply values, and costs) must
     68  /// \warning Both number types must be signed and all input data must
    7069  /// be integer.
    71   /// \warning This algorithm does not support negative costs for
    72   /// arcs having infinite upper bound.
     70  /// \warning This algorithm does not support negative costs for such
     71  /// arcs that have infinite upper bound.
    7372  ///
    7473  /// \note For more information about the three available methods,
     
    118117    /// \ref CycleCanceling provides three different cycle-canceling
    119118    /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
    120     /// is used, which is by far the most efficient and the most robust.
     119    /// is used, which proved to be the most efficient and the most robust
     120    /// on various test inputs.
    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 a map in which \c k is assigned to \c s, \c -k is
     352    /// with such 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

    r919 r877  
    3737  ///Euler tour iterator for digraphs.
    3838
    39   /// \ingroup graph_properties
     39  /// \ingroup graph_prop
    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

    r915 r877  
    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. A notable
    57   /// use of this algorithm is testing network reliability.
     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.
    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.
    917915    ///
    918916    /// \pre \ref init() must be called before using this function.
     
    927925    /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with
    928926    /// \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.
    931927    ///
    932928    /// \pre \ref init() must be called before using this function.
     
    938934    /// \brief Run the algorithm.
    939935    ///
    940     /// This function runs the algorithm. It chooses source node,
    941     /// then calls \ref init(), \ref calculateOut()
     936    /// This function runs the algorithm. It finds nodes \c source and
     937    /// \c target arbitrarily and then calls \ref init(), \ref calculateOut()
    942938    /// and \ref calculateIn().
    943939    void run() {
     
    949945    /// \brief Run the algorithm.
    950946    ///
    951     /// This function runs the algorithm. It calls \ref init(),
    952     /// \ref calculateOut() and \ref calculateIn() with the given
    953     /// source node.
     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().
    954950    void run(const Node& s) {
    955951      init(s);
     
    970966    /// \brief Return the value of the minimum cut.
    971967    ///
    972     /// This function returns the value of the best cut found by the
    973     /// previously called \ref run(), \ref calculateOut() or \ref
    974     /// calculateIn().
     968    /// This function returns the value of the minimum cut.
    975969    ///
    976970    /// \pre \ref run(), \ref calculateOut() or \ref calculateIn()
     
    983977    /// \brief Return a minimum cut.
    984978    ///
    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
     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
    992982    /// for the nodes of \f$ X \f$).
    993983    ///
  • lemon/kruskal.h

    r921 r584  
    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///
    3336
    3437namespace lemon {
  • lemon/network_simplex.h

    r922 r889  
    4848  /// flow problem.
    4949  ///
    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.
     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.
    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 \c V and \c C must be signed number types.
    67   /// \warning All input data (capacities, supply values, and costs) must
     66  /// \warning Both number types must be signed and all input data must
    6867  /// be integer.
    6968  ///
     
    127126    /// of the algorithm.
    128127    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
    129     /// turend out to be the most efficient and the most robust on various
     128    /// proved to be the most efficient and the most robust on various
    130129    /// test inputs.
    131130    /// However, another pivot rule can be selected using the \ref run()
     
    168167    typedef std::vector<Value> ValueVector;
    169168    typedef std::vector<Cost> CostVector;
    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
     169    typedef std::vector<char> BoolVector;
     170    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    173171
    174172    // State constants for arcs
     
    179177    };
    180178
    181     // Direction constants for tree arcs
    182     enum ArcDirection {
    183       DIR_DOWN = -1,
    184       DIR_UP   =  1
    185     };
     179    typedef std::vector<signed char> StateVector;
     180    // Note: vector<signed char> is used instead of vector<ArcState> for
     181    // efficiency reasons
    186182
    187183  private:
     
    222218    IntVector _succ_num;
    223219    IntVector _last_succ;
    224     CharVector _pred_dir;
    225     CharVector _state;
    226220    IntVector _dirty_revs;
     221    BoolVector _forward;
     222    StateVector _state;
    227223    int _root;
    228224
    229225    // Temporary data used in the current pivot iteration
    230226    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;
    231229    Value delta;
    232230
     
    253251      const IntVector  &_target;
    254252      const CostVector &_cost;
    255       const CharVector &_state;
     253      const StateVector &_state;
    256254      const CostVector &_pi;
    257255      int &_in_arc;
     
    305303      const IntVector  &_target;
    306304      const CostVector &_cost;
    307       const CharVector &_state;
     305      const StateVector &_state;
    308306      const CostVector &_pi;
    309307      int &_in_arc;
     
    344342      const IntVector  &_target;
    345343      const CostVector &_cost;
    346       const CharVector &_state;
     344      const StateVector &_state;
    347345      const CostVector &_pi;
    348346      int &_in_arc;
     
    417415      const IntVector  &_target;
    418416      const CostVector &_cost;
    419       const CharVector &_state;
     417      const StateVector &_state;
    420418      const CostVector &_pi;
    421419      int &_in_arc;
     
    520518      const IntVector  &_target;
    521519      const CostVector &_cost;
    522       const CharVector &_state;
     520      const StateVector &_state;
    523521      const CostVector &_pi;
    524522      int &_in_arc;
     
    573571        // Check the current candidate list
    574572        int e;
    575         Cost c;
    576573        for (int i = 0; i != _curr_length; ++i) {
    577574          e = _candidates[i];
    578           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    579           if (c < 0) {
    580             _cand_cost[e] = c;
    581           } else {
     575          _cand_cost[e] = _state[e] *
     576            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     577          if (_cand_cost[e] >= 0) {
    582578            _candidates[i--] = _candidates[--_curr_length];
    583579          }
     
    589585
    590586        for (e = _next_arc; e != _search_arc_num; ++e) {
    591           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    592           if (c < 0) {
    593             _cand_cost[e] = c;
     587          _cand_cost[e] = _state[e] *
     588            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     589          if (_cand_cost[e] < 0) {
    594590            _candidates[_curr_length++] = e;
    595591          }
     
    638634    ///
    639635    /// \param graph The digraph the algorithm runs on.
    640     /// \param arc_mixing Indicate if the arcs will be stored in a
     636    /// \param arc_mixing Indicate if the arcs have to be stored in a
    641637    /// mixed order in the internal data structure.
    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) :
     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) :
    646641      _graph(graph), _node_id(graph), _arc_id(graph),
    647642      _arc_mixing(arc_mixing),
     
    736731    ///
    737732    /// \return <tt>(*this)</tt>
    738     ///
    739     /// \sa supplyType()
    740733    template<typename SupplyMap>
    741734    NetworkSimplex& supplyMap(const SupplyMap& map) {
     
    754747    ///
    755748    /// Using this function has the same effect as using \ref supplyMap()
    756     /// with a map in which \c k is assigned to \c s, \c -k is
     749    /// with such a map in which \c k is assigned to \c s, \c -k is
    757750    /// assigned to \c t and all other nodes have zero supply value.
    758751    ///
     
    921914      _parent.resize(all_node_num);
    922915      _pred.resize(all_node_num);
    923       _pred_dir.resize(all_node_num);
     916      _forward.resize(all_node_num);
    924917      _thread.resize(all_node_num);
    925918      _rev_thread.resize(all_node_num);
     
    935928      if (_arc_mixing) {
    936929        // Store the arcs in a mixed order
    937         const int skip = std::max(_arc_num / _node_num, 3);
     930        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
    938931        int i = 0, j = 0;
    939932        for (ArcIt a(_graph); a != INVALID; ++a) {
     
    941934          _source[i] = _node_id[_graph.source(a)];
    942935          _target[i] = _node_id[_graph.target(a)];
    943           if ((i += skip) >= _arc_num) i = ++j;
     936          if ((i += k) >= _arc_num) i = ++j;
    944937        }
    945938      } else {
     
    11241117          _state[e] = STATE_TREE;
    11251118          if (_supply[u] >= 0) {
    1126             _pred_dir[u] = DIR_UP;
     1119            _forward[u] = true;
    11271120            _pi[u] = 0;
    11281121            _source[e] = u;
     
    11311124            _cost[e] = 0;
    11321125          } else {
    1133             _pred_dir[u] = DIR_DOWN;
     1126            _forward[u] = false;
    11341127            _pi[u] = ART_COST;
    11351128            _source[e] = _root;
     
    11511144          _last_succ[u] = u;
    11521145          if (_supply[u] >= 0) {
    1153             _pred_dir[u] = DIR_UP;
     1146            _forward[u] = true;
    11541147            _pi[u] = 0;
    11551148            _pred[u] = e;
     
    11611154            _state[e] = STATE_TREE;
    11621155          } else {
    1163             _pred_dir[u] = DIR_DOWN;
     1156            _forward[u] = false;
    11641157            _pi[u] = ART_COST;
    11651158            _pred[u] = f;
     
    11921185          _last_succ[u] = u;
    11931186          if (_supply[u] <= 0) {
    1194             _pred_dir[u] = DIR_DOWN;
     1187            _forward[u] = false;
    11951188            _pi[u] = 0;
    11961189            _pred[u] = e;
     
    12021195            _state[e] = STATE_TREE;
    12031196          } else {
    1204             _pred_dir[u] = DIR_UP;
     1197            _forward[u] = true;
    12051198            _pi[u] = -ART_COST;
    12061199            _pred[u] = f;
     
    12451238      // Initialize first and second nodes according to the direction
    12461239      // of the cycle
    1247       int first, second;
    12481240      if (_state[in_arc] == STATE_LOWER) {
    12491241        first  = _source[in_arc];
     
    12551247      delta = _cap[in_arc];
    12561248      int result = 0;
    1257       Value c, d;
     1249      Value d;
    12581250      int e;
    12591251
    1260       // Search the cycle form the first node to the join node
     1252      // Search the cycle along the path form the first node to the root
    12611253      for (int u = first; u != join; u = _parent[u]) {
    12621254        e = _pred[u];
    1263         d = _flow[e];
    1264         if (_pred_dir[u] == DIR_DOWN) {
    1265           c = _cap[e];
    1266           d = c >= MAX ? INF : c - d;
    1267         }
     1255        d = _forward[u] ?
     1256          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
    12681257        if (d < delta) {
    12691258          delta = d;
     
    12721261        }
    12731262      }
    1274 
    1275       // Search the cycle form the second node to the join node
     1263      // Search the cycle along the path form the second node to the root
    12761264      for (int u = second; u != join; u = _parent[u]) {
    12771265        e = _pred[u];
    1278         d = _flow[e];
    1279         if (_pred_dir[u] == DIR_UP) {
    1280           c = _cap[e];
    1281           d = c >= MAX ? INF : c - d;
    1282         }
     1266        d = _forward[u] ?
     1267          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
    12831268        if (d <= delta) {
    12841269          delta = d;
     
    13051290        _flow[in_arc] += val;
    13061291        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
    1307           _flow[_pred[u]] -= _pred_dir[u] * val;
     1292          _flow[_pred[u]] += _forward[u] ? -val : val;
    13081293        }
    13091294        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
    1310           _flow[_pred[u]] += _pred_dir[u] * val;
     1295          _flow[_pred[u]] += _forward[u] ? val : -val;
    13111296        }
    13121297      }
     
    13231308    // Update the tree structure
    13241309    void updateTreeStructure() {
     1310      int u, w;
    13251311      int old_rev_thread = _rev_thread[u_out];
    13261312      int old_succ_num = _succ_num[u_out];
     
    13281314      v_out = _parent[u_out];
    13291315
    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         }
     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]];
    13481323      } else {
    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;
     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;
    14161398      }
    14171399
    14181400      // Update _last_succ from v_in towards the root
    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 
     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      }
    14251405      // Update _last_succ from v_out towards the root
    14261406      if (join != old_rev_thread && v_in != old_rev_thread) {
    1427         for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
     1407        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
    14281408             u = _parent[u]) {
    14291409          _last_succ[u] = old_rev_thread;
    14301410        }
    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;
     1411      } else {
     1412        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
    14341413             u = _parent[u]) {
    1435           _last_succ[u] = last_succ_out;
     1414          _last_succ[u] = _last_succ[u_out];
    14361415        }
    14371416      }
    14381417
    14391418      // Update _succ_num from v_in to join
    1440       for (int u = v_in; u != join; u = _parent[u]) {
     1419      for (u = v_in; u != join; u = _parent[u]) {
    14411420        _succ_num[u] += old_succ_num;
    14421421      }
    14431422      // Update _succ_num from v_out to join
    1444       for (int u = v_out; u != join; u = _parent[u]) {
     1423      for (u = v_out; u != join; u = _parent[u]) {
    14451424        _succ_num[u] -= old_succ_num;
    14461425      }
    14471426    }
    14481427
    1449     // Update potentials in the subtree that has been moved
     1428    // Update potentials
    14501429    void updatePotential() {
    1451       Cost sigma = _pi[v_in] - _pi[u_in] -
    1452                    _pred_dir[u_in] * _cost[in_arc];
     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
    14531434      int end = _thread[_last_succ[u_in]];
    14541435      for (int u = u_in; u != end; u = _thread[u]) {
  • lemon/path.h

    r920 r877  
    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 n-th arc.
     138    /// \brief The nth 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 n-th arc
     146    /// \brief Initialize arc iterator to point to the nth 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 n-th arc.
     330    /// \brief The nth 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 n-th arc.
     337    /// \brief  Initializes arc iterator to point to the nth 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 n-th arc.
    508     ///
    509     /// This function looks for the n-th arc in O(n) time.
     507    /// \brief The nth arc.
     508    ///
     509    /// This function looks for the nth 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 n-th arc.
     519    /// \brief Initializes arc iterator to point to the nth 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 n-th arc.
     834    /// \brief The nth 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 n-th arc.
     841    /// \brief The arc iterator pointing to the nth 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

    r966 r964  
    3737  maps_test
    3838  matching_test
    39   max_cardinality_search_test
    40   max_clique_test
    4139  min_cost_arborescence_test
    4240  min_cost_flow_test
    4341  min_mean_cycle_test
    44   nagamochi_ibaraki_test
    4542  path_test
    4643  planarity_test
  • test/Makefile.am

    r966 r964  
    3535        test/maps_test \
    3636        test/matching_test \
    37         test/max_cardinality_search_test \
    38         test/max_clique_test \
    3937        test/min_cost_arborescence_test \
    4038        test/min_cost_flow_test \
    4139        test/min_mean_cycle_test \
    42         test/nagamochi_ibaraki_test \
    4340        test/path_test \
    4441        test/planarity_test \
     
    8279test_graph_test_SOURCES = test/graph_test.cc
    8380test_graph_utils_test_SOURCES = test/graph_utils_test.cc
    84 test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
    8581test_heap_test_SOURCES = test/heap_test.cc
    8682test_kruskal_test_SOURCES = test/kruskal_test.cc
     83test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
    8784test_lgf_test_SOURCES = test/lgf_test.cc
    8885test_lp_test_SOURCES = test/lp_test.cc
     
    9087test_mip_test_SOURCES = test/mip_test.cc
    9188test_matching_test_SOURCES = test/matching_test.cc
    92 test_max_cardinality_search_test_SOURCES = test/max_cardinality_search_test.cc
    93 test_max_clique_test_SOURCES = test/max_clique_test.cc
    9489test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
    9590test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
    9691test_min_mean_cycle_test_SOURCES = test/min_mean_cycle_test.cc
    97 test_nagamochi_ibaraki_test_SOURCES = test/nagamochi_ibaraki_test.cc
    9892test_path_test_SOURCES = test/path_test.cc
    9993test_planarity_test_SOURCES = test/planarity_test.cc
  • test/graph_copy_test.cc

    r894 r893  
    1919#include <lemon/smart_graph.h>
    2020#include <lemon/list_graph.h>
    21 #include <lemon/static_graph.h>
    2221#include <lemon/lgf_reader.h>
    2322#include <lemon/error.h>
     
    2827using namespace lemon;
    2928
    30 template <typename GR>
    3129void digraph_copy_test() {
    3230  const int nn = 10;
     
    5452    }
    5553  }
    56  
     54
    5755  // Test digraph copy
    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);
     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);
    6967
    7068  digraphCopy(from, to).
     
    8987  }
    9088
    91   for (typename GR::NodeIt it(to); it != INVALID; ++it) {
     89  for (ListDigraph::NodeIt it(to); it != INVALID; ++it) {
    9290    check(nr[ncr[it]] == it, "Wrong copy.");
    9391  }
    9492
    95   for (typename GR::ArcIt it(to); it != INVALID; ++it) {
     93  for (ListDigraph::ArcIt it(to); it != INVALID; ++it) {
    9694    check(er[ecr[it]] == it, "Wrong copy.");
    9795  }
     
    106104}
    107105
    108 template <typename GR>
    109106void graph_copy_test() {
    110107  const int nn = 10;
     
    139136
    140137  // Test graph copy
    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);
     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);
    156153
    157154  graphCopy(from, to).
     
    188185  }
    189186
    190   for (typename GR::NodeIt it(to); it != INVALID; ++it) {
     187  for (ListGraph::NodeIt it(to); it != INVALID; ++it) {
    191188    check(nr[ncr[it]] == it, "Wrong copy.");
    192189  }
    193190
    194   for (typename GR::ArcIt it(to); it != INVALID; ++it) {
     191  for (ListGraph::ArcIt it(to); it != INVALID; ++it) {
    195192    check(ar[acr[it]] == it, "Wrong copy.");
    196193  }
    197   for (typename GR::EdgeIt it(to); it != INVALID; ++it) {
     194  for (ListGraph::EdgeIt it(to); it != INVALID; ++it) {
    198195    check(er[ecr[it]] == it, "Wrong copy.");
    199196  }
     
    212209
    213210int main() {
    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>();
     211  digraph_copy_test();
     212  graph_copy_test();
    219213
    220214  return 0;
Note: See TracChangeset for help on using the changeset viewer.