lemon-0.x: comparison lemon/nagamochi

equal deleted inserted replaced

-:7b08652f0c58
+:17fec95e6617
 * express or implied, and with no claim as to its suitability for any
 * purpose.
 *
 */
-#ifndef LEMON_MIN_CUT_H
+#ifndef LEMON_NAGAMOCHI_IBARAKI_H
-#define LEMON_MIN_CUT_H
+#define LEMON_NAGAMOCHI_IBARAKI_H
 /// \ingroup topology
 /// \file
-/// \brief Maximum cardinality search and min cut in undirected graphs.
+/// \brief Maximum cardinality search and minimum cut in undirected
+/// graphs.
 #include <lemon/list_graph.h>
 #include <lemon/bin_heap.h>
 #include <lemon/bucket_heap.h>
+#include <lemon/unionfind.h>
+#include <lemon/topology.h>
 #include <lemon/bits/invalid.h>
 #include <lemon/error.h>
 #include <lemon/maps.h>
 #include <functional>
+#include <lemon/graph_writer.h>
+#include <lemon/time_measure.h>
 namespace lemon {
 namespace _min_cut_bits {
 /// \param graph is the graph, to which we would like to define the \ref
 /// CardinalityMap
 static CardinalityMap *createCardinalityMap(const Graph &graph) {
 return new CardinalityMap(graph);
 }
 };
 /// \ingroup topology
 ///
 /// \pre \ref run() must be called before using this function.
 /// \warning If node \c v in unreachable from the root the return value
 /// of this funcion is undefined.
 Value cardinality(Node node) const { return (*_cardinality)[node]; }
+/// \brief The current cardinality of a node.
+///
+/// Returns the current cardinality of a node.
+/// \pre the given node should be reached but not processed
+Value currentCardinality(Node node) const { return (*_heap)[node]; }
 /// \brief Returns a reference to the NodeMap of cardinalities.
 ///
 /// Returns a reference to the NodeMap of cardinalities. \pre \ref run()
 /// must be called before using this function.
 const CardinalityMap &cardinalityMap() const { return *_cardinality;}
 #endif
 {
 throw UninitializedParameter();
 }
+/// \brief The CutValueMap type
+///
+/// The type of the map that stores the cut value of a node.
+typedef AuxGraph::NodeMap<Value> AuxCutValueMap;
+/// \brief Instantiates a AuxCutValueMap.
+///
+/// This function instantiates a \ref AuxCutValueMap.
+static AuxCutValueMap *createAuxCutValueMap(const AuxGraph& graph) {
+return new AuxCutValueMap(graph);
+}
 /// \brief The AuxCapacityMap type
 ///
-/// The type of the map that stores the auxing edge capacities.
+/// The type of the map that stores the auxiliary edge capacities.
 typedef AuxGraph::UEdgeMap<Value> AuxCapacityMap;
 /// \brief Instantiates a AuxCapacityMap.
 ///
 /// This function instantiates a \ref AuxCapacityMap.
 return new ListRefMap(graph);
 }
 };
-namespace _min_cut_bits {
-template <typename _Key>
-class LastTwoMap {
-public:
-typedef _Key Key;
-typedef bool Value;
-LastTwoMap(int _num) : num(_num) {}
-void set(const Key& key, bool val) {
-if (!val) return;
---num;
-if (num > 1) return;
-keys[num] = key;
-}
-Key operator[](int index) const { return keys[index]; }
-private:
-Key keys[2];
-int num;
-};
-}
 /// \ingroup topology
 ///
 /// \brief Calculates the minimum cut in an undirected graph.
 ///
 /// the network reliability specifically to test how many links have
 /// to be destroyed in the network to split it at least two
 /// distinict subnetwork.
 ///
 /// The complexity of the algorithm is \f$ O(ne\log(n)) \f$ but with
-/// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n))
+/// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n)) \f$.
-/// \f$. When capacity map is neutral then it uses BucketHeap which
+/// When capacity map is neutral then it uses BucketHeap which
 /// results \f$ O(ne) \f$ time complexity.
+///
+/// \warning The value type of the capacity map should be able to hold
+/// any cut value of the graph, otherwise the result can overflow.
 #ifdef DOXYGEN
 template <typename _Graph, typename _CapacityMap, typename _Traits>
 #else
 template <typename _Graph = ListUGraph,
 	    typename _CapacityMap = typename _Graph::template UEdgeMap<int>,
 typedef typename Traits::CapacityMap CapacityMap;
 /// The type of the aux graph
 typedef typename Traits::AuxGraph AuxGraph;
 /// The type of the aux capacity map
 typedef typename Traits::AuxCapacityMap AuxCapacityMap;
