/* -*- C++ -*-

  *

  * This file is a part of LEMON, a generic C++ optimization library

  *

  * Copyright (C) 2003-2006

  * 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_MIN_COST_ARBORESCENCE_H

 #define LEMON_MIN_COST_ARBORESCENCE_H

 ///\ingroup spantree

 ///\file

 ///\brief Minimum Cost Arborescence algorithm.

 #include <vector>

 #include <lemon/list_graph.h>

 namespace lemon {

   /// \brief Default traits class of MinCostArborescence class.

   /// 

   /// Default traits class of MinCostArborescence class.

   /// \param _Graph Graph type.

   /// \param _CostMap Type of cost map.

   template <class _Graph, class _CostMap>

   struct MinCostArborescenceDefaultTraits{

     /// \brief The graph type the algorithm runs on. 

     typedef _Graph Graph;

     /// \brief The type of the map that stores the edge costs.

     ///

     /// The type of the map that stores the edge costs.

     /// It must meet the \ref concept::ReadMap "ReadMap" concept.

     typedef _CostMap CostMap;

     /// \brief The value type of the costs.

     ///

     /// The value type of the costs.

     typedef typename CostMap::Value Value;

     /// \brief The type of the map that stores which edges are 

     /// in the arborescence.

     ///

     /// The type of the map that stores which edges are in the arborescence.

     /// It must meet the \ref concept::ReadWriteMap "ReadWriteMap" concept.

     /// Initially it will be setted to false on each edge. The algorithm

     /// may set each value one time to true and maybe after it to false again.

     /// Therefore you cannot use maps like BackInserteBoolMap with this

     /// algorithm.   

     typedef typename Graph::template EdgeMap<bool> ArborescenceMap; 

     /// \brief Instantiates a ArborescenceMap.

     ///

     /// This function instantiates a \ref ArborescenceMap. 

     /// \param _graph is the graph, to which we would like to define the 

     /// ArborescenceMap.

     static ArborescenceMap *createArborescenceMap(const Graph &_graph){

       return new ArborescenceMap(_graph);

     }

   };

   /// \ingroup spantree

   ///

   /// \brief %MinCostArborescence algorithm class.

   ///

   /// This class provides an efficient implementation of 

   /// %MinCostArborescence algorithm. The arborescence is a tree 

   /// which is directed from a given source node of the graph. One or

   /// more sources should be given for the algorithm and it will calculate

   /// the minimum cost subgraph which are union of arborescences with the

   /// given sources and spans all the nodes which are reachable from the

   /// sources. The time complexity of the algorithm is O(n^2 + e).

   ///

   /// \param _Graph The graph type the algorithm runs on. The default value

   /// is \ref ListGraph. The value of _Graph is not used directly by

   /// MinCostArborescence, it is only passed to 

   /// \ref MinCostArborescenceDefaultTraits.

   /// \param _CostMap This read-only EdgeMap determines the costs of the

   /// edges. It is read once for each edge, so the map may involve in

   /// relatively time consuming process to compute the edge cost if

   /// it is necessary. The default map type is \ref

   /// concept::StaticGraph::EdgeMap "Graph::EdgeMap<int>".  The value

   /// of _CostMap is not used directly by MinCostArborescence, 

   /// it is only passed to \ref MinCostArborescenceDefaultTraits.  

   /// \param _Traits Traits class to set various data types used 

   /// by the algorithm.  The default traits class is 

   /// \ref MinCostArborescenceDefaultTraits

   /// "MinCostArborescenceDefaultTraits<_Graph,_CostMap>".  See \ref

   /// MinCostArborescenceDefaultTraits for the documentation of a 

   /// MinCostArborescence traits class.

   ///

   /// \author Balazs Dezso

 #ifndef DOXYGEN

   template <typename _Graph = ListGraph, 

             typename _CostMap = typename _Graph::template EdgeMap<int>,

             typename _Traits = 

             MinCostArborescenceDefaultTraits<_Graph, _CostMap> >

 #else 

   template <typename _Graph, typename _CostMap, typedef _Traits>

 #endif

   class MinCostArborescence {

   public:

     /// \brief \ref Exception for uninitialized parameters.

     ///

     /// This error represents problems in the initialization

     /// of the parameters of the algorithms.    

     class UninitializedParameter : public lemon::UninitializedParameter {

     public:

       virtual const char* exceptionName() const {

 	return "lemon::MinCostArborescence::UninitializedParameter";

       }

     };

     /// The traits.

     typedef _Traits Traits;

     /// The type of the underlying graph.

     typedef typename Traits::Graph Graph;

     /// The type of the map that stores the edge costs.

     typedef typename Traits::CostMap CostMap;

     ///The type of the costs of the edges.

     typedef typename Traits::Value Value;

     ///The type of the map that stores which edges are in the arborescence.

     typedef typename Traits::ArborescenceMap ArborescenceMap;

   protected:

     typedef typename Graph::Node Node;

     typedef typename Graph::Edge Edge;

     typedef typename Graph::NodeIt NodeIt;

     typedef typename Graph::EdgeIt EdgeIt;

     typedef typename Graph::InEdgeIt InEdgeIt;

     typedef typename Graph::OutEdgeIt OutEdgeIt;

     struct CostEdge {

       Edge edge;

       Value value;

       CostEdge() {}

       CostEdge(Edge _edge, Value _value) : edge(_edge), value(_value) {}

     };

     const Graph* graph;

     const CostMap* cost;

     ArborescenceMap* _arborescence_map;

     bool local_arborescence_map;

     typedef typename Graph::template NodeMap<int> LevelMap;

     LevelMap *_level;

     typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;

     CostEdgeMap *_cost_edges; 

     struct StackLevel {

       std::vector<CostEdge> edges;

       int node_level;

     };

     std::vector<StackLevel> level_stack;    

     std::vector<Node> queue;

     int node_counter;

   public:

     /// \name Named template parameters

     /// @{

     template <class T>

     struct DefArborescenceMapTraits : public Traits {

       typedef T ArborescenceMap;

       static ArborescenceMap *createArborescenceMap(const Graph &)

       {

 	throw UninitializedParameter();

       }

     };

     /// \brief \ref named-templ-param "Named parameter" for 

     /// setting ArborescenceMap type

     ///

     /// \ref named-templ-param "Named parameter" for setting 

     /// ArborescenceMap type

     template <class T>

     struct DefArborescenceMap 

       : public MinCostArborescence<Graph, CostMap,

                                    DefArborescenceMapTraits<T> > {

       typedef MinCostArborescence<Graph, CostMap, 

                                    DefArborescenceMapTraits<T> > Create;

     };

     /// @}

     /// \brief Constructor.

     ///

     /// \param _graph The graph the algorithm will run on.

     /// \param _cost The cost map used by the algorithm.

     MinCostArborescence(const Graph& _graph, const CostMap& _cost) 

       : graph(&_graph), cost(&_cost),

         _arborescence_map(0), local_arborescence_map(false), 

         _level(0), _cost_edges(0) {}

     /// \brief Destructor.

     ~MinCostArborescence() {

       destroyStructures();

     }

     /// \brief Sets the arborescence map.

     /// 

     /// Sets the arborescence map.

     /// \return \c (*this)

     MinCostArborescence& arborescenceMap(ArborescenceMap& m) {

       _arborescence_map = &m;

       return *this;

     }

     /// \name Query Functions

     /// The result of the %MinCostArborescence algorithm can be obtained 

     /// using these functions.\n

     /// Before the use of these functions,

     /// either run() or start() must be called.

     /// @{

     /// \brief Returns a reference to the arborescence map.

     ///

     /// Returns a reference to the arborescence map.

     const ArborescenceMap& arborescenceMap() const {

       return *_arborescence_map;

     }

     /// \brief Returns true if the edge is in the arborescence.

     ///

     /// Returns true if the edge is in the arborescence.

     /// \param edge The edge of the graph.

     /// \pre \ref run() must be called before using this function.

     bool arborescenceEdge(Edge edge) const {

       return (*_arborescence_map)[edge];

     }

     /// \brief Returns the cost of the arborescence.

     ///

     /// Returns the cost of the arborescence.

    Value arborescenceCost() const {

       Value sum = 0;

       for (EdgeIt it(*graph); it != INVALID; ++it) {

         if (arborescenceEdge(it)) {

           sum += (*cost)[it];

         }

       }

       return sum;

     }

     /// @}

     /// \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(), then you can add several 

     /// source nodes with \ref addSource().

     /// Finally \ref start() will perform the actual path

     /// computation.

     ///@{

     /// \brief Initializes the internal data structures.

     ///

     /// Initializes the internal data structures.

     ///

     void init() {

       initStructures();

       for (NodeIt it(*graph); it != INVALID; ++it) {

         (*_cost_edges)[it].edge = INVALID;

         (*_level)[it] = -3; 

       }

       for (EdgeIt it(*graph); it != INVALID; ++it) {

         _arborescence_map->set(it, false);

       }

     }

     /// \brief Adds a new source node.

     ///

     /// Adds a new source node to the algorithm.

     void addSource(Node source) {

       std::vector<Node> nodes;

       nodes.push_back(source);

       while (!nodes.empty()) {

         Node node = nodes.back();

         nodes.pop_back();

         for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {

           if ((*_level)[graph->target(it)] == -3) {

             (*_level)[graph->target(it)] = -2;

             nodes.push_back(graph->target(it));

             queue.push_back(graph->target(it));

           }

         }

       }

       (*_level)[source] = -1;

     }

     /// \brief Processes the next node in the priority queue.

     ///

     /// Processes the next node in the priority queue.

     ///

     /// \return The processed node.

     ///

     /// \warning The queue must not be empty!

     Node processNextNode() {

       node_counter = 0;

       Node node = queue.back();

       queue.pop_back();

       if ((*_level)[node] == -2) {

         Edge edge = prepare(node);

         while ((*_level)[graph->source(edge)] != -1) {

           if ((*_level)[graph->source(edge)] >= 0) {

             edge = contract(bottom((*_level)[graph->source(edge)]));

           } else {

             edge = prepare(graph->source(edge));

           }

         }

         finalize(graph->target(edge));

         level_stack.clear();        

       }

       return node;

     }

     /// \brief Returns the number of the nodes to be processed.

     ///

     /// Returns the number of the nodes to be processed.

     int queueSize() const {

       return queue.size();

     }

     /// \brief Returns \c false if there are nodes to be processed.

     ///

     /// Returns \c false if there are nodes to be processed.

     bool emptyQueue() const {

       return queue.empty();

     }

     /// \brief Executes the algorithm.

     ///

     /// Executes the algorithm.

     ///

     /// \pre init() must be called and at least one node should be added

     /// with addSource() before using this function.

     ///

     ///\note mca.start() is just a shortcut of the following code.

     ///\code

     ///while (!mca.emptyQueue()) {

     ///  mca.processNextNode();

     ///}

     ///\endcode

     void start() {

       while (!emptyQueue()) {

         processNextNode();

       }

     }

     /// \brief Runs %MinCostArborescence algorithm from node \c s.

     /// 

     /// This method runs the %MinCostArborescence algorithm from 

     /// a root node \c s.

     ///

     ///\note mca.run(s) is just a shortcut of the following code.

     ///\code

     ///mca.init();

     ///mca.addSource(s);

     ///mca.start();

     ///\endcode

     void run(Node node) {

       init();

       addSource(node);

       start();

     }

     ///@}

   protected:

     void initStructures() {

       if (!_arborescence_map) {

         local_arborescence_map = true;

         _arborescence_map = Traits::createArborescenceMap(*graph);

       }

       if (!_level) {

         _level = new LevelMap(*graph);

       }

       if (!_cost_edges) {

         _cost_edges = new CostEdgeMap(*graph);

       }

     }

     void destroyStructures() {

       if (_level) {

         delete _level;

       }

       if (!_cost_edges) {

         delete _cost_edges;

       }

       if (local_arborescence_map) {

         delete _arborescence_map;

       }

     }

     Edge prepare(Node node) {

       std::vector<Node> nodes;

       (*_level)[node] = node_counter;

       for (InEdgeIt it(*graph, node); it != INVALID; ++it) {

         Edge edge = it;

         Value value = (*cost)[it];

         if (graph->source(edge) == node || 

             (*_level)[graph->source(edge)] == -3) continue;

         if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {

           (*_cost_edges)[graph->source(edge)].edge = edge;

           (*_cost_edges)[graph->source(edge)].value = value;

           nodes.push_back(graph->source(edge));

         } else {

           if ((*_cost_edges)[graph->source(edge)].value > value) {

             (*_cost_edges)[graph->source(edge)].edge = edge;

             (*_cost_edges)[graph->source(edge)].value = value;

           }

         }      

       }

       CostEdge minimum = (*_cost_edges)[nodes[0]]; 

       for (int i = 1; i < (int)nodes.size(); ++i) {

         if ((*_cost_edges)[nodes[i]].value < minimum.value) {

           minimum = (*_cost_edges)[nodes[i]];

         }

       }

       StackLevel level;

       level.node_level = node_counter;

       for (int i = 0; i < (int)nodes.size(); ++i) {

         (*_cost_edges)[nodes[i]].value -= minimum.value;

         level.edges.push_back((*_cost_edges)[nodes[i]]);

         (*_cost_edges)[nodes[i]].edge = INVALID;

       }

       level_stack.push_back(level);

       ++node_counter;

       _arborescence_map->set(minimum.edge, true);

       return minimum.edge;

     }

     Edge contract(int node_bottom) {

       std::vector<Node> nodes;

       while (!level_stack.empty() && 

              level_stack.back().node_level >= node_bottom) {

         for (int i = 0; i < (int)level_stack.back().edges.size(); ++i) {

           Edge edge = level_stack.back().edges[i].edge;

           Value value = level_stack.back().edges[i].value;

           if ((*_level)[graph->source(edge)] >= node_bottom) continue;

           if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {

             (*_cost_edges)[graph->source(edge)].edge = edge;

             (*_cost_edges)[graph->source(edge)].value = value;

             nodes.push_back(graph->source(edge));

           } else {

             if ((*_cost_edges)[graph->source(edge)].value > value) {

               (*_cost_edges)[graph->source(edge)].edge = edge;

               (*_cost_edges)[graph->source(edge)].value = value;

             }

           }

         }

         level_stack.pop_back();

       }

       CostEdge minimum = (*_cost_edges)[nodes[0]]; 

       for (int i = 1; i < (int)nodes.size(); ++i) {

         if ((*_cost_edges)[nodes[i]].value < minimum.value) {

           minimum = (*_cost_edges)[nodes[i]];

         }

       }

       StackLevel level;

       level.node_level = node_bottom;

       for (int i = 0; i < (int)nodes.size(); ++i) {

         (*_cost_edges)[nodes[i]].value -= minimum.value;

         level.edges.push_back((*_cost_edges)[nodes[i]]);

         (*_cost_edges)[nodes[i]].edge = INVALID;

       }

       level_stack.push_back(level);

       _arborescence_map->set(minimum.edge, true);

       return minimum.edge;

     }

     int bottom(int level) {

       int k = level_stack.size() - 1;

       while (level_stack[k].node_level > level) {

         --k;

       }

       return level_stack[k].node_level;

     }

     void finalize(Node source) {

       std::vector<Node> nodes;

       nodes.push_back(source);

       while (!nodes.empty()) {

         Node node = nodes.back();

         nodes.pop_back();

         for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {

           if ((*_level)[graph->target(it)] >= 0 && (*_arborescence_map)[it]) {

             (*_level)[graph->target(it)] = -1;

             nodes.push_back(graph->target(it));

           } else {

             _arborescence_map->set(it, false);

           }

         }

       }

       (*_level)[source] = -1;      

     }

   };

   /// \ingroup spantree

   ///

   /// \brief Function type interface for MinCostArborescence algorithm.

   ///

   /// Function type interface for MinCostArborescence algorithm.

   /// \param graph The Graph that the algorithm runs on.

   /// \param cost The CostMap of the edges.

   /// \param source The source of the arborescence.

   /// \retval arborescence The bool EdgeMap which stores the arborescence.

   /// \return The cost of the arborescence. 

   ///

   /// \sa MinCostArborescence

   template <typename Graph, typename CostMap, typename ArborescenceMap>

   typename CostMap::Value minCostArborescence(const Graph& graph, 

                                               const CostMap& cost,

                                               typename Graph::Node source,

                                               ArborescenceMap& arborescence) {

     typename MinCostArborescence<Graph, CostMap>

       ::template DefArborescenceMap<ArborescenceMap>

       ::Create mca(graph, cost);

     mca.arborescenceMap(arborescence);

     mca.run(source);

     return mca.arborescenceCost();

   }

 }

 #endif

 // Hilbert - Huang