diff options
Diffstat (limited to 'extern/quadriflow/3rd/lemon-1.3.1/lemon/nagamochi_ibaraki.h')
| -rw-r--r-- | extern/quadriflow/3rd/lemon-1.3.1/lemon/nagamochi_ibaraki.h | 702 |
1 files changed, 702 insertions, 0 deletions
diff --git a/extern/quadriflow/3rd/lemon-1.3.1/lemon/nagamochi_ibaraki.h b/extern/quadriflow/3rd/lemon-1.3.1/lemon/nagamochi_ibaraki.h new file mode 100644 index 00000000000..57a6ba64d2b --- /dev/null +++ b/extern/quadriflow/3rd/lemon-1.3.1/lemon/nagamochi_ibaraki.h @@ -0,0 +1,702 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2013 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_NAGAMOCHI_IBARAKI_H +#define LEMON_NAGAMOCHI_IBARAKI_H + + +/// \ingroup min_cut +/// \file +/// \brief Implementation of the Nagamochi-Ibaraki algorithm. + +#include <lemon/core.h> +#include <lemon/bin_heap.h> +#include <lemon/bucket_heap.h> +#include <lemon/maps.h> +#include <lemon/radix_sort.h> +#include <lemon/unionfind.h> + +#include <cassert> + +namespace lemon { + + /// \brief Default traits class for NagamochiIbaraki class. + /// + /// Default traits class for NagamochiIbaraki class. + /// \param GR The undirected graph type. + /// \param CM Type of capacity map. + template <typename GR, typename CM> + struct NagamochiIbarakiDefaultTraits { + /// The type of the capacity map. + typedef typename CM::Value Value; + + /// The undirected graph type the algorithm runs on. + typedef GR Graph; + + /// \brief The type of the map that stores the edge capacities. + /// + /// The type of the map that stores the edge capacities. + /// It must meet the \ref concepts::ReadMap "ReadMap" concept. + typedef CM CapacityMap; + + /// \brief Instantiates a CapacityMap. + /// + /// This function instantiates a \ref CapacityMap. +#ifdef DOXYGEN + static CapacityMap *createCapacityMap(const Graph& graph) +#else + static CapacityMap *createCapacityMap(const Graph&) +#endif + { + LEMON_ASSERT(false, "CapacityMap is not initialized"); + return 0; // ignore warnings + } + + /// \brief The cross reference type used by heap. + /// + /// The cross reference type used by heap. + /// Usually \c Graph::NodeMap<int>. + typedef typename Graph::template NodeMap<int> HeapCrossRef; + + /// \brief Instantiates a HeapCrossRef. + /// + /// This function instantiates a \ref HeapCrossRef. + /// \param g is the graph, to which we would like to define the + /// \ref HeapCrossRef. + static HeapCrossRef *createHeapCrossRef(const Graph& g) { + return new HeapCrossRef(g); + } + + /// \brief The heap type used by NagamochiIbaraki algorithm. + /// + /// The heap type used by NagamochiIbaraki algorithm. It has to + /// maximize the priorities. + /// + /// \sa BinHeap + /// \sa NagamochiIbaraki + typedef BinHeap<Value, HeapCrossRef, std::greater<Value> > Heap; + + /// \brief Instantiates a Heap. + /// + /// This function instantiates a \ref Heap. + /// \param r is the cross reference of the heap. + static Heap *createHeap(HeapCrossRef& r) { + return new Heap(r); + } + }; + + /// \ingroup min_cut + /// + /// \brief Calculates the minimum cut in an undirected graph. + /// + /// Calculates the minimum cut in an undirected graph with the + /// Nagamochi-Ibaraki algorithm. The algorithm separates the graph's + /// nodes into two partitions with the minimum sum of edge capacities + /// between the two partitions. The algorithm can be used to test + /// the network reliability, especially to test how many links have + /// to be destroyed in the network to split it to at least two + /// distinict subnetworks. + /// + /// The complexity of the algorithm is \f$ O(nm\log(n)) \f$ but with + /// \ref FibHeap "Fibonacci heap" it can be decreased to + /// \f$ O(nm+n^2\log(n)) \f$. When the edges have unit capacities, + /// \c BucketHeap can be used which yields \f$ O(nm) \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. + /// \note This capacity is supposed to be integer type. +#ifdef DOXYGEN + template <typename GR, typename CM, typename TR> +#else + template <typename GR, + typename CM = typename GR::template EdgeMap<int>, + typename TR = NagamochiIbarakiDefaultTraits<GR, CM> > +#endif + class NagamochiIbaraki { + public: + + typedef TR Traits; + /// The type of the underlying graph. + typedef typename Traits::Graph Graph; + + /// The type of the capacity map. + typedef typename Traits::CapacityMap CapacityMap; + /// The value type of the capacity map. + typedef typename Traits::CapacityMap::Value Value; + + /// The heap type used by the algorithm. + typedef typename Traits::Heap Heap; + /// The cross reference type used for the heap. + typedef typename Traits::HeapCrossRef HeapCrossRef; + + ///\name Named template parameters + + ///@{ + + struct SetUnitCapacityTraits : public Traits { + typedef ConstMap<typename Graph::Edge, Const<int, 1> > CapacityMap; + static CapacityMap *createCapacityMap(const Graph&) { + return new CapacityMap(); + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// the capacity map to a constMap<Edge, int, 1>() instance + /// + /// \ref named-templ-param "Named parameter" for setting + /// the capacity map to a constMap<Edge, int, 1>() instance + struct SetUnitCapacity + : public NagamochiIbaraki<Graph, CapacityMap, + SetUnitCapacityTraits> { + typedef NagamochiIbaraki<Graph, CapacityMap, + SetUnitCapacityTraits> Create; + }; + + + template <class H, class CR> + struct SetHeapTraits : public Traits { + typedef CR HeapCrossRef; + typedef H Heap; + static HeapCrossRef *createHeapCrossRef(int num) { + LEMON_ASSERT(false, "HeapCrossRef is not initialized"); + return 0; // ignore warnings + } + static Heap *createHeap(HeapCrossRef &) { + LEMON_ASSERT(false, "Heap is not initialized"); + return 0; // ignore warnings + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// heap and cross reference type + /// + /// \ref named-templ-param "Named parameter" for setting heap and + /// cross reference type. The heap has to maximize the priorities. + template <class H, class CR = RangeMap<int> > + struct SetHeap + : public NagamochiIbaraki<Graph, CapacityMap, SetHeapTraits<H, CR> > { + typedef NagamochiIbaraki< Graph, CapacityMap, SetHeapTraits<H, CR> > + Create; + }; + + template <class H, class CR> + struct SetStandardHeapTraits : public Traits { + typedef CR HeapCrossRef; + typedef H Heap; + static HeapCrossRef *createHeapCrossRef(int size) { + return new HeapCrossRef(size); + } + static Heap *createHeap(HeapCrossRef &crossref) { + return new Heap(crossref); + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// heap and cross reference type with automatic allocation + /// + /// \ref named-templ-param "Named parameter" for setting heap and + /// cross reference type with automatic allocation. They should + /// have standard constructor interfaces to be able to + /// automatically created by the algorithm (i.e. the graph should + /// be passed to the constructor of the cross reference and the + /// cross reference should be passed to the constructor of the + /// heap). However, external heap and cross reference objects + /// could also be passed to the algorithm using the \ref heap() + /// function before calling \ref run() or \ref init(). The heap + /// has to maximize the priorities. + /// \sa SetHeap + template <class H, class CR = RangeMap<int> > + struct SetStandardHeap + : public NagamochiIbaraki<Graph, CapacityMap, + SetStandardHeapTraits<H, CR> > { + typedef NagamochiIbaraki<Graph, CapacityMap, + SetStandardHeapTraits<H, CR> > Create; + }; + + ///@} + + + private: + + const Graph &_graph; + const CapacityMap *_capacity; + bool _local_capacity; // unit capacity + + struct ArcData { + typename Graph::Node target; + int prev, next; + }; + struct EdgeData { + Value capacity; + Value cut; + }; + + struct NodeData { + int first_arc; + typename Graph::Node prev, next; + int curr_arc; + typename Graph::Node last_rep; + Value sum; + }; + + typename Graph::template NodeMap<NodeData> *_nodes; + std::vector<ArcData> _arcs; + std::vector<EdgeData> _edges; + + typename Graph::Node _first_node; + int _node_num; + + Value _min_cut; + + HeapCrossRef *_heap_cross_ref; + bool _local_heap_cross_ref; + Heap *_heap; + bool _local_heap; + + typedef typename Graph::template NodeMap<typename Graph::Node> NodeList; + NodeList *_next_rep; + + typedef typename Graph::template NodeMap<bool> MinCutMap; + MinCutMap *_cut_map; + + void createStructures() { + if (!_nodes) { + _nodes = new (typename Graph::template NodeMap<NodeData>)(_graph); + } + if (!_capacity) { + _local_capacity = true; + _capacity = Traits::createCapacityMap(_graph); + } + if (!_heap_cross_ref) { + _local_heap_cross_ref = true; + _heap_cross_ref = Traits::createHeapCrossRef(_graph); + } + if (!_heap) { + _local_heap = true; + _heap = Traits::createHeap(*_heap_cross_ref); + } + if (!_next_rep) { + _next_rep = new NodeList(_graph); + } + if (!_cut_map) { + _cut_map = new MinCutMap(_graph); + } + } + + protected: + //This is here to avoid a gcc-3.3 compilation error. + //It should never be called. + NagamochiIbaraki() {} + + public: + + typedef NagamochiIbaraki Create; + + + /// \brief Constructor. + /// + /// \param graph The graph the algorithm runs on. + /// \param capacity The capacity map used by the algorithm. + NagamochiIbaraki(const Graph& graph, const CapacityMap& capacity) + : _graph(graph), _capacity(&capacity), _local_capacity(false), + _nodes(0), _arcs(), _edges(), _min_cut(), + _heap_cross_ref(0), _local_heap_cross_ref(false), + _heap(0), _local_heap(false), + _next_rep(0), _cut_map(0) {} + + /// \brief Constructor. + /// + /// This constructor can be used only when the Traits class + /// defines how can the local capacity map be instantiated. + /// If the SetUnitCapacity used the algorithm automatically + /// constructs the capacity map. + /// + ///\param graph The graph the algorithm runs on. + NagamochiIbaraki(const Graph& graph) + : _graph(graph), _capacity(0), _local_capacity(false), + _nodes(0), _arcs(), _edges(), _min_cut(), + _heap_cross_ref(0), _local_heap_cross_ref(false), + _heap(0), _local_heap(false), + _next_rep(0), _cut_map(0) {} + + /// \brief Destructor. + /// + /// Destructor. + ~NagamochiIbaraki() { + if (_local_capacity) delete _capacity; + if (_nodes) delete _nodes; + if (_local_heap) delete _heap; + if (_local_heap_cross_ref) delete _heap_cross_ref; + if (_next_rep) delete _next_rep; + if (_cut_map) delete _cut_map; + } + + /// \brief Sets the heap and the cross reference used by algorithm. + /// + /// Sets the heap and the cross reference used by algorithm. + /// If you don't use this function before calling \ref run(), + /// it will allocate one. The destuctor deallocates this + /// automatically allocated heap and cross reference, of course. + /// \return <tt> (*this) </tt> + NagamochiIbaraki &heap(Heap& hp, HeapCrossRef &cr) + { + if (_local_heap_cross_ref) { + delete _heap_cross_ref; + _local_heap_cross_ref = false; + } + _heap_cross_ref = &cr; + if (_local_heap) { + delete _heap; + _local_heap = false; + } + _heap = &hp; + return *this; + } + + /// \name Execution control + /// The simplest way to execute the algorithm is to use + /// one of the member functions called \c run(). + /// \n + /// If you need more control on the execution, + /// first you must call \ref init() and then call the start() + /// or proper times the processNextPhase() member functions. + + ///@{ + + /// \brief Initializes the internal data structures. + /// + /// Initializes the internal data structures. + void init() { + createStructures(); + + int edge_num = countEdges(_graph); + _edges.resize(edge_num); + _arcs.resize(2 * edge_num); + + typename Graph::Node prev = INVALID; + _node_num = 0; + for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) { + (*_cut_map)[n] = false; + (*_next_rep)[n] = INVALID; + (*_nodes)[n].last_rep = n; + (*_nodes)[n].first_arc = -1; + (*_nodes)[n].curr_arc = -1; + (*_nodes)[n].prev = prev; + if (prev != INVALID) { + (*_nodes)[prev].next = n; + } + (*_nodes)[n].next = INVALID; + (*_nodes)[n].sum = 0; + prev = n; + ++_node_num; + } + + _first_node = typename Graph::NodeIt(_graph); + + int index = 0; + for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) { + for (typename Graph::OutArcIt a(_graph, n); a != INVALID; ++a) { + typename Graph::Node m = _graph.target(a); + + if (!(n < m)) continue; + + (*_nodes)[n].sum += (*_capacity)[a]; + (*_nodes)[m].sum += (*_capacity)[a]; + + int c = (*_nodes)[m].curr_arc; + if (c != -1 && _arcs[c ^ 1].target == n) { + _edges[c >> 1].capacity += (*_capacity)[a]; + } else { + _edges[index].capacity = (*_capacity)[a]; + + _arcs[index << 1].prev = -1; + if ((*_nodes)[n].first_arc != -1) { + _arcs[(*_nodes)[n].first_arc].prev = (index << 1); + } + _arcs[index << 1].next = (*_nodes)[n].first_arc; + (*_nodes)[n].first_arc = (index << 1); + _arcs[index << 1].target = m; + + (*_nodes)[m].curr_arc = (index << 1); + + _arcs[(index << 1) | 1].prev = -1; + if ((*_nodes)[m].first_arc != -1) { + _arcs[(*_nodes)[m].first_arc].prev = ((index << 1) | 1); + } + _arcs[(index << 1) | 1].next = (*_nodes)[m].first_arc; + (*_nodes)[m].first_arc = ((index << 1) | 1); + _arcs[(index << 1) | 1].target = n; + + ++index; + } + } + } + + typename Graph::Node cut_node = INVALID; + _min_cut = std::numeric_limits<Value>::max(); + + for (typename Graph::Node n = _first_node; + n != INVALID; n = (*_nodes)[n].next) { + if ((*_nodes)[n].sum < _min_cut) { + cut_node = n; + _min_cut = (*_nodes)[n].sum; + } + } + (*_cut_map)[cut_node] = true; + if (_min_cut == 0) { + _first_node = INVALID; + } + } + + public: + + /// \brief Processes the next phase + /// + /// Processes the next phase in the algorithm. It must be called + /// at most one less the number of the nodes in the graph. + /// + ///\return %True when the algorithm finished. + bool processNextPhase() { + if (_first_node == INVALID) return true; + + _heap->clear(); + for (typename Graph::Node n = _first_node; + n != INVALID; n = (*_nodes)[n].next) { + (*_heap_cross_ref)[n] = Heap::PRE_HEAP; + } + + std::vector<typename Graph::Node> order; + order.reserve(_node_num); + int sep = 0; + + Value alpha = 0; + Value pmc = std::numeric_limits<Value>::max(); + + _heap->push(_first_node, static_cast<Value>(0)); + while (!_heap->empty()) { + typename Graph::Node n = _heap->top(); + Value v = _heap->prio(); + + _heap->pop(); + for (int a = (*_nodes)[n].first_arc; a != -1; a = _arcs[a].next) { + switch (_heap->state(_arcs[a].target)) { + case Heap::PRE_HEAP: + { + Value nv = _edges[a >> 1].capacity; + _heap->push(_arcs[a].target, nv); + _edges[a >> 1].cut = nv; + } break; + case Heap::IN_HEAP: + { + Value nv = _edges[a >> 1].capacity + (*_heap)[_arcs[a].target]; + _heap->decrease(_arcs[a].target, nv); + _edges[a >> 1].cut = nv; + } break; + case Heap::POST_HEAP: + break; + } + } + + alpha += (*_nodes)[n].sum; + alpha -= 2 * v; + + order.push_back(n); + if (!_heap->empty()) { + if (alpha < pmc) { + pmc = alpha; + sep = order.size(); + } + } + } + + if (static_cast<int>(order.size()) < _node_num) { + _first_node = INVALID; + for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) { + (*_cut_map)[n] = false; + } + for (int i = 0; i < static_cast<int>(order.size()); ++i) { + typename Graph::Node n = order[i]; + while (n != INVALID) { + (*_cut_map)[n] = true; + n = (*_next_rep)[n]; + } + } + _min_cut = 0; + return true; + } + + if (pmc < _min_cut) { + for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) { + (*_cut_map)[n] = false; + } + for (int i = 0; i < sep; ++i) { + typename Graph::Node n = order[i]; + while (n != INVALID) { + (*_cut_map)[n] = true; + n = (*_next_rep)[n]; + } + } + _min_cut = pmc; + } + + for (typename Graph::Node n = _first_node; + n != INVALID; n = (*_nodes)[n].next) { + bool merged = false; + for (int a = (*_nodes)[n].first_arc; a != -1; a = _arcs[a].next) { + if (!(_edges[a >> 1].cut < pmc)) { + if (!merged) { + for (int b = (*_nodes)[n].first_arc; b != -1; b = _arcs[b].next) { + (*_nodes)[_arcs[b].target].curr_arc = b; + } + merged = true; + } + typename Graph::Node m = _arcs[a].target; + int nb = 0; + for (int b = (*_nodes)[m].first_arc; b != -1; b = nb) { + nb = _arcs[b].next; + if ((b ^ a) == 1) continue; + typename Graph::Node o = _arcs[b].target; + int c = (*_nodes)[o].curr_arc; + if (c != -1 && _arcs[c ^ 1].target == n) { + _edges[c >> 1].capacity += _edges[b >> 1].capacity; + (*_nodes)[n].sum += _edges[b >> 1].capacity; + if (_edges[b >> 1].cut < _edges[c >> 1].cut) { + _edges[b >> 1].cut = _edges[c >> 1].cut; + } + if (_arcs[b ^ 1].prev != -1) { + _arcs[_arcs[b ^ 1].prev].next = _arcs[b ^ 1].next; + } else { + (*_nodes)[o].first_arc = _arcs[b ^ 1].next; + } + if (_arcs[b ^ 1].next != -1) { + _arcs[_arcs[b ^ 1].next].prev = _arcs[b ^ 1].prev; + } + } else { + if (_arcs[a].next != -1) { + _arcs[_arcs[a].next].prev = b; + } + _arcs[b].next = _arcs[a].next; + _arcs[b].prev = a; + _arcs[a].next = b; + _arcs[b ^ 1].target = n; + + (*_nodes)[n].sum += _edges[b >> 1].capacity; + (*_nodes)[o].curr_arc = b; + } + } + + if (_arcs[a].prev != -1) { + _arcs[_arcs[a].prev].next = _arcs[a].next; + } else { + (*_nodes)[n].first_arc = _arcs[a].next; + } + if (_arcs[a].next != -1) { + _arcs[_arcs[a].next].prev = _arcs[a].prev; + } + + (*_nodes)[n].sum -= _edges[a >> 1].capacity; + (*_next_rep)[(*_nodes)[n].last_rep] = m; + (*_nodes)[n].last_rep = (*_nodes)[m].last_rep; + + if ((*_nodes)[m].prev != INVALID) { + (*_nodes)[(*_nodes)[m].prev].next = (*_nodes)[m].next; + } else{ + _first_node = (*_nodes)[m].next; + } + if ((*_nodes)[m].next != INVALID) { + (*_nodes)[(*_nodes)[m].next].prev = (*_nodes)[m].prev; + } + --_node_num; + } + } + } + + if (_node_num == 1) { + _first_node = INVALID; + return true; + } + + return false; + } + + /// \brief Executes the algorithm. + /// + /// Executes the algorithm. + /// + /// \pre init() must be called + void start() { + while (!processNextPhase()) {} + } + + + /// \brief Runs %NagamochiIbaraki algorithm. + /// + /// This method runs the %Min cut algorithm + /// + /// \note mc.run(s) is just a shortcut of the following code. + ///\code + /// mc.init(); + /// mc.start(); + ///\endcode + void run() { + init(); + start(); + } + + ///@} + + /// \name Query Functions + /// + /// The result of the %NagamochiIbaraki + /// algorithm can be obtained using these functions.\n + /// Before the use of these functions, either run() or start() + /// must be called. + + ///@{ + + /// \brief Returns the min cut value. + /// + /// Returns the min cut value if the algorithm finished. + /// After the first processNextPhase() it is a value of a + /// valid cut in the graph. + Value minCutValue() const { + return _min_cut; + } + + /// \brief Returns a min cut in a NodeMap. + /// + /// It sets the nodes of one of the two partitions to true and + /// the other partition to false. + /// \param cutMap A \ref concepts::WriteMap "writable" node map with + /// \c bool (or convertible) value type. + template <typename CutMap> + Value minCutMap(CutMap& cutMap) const { + for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) { + cutMap.set(n, (*_cut_map)[n]); + } + return minCutValue(); + } + + ///@} + + }; +} + +#endif |