+/// The type of the aux cut value map
+typedef typename Traits::AuxCutValueMap AuxCutValueMap;
 /// The cross reference type used for the current heap.
 typedef typename Traits::HeapCrossRef HeapCrossRef;
 /// The heap type used by the max cardinality algorithm.
 typedef typename Traits::Heap Heap;
 /// The node refrefernces between the original and aux graph type.
 /// Pointer to the aux capacity map
 AuxCapacityMap *_aux_capacity;
 /// \brief Indicates if \ref _aux_capacity is locally allocated
 /// (\c true) or not.
 bool local_aux_capacity;
+/// Pointer to the aux cut value map
+AuxCutValueMap *_aux_cut_value;
+/// \brief Indicates if \ref _aux_cut_value is locally allocated
+/// (\c true) or not.
+bool local_aux_cut_value;
 /// Pointer to the heap cross references.
 HeapCrossRef *_heap_cross_ref;
 /// \brief Indicates if \ref _heap_cross_ref is locally allocated
 /// (\c true) or not.
 bool local_heap_cross_ref;
 /// The min cut value.
 Value _min_cut;
 /// The number of the nodes of the aux graph.
 int _node_num;
-/// The first and last node of the min cut in the next list;
+/// The first and last node of the min cut in the next list.
-typename Graph::Node _first_node, _last_node;
+std::vector<typename Graph::Node> _cut;
 /// \brief The first and last element in the list associated
 /// to the aux graph node.
 NodeRefMap *_first, *_last;
 /// \brief The next node in the node lists.
 ListRefMap *_next;
