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