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