-void create_structures() {
+void createStructures() {
 if (!_capacity) {
 local_capacity = true;
 _capacity = Traits::createCapacityMap(*_graph);
 }
 if(!_aux_graph) {
 }
 if(!_aux_capacity) {
 	local_aux_capacity = true;
 	_aux_capacity = Traits::createAuxCapacityMap(*_aux_graph);
 }
+if(!_aux_cut_value) {
+	local_aux_cut_value = true;
+	_aux_cut_value = Traits::createAuxCutValueMap(*_aux_graph);
+}
 _first = Traits::createNodeRefMap(*_aux_graph);
 _last = Traits::createNodeRefMap(*_aux_graph);
 _next = Traits::createListRefMap(*_graph);
-typename Graph::template NodeMap<typename AuxGraph::Node> ref(*_graph);
-for (typename Graph::NodeIt it(*_graph); it != INVALID; ++it) {
-_next->set(it, INVALID);
-typename AuxGraph::Node node = _aux_graph->addNode();
-_first->set(node, it);
-_last->set(node, it);
-ref.set(it, node);
-}
-for (typename Graph::UEdgeIt it(*_graph); it != INVALID; ++it) {
-if (_graph->source(it) == _graph->target(it)) continue;
-typename AuxGraph::UEdge uedge =
-_aux_graph->addEdge(ref[_graph->source(it)],
-ref[_graph->target(it)]);
-_aux_capacity->set(uedge, (*_capacity)[it]);
-}
 if (!_heap_cross_ref) {
 	local_heap_cross_ref = true;
 	_heap_cross_ref = Traits::createHeapCrossRef(*_aux_graph);
 }
 if (!_heap) {
 	local_heap = true;
 	_heap = Traits::createHeap(*_heap_cross_ref);
 }
 }
+void createAuxGraph() {
+typename Graph::template NodeMap<typename AuxGraph::Node> ref(*_graph);
+for (typename Graph::NodeIt n(*_graph); n != INVALID; ++n) {
+_next->set(n, INVALID);
+typename AuxGraph::Node node = _aux_graph->addNode();
+_first->set(node, n);
+_last->set(node, n);
+ref.set(n, node);
+}
+typename AuxGraph::template NodeMap<typename AuxGraph::UEdge>
+edges(*_aux_graph, INVALID);
+for (typename Graph::NodeIt n(*_graph); n != INVALID; ++n) {
+for (typename Graph::IncEdgeIt e(*_graph, n); e != INVALID; ++e) {
+typename Graph::Node tn = _graph->runningNode(e);
+if (n < tn || n == tn) continue;
+if (edges[ref[tn]] != INVALID) {
+Value value =
+(*_aux_capacity)[edges[ref[tn]]] + (*_capacity)[e];
+_aux_capacity->set(edges[ref[tn]], value);
+} else {
+edges.set(ref[tn], _aux_graph->addEdge(ref[n], ref[tn]));
+Value value = (*_capacity)[e];
+_aux_capacity->set(edges[ref[tn]], value);
+}
+}
+for (typename Graph::IncEdgeIt e(*_graph, n); e != INVALID; ++e) {
+typename Graph::Node tn = _graph->runningNode(e);
+edges.set(ref[tn], INVALID);
+}
+}
+_cut.resize(1, INVALID);
+_min_cut = std::numeric_limits<Value>::max();
+for (typename AuxGraph::NodeIt n(*_aux_graph); n != INVALID; ++n) {
+Value value = 0;
+for (typename AuxGraph::IncEdgeIt e(*_aux_graph, n);
+e != INVALID; ++e) {
+value += (*_aux_capacity)[e];
+}
+if (_min_cut > value) {
+_min_cut = value;
+_cut[0] = (*_first)[n];
+}
+(*_aux_cut_value)[n] = value;
+}
+}
 public :
 typedef NagamochiIbaraki Create;
 NagamochiIbaraki(const Graph& graph, const CapacityMap& capacity)
 : _graph(&graph),
 _capacity(&capacity), local_capacity(false),
 _aux_graph(0), local_aux_graph(false),
 _aux_capacity(0), local_aux_capacity(false),
+_aux_cut_value(0), local_aux_cut_value(false),
 _heap_cross_ref(0), local_heap_cross_ref(false),
 _heap(0), local_heap(false),
 _first(0), _last(0), _next(0) {}
 /// \brief Constructor.
 NagamochiIbaraki(const Graph& graph)
 : _graph(&graph),
 _capacity(0), local_capacity(false),
 _aux_graph(0), local_aux_graph(false),
 _aux_capacity(0), local_aux_capacity(false),
+_aux_cut_value(0), local_aux_cut_value(false),
 _heap_cross_ref(0), local_heap_cross_ref(false),
 _heap(0), local_heap(false),
 _first(0), _last(0), _next(0) {}
 /// \brief Destructor.
 if (local_heap_cross_ref) delete _heap_cross_ref;
 if (_first) delete _first;
 if (_last) delete _last;
 if (_next) delete _next;
 if (local_aux_capacity) delete _aux_capacity;
+if (local_aux_cut_value) delete _aux_cut_value;
 if (local_aux_graph) delete _aux_graph;
 if (local_capacity) delete _capacity;
 }
 /// \brief Sets the heap and the cross reference used by algorithm.
 /// \brief Initializes the internal data structures.
 ///
 /// Initializes the internal data structures.
 void init() {
-create_structures();
-_first_node = _last_node = INVALID;
 _node_num = countNodes(*_graph);
-}
+createStructures();
+createAuxGraph();
+}
+private:
+struct EdgeInfo {
+typename AuxGraph::Node source, target;
+Value capacity;
+};
+struct EdgeInfoLess {
+bool operator()(const EdgeInfo& left, const EdgeInfo& right) const {
+return (left.source < right.source) ||
+(left.source == right.source && left.target < right.target);
+}
+};
+public:
 /// \brief Processes the next phase
 ///
-/// Processes the next phase in the algorithm. The function
+/// Processes the next phase in the algorithm. The function should
-/// should be called countNodes(graph) - 1 times to get
+/// be called at most countNodes(graph) - 1 times to get surely
-/// surely the min cut in the graph. The
+/// the min cut in the graph.
 ///
 ///\return %True when the algorithm finished.
 bool processNextPhase() {
 if (_node_num <= 1) return true;
-using namespace _min_cut_bits;
 typedef typename AuxGraph::Node Node;
 typedef typename AuxGraph::NodeIt NodeIt;
 typedef typename AuxGraph::UEdge UEdge;
+typedef typename AuxGraph::UEdgeIt UEdgeIt;
 typedef typename AuxGraph::IncEdgeIt IncEdgeIt;
-typedef typename MaxCardinalitySearch<AuxGraph, AuxCapacityMap>::
+typename AuxGraph::template UEdgeMap<Value> _edge_cut(*_aux_graph);
-template DefHeap<Heap, HeapCrossRef>::
-template DefCardinalityMap<NullMap<Node, Value> >::
-template DefProcessedMap<LastTwoMap<Node> >::
+for (NodeIt n(*_aux_graph); n != INVALID; ++n) {
-Create MaxCardinalitySearch;
+_heap_cross_ref->set(n, Heap::PRE_HEAP);
+}
+std::vector<Node> nodes;
+nodes.reserve(_node_num);
+int sep = 0;
+Value alpha = 0;
+Value emc = std::numeric_limits<Value>::max();
+_heap->push(typename AuxGraph::NodeIt(*_aux_graph), 0);
+while (!_heap->empty()) {
+Node node = _heap->top();
+Value value = _heap->prio();
+_heap->pop();
+for (typename AuxGraph::IncEdgeIt e(*_aux_graph, node);
+e != INVALID; ++e) {
+Node tn = _aux_graph->runningNode(e);
+switch (_heap->state(tn)) {
+case Heap::PRE_HEAP:
+_heap->push(tn, (*_aux_capacity)[e]);
+_edge_cut[e] = (*_heap)[tn];
+break;
+case Heap::IN_HEAP:
+_heap->decrease(tn, (*_aux_capacity)[e] + (*_heap)[tn]);
+_edge_cut[e] = (*_heap)[tn];
+break;
+case Heap::POST_HEAP:
+break;
+}
+}
+alpha += (*_aux_cut_value)[node];
+alpha -= 2 * value;
+nodes.push_back(node);
+if (!_heap->empty()) {
+if (alpha < emc) {
+emc = alpha;
+sep = nodes.size();
+}
+}
+}
+if ((int)nodes.size() < _node_num) {
+_aux_graph->clear();
+_node_num = 1;
+_cut.clear();
+for (int i = 0; i < (int)nodes.size(); ++i) {
+typename Graph::Node n = (*_first)[nodes[i]];
+while (n != INVALID) {
+_cut.push_back(n);
+n = (*_next)[n];
+}
+}
+_min_cut = 0;
+return true;
+}
+if (emc < _min_cut) {
+_cut.clear();
+for (int i = 0; i < sep; ++i) {
+typename Graph::Node n = (*_first)[nodes[i]];
+while (n != INVALID) {
+_cut.push_back(n);
+n = (*_next)[n];
+}
+}
+_min_cut = emc;
+}
+typedef typename AuxGraph::template NodeMap<int> UfeCr;
+UfeCr ufecr(*_aux_graph);
+typedef UnionFindEnum<UfeCr> Ufe;
+Ufe ufe(ufecr);
+for (typename AuxGraph::NodeIt n(*_aux_graph); n != INVALID; ++n) {
+ufe.insert(n);
+}
+for (typename AuxGraph::UEdgeIt e(*_aux_graph); e != INVALID; ++e) {
+if (_edge_cut[e] >= emc) {
+ufe.join(_aux_graph->source(e), _aux_graph->target(e));
+}
+}
+if (ufe.size((typename Ufe::ClassIt)(ufe)) == _node_num) {
+_aux_graph->clear();
+_node_num = 1;
+return true;
+}
-MaxCardinalitySearch mcs(*_aux_graph, *_aux_capacity);
+std::vector<typename AuxGraph::Node> remnodes;
-for (NodeIt it(*_aux_graph); it != INVALID; ++it) {
-_heap_cross_ref->set(it, Heap::PRE_HEAP);
-}
-mcs.heap(*_heap, *_heap_cross_ref);
-LastTwoMap<Node> last_two_nodes(_node_num);
-mcs.processedMap(last_two_nodes);
-NullMap<Node, Value> cardinality;
-mcs.cardinalityMap(cardinality);
-mcs.run();
-Node new_node = _aux_graph->addNode();
 typename AuxGraph::template NodeMap<UEdge> edges(*_aux_graph, INVALID);
+for (typename Ufe::ClassIt c(ufe); c != INVALID; ++c) {
-Node first_node = last_two_nodes[0];
+if (ufe.size(c) == 1) continue;
-Node second_node = last_two_nodes[1];
+for (typename Ufe::ItemIt r(ufe, c); r != INVALID; ++r) {
+if ((Node)r == (Node)c) continue;
-_next->set((*_last)[first_node], (*_first)[second_node]);
+_next->set((*_last)[c], (*_first)[r]);
-_first->set(new_node, (*_first)[first_node]);
+_last->set(c, (*_last)[r]);
-_last->set(new_node, (*_last)[second_node]);
+remnodes.push_back(r);
+--_node_num;
-Value current_cut = 0;
+}
-for (IncEdgeIt it(*_aux_graph, first_node); it != INVALID; ++it) {
+}
-Node node = _aux_graph->runningNode(it);
-current_cut += (*_aux_capacity)[it];
+std::vector<EdgeInfo> addedges;
-if (node == second_node) continue;
+std::vector<UEdge> remedges;
-if (edges[node] == INVALID) {
-edges[node] = _aux_graph->addEdge(new_node, node);
+for (typename AuxGraph::UEdgeIt e(*_aux_graph);
-(*_aux_capacity)[edges[node]] = (*_aux_capacity)[it];
+e != INVALID; ++e) {
+Node sn = ufe.find(_aux_graph->source(e));
+Node tn = ufe.find(_aux_graph->target(e));
+if ((ufe.size(sn) == 1 && ufe.size(tn) == 1)) {
+continue;
+}
+if (sn == tn) {
+remedges.push_back(e);
+continue;
+}
+EdgeInfo info;
+if (sn < tn) {
+info.source = sn;
+info.target = tn;
 } else {
-(*_aux_capacity)[edges[node]] += (*_aux_capacity)[it];
+info.source = tn;
+info.target = sn;
 }
-}
+info.capacity = (*_aux_capacity)[e];
+addedges.push_back(info);
-if (_first_node == INVALID || current_cut < _min_cut) {
+remedges.push_back(e);
-_first_node = (*_first)[first_node];
+}
-_last_node = (*_last)[first_node];
-_min_cut = current_cut;
+for (int i = 0; i < (int)remedges.size(); ++i) {
-}
+_aux_graph->erase(remedges[i]);
+}
-_aux_graph->erase(first_node);
+sort(addedges.begin(), addedges.end(), EdgeInfoLess());
-for (IncEdgeIt it(*_aux_graph, second_node); it != INVALID; ++it) {
-Node node = _aux_graph->runningNode(it);
+{
-if (edges[node] == INVALID) {
+int i = 0;
-edges[node] = _aux_graph->addEdge(new_node, node);
+while (i < (int)addedges.size()) {
-(*_aux_capacity)[edges[node]] = (*_aux_capacity)[it];
+Node sn = addedges[i].source;
-} else {
+Node tn = addedges[i].target;
-(*_aux_capacity)[edges[node]] += (*_aux_capacity)[it];
+Value ec = addedges[i].capacity;
+++i;
+while (i < (int)addedges.size() &&
+sn == addedges[i].source && tn == addedges[i].target) {
+ec += addedges[i].capacity;
+++i;
+}
+typename AuxGraph::UEdge ne = _aux_graph->addEdge(sn, tn);
+(*_aux_capacity)[ne] = ec;
 }
 }
-_aux_graph->erase(second_node);
+for (typename Ufe::ClassIt c(ufe); c != INVALID; ++c) {
---_node_num;
+if (ufe.size(c) == 1) continue;
+Value cutvalue = 0;
+for (typename AuxGraph::IncEdgeIt e(*_aux_graph, c);
+e != INVALID; ++e) {
+cutvalue += (*_aux_capacity)[e];
+}
+(*_aux_cut_value)[c] = cutvalue;
+}
+for (int i = 0; i < (int)remnodes.size(); ++i) {
+_aux_graph->erase(remnodes[i]);
+}
 return _node_num == 1;
 }
 /// \brief Executes the algorithm.
 ///
 /// It sets the nodes of one of the two partitions to true in
 /// the given BoolNodeMap. The map contains a valid cut if the
 /// map have been set false previously.
 template <typename NodeMap>
 Value quickMinCut(NodeMap& nodeMap) const {
-for (typename Graph::Node it = _first_node;
+for (int i = 0; i < (int)_cut.size(); ++i) {
-it != _last_node; it = (*_next)[it]) {
+nodeMap.set(_cut[i], true);
-nodeMap.set(it, true);
+}
-}
-nodeMap.set(_last_node, true);
 return minCut();
 }
 /// \brief Returns a min cut in a NodeMap.
 ///
 return minCut();
 }
 ///@}
+private:
 };
 }
 #endif

changeset 2343	21587bc5922b
parent 2284	05ff57dc401d
child 2376	0ed45a6c74b1