Index: lemon/steiner.h
===================================================================
--- lemon/steiner.h	(revision 2382)
+++ lemon/steiner.h	(revision 2382)
@@ -0,0 +1,277 @@
+/* -*- 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_STEINER_H
+#define LEMON_STEINER_H
+
+///\ingroup approx
+///\file
+///\brief Algorithm for the 2-approximation of Steiner Tree problem.
+///
+
+#include <lemon/smart_graph.h>
+#include <lemon/graph_utils.h>
+#include <lemon/error.h>
+
+#include <lemon/ugraph_adaptor.h>
+#include <lemon/maps.h>
+
+#include <lemon/dijkstra.h>
+#include <lemon/prim.h>
+
+
+namespace lemon {
+
+  /// \ingroup approx
+  
+  /// \brief Algorithm for the 2-approximation of Steiner Tree problem
+  ///
+  /// The Steiner-tree problem is the next: Given a connected
+  /// undirected graph, a cost function on the edges and a subset of
+  /// the nodes. Construct a tree with minimum cost which covers the
+  /// given subset of the nodes. The problem is NP-hard moreover
+  /// it is APX-complete too.
+  ///
+  /// Mehlhorn's approximation algorithm is implemented in this class,
+  /// which gives a 2-approximation for the Steiner-tree problem. The
+  /// algorithm's time complexity is O(nlog(n)+e).
+  template <typename UGraph,
+            typename CostMap = typename UGraph:: template UEdgeMap<double> >
+  class SteinerTree {
+  public:
+    
+    UGRAPH_TYPEDEFS(typename UGraph)
+
+    typedef typename CostMap::Value Value;
+    
+  private:
+
+    class CompMap {
+    public:
+      typedef Node Key;
+      typedef Edge Value;
+
+    private:
+      const UGraph& _graph;
+      typename UGraph::template NodeMap<int> _comp;
+
+    public:
+      CompMap(const UGraph& graph) : _graph(graph), _comp(graph) {}
+
+      void set(const Node& node, const Edge& edge) {
+        if (edge != INVALID) {
+          _comp.set(node, _comp[_graph.source(edge)]);
+        } else {
+          _comp.set(node, -1);
+        }
+      }
+
+      int comp(const Node& node) const { return _comp[node]; }
+      void comp(const Node& node, int value) { _comp.set(node, value); }
+    };
+
+    typedef typename UGraph::template NodeMap<Edge> PredMap;
+
+    typedef ForkWriteMap<PredMap, CompMap> ForkedMap;
+
+
+    struct External {
+      int source, target;
+      UEdge uedge;
+      Value value;
+
+      External(int s, int t, const UEdge& e, const Value& v)
+        : source(s), target(t), uedge(e), value(v) {}
+    };
+
+    struct ExternalLess {
+      bool operator()(const External& left, const External& right) const {
+        return (left.source < right.source) || 
+          (left.source == right.source && left.target < right.target);
+      }
+    };
+
+
+    typedef typename UGraph::template NodeMap<bool> FilterMap;
+
+    typedef typename UGraph::template UEdgeMap<bool> TreeMap;
+
+    const UGraph& _graph;
+    const CostMap& _cost;
+
+    typename Dijkstra<UGraph, CostMap>::
+    template DefPredMap<ForkedMap>::Create _dijkstra;
+
+    PredMap* _pred;
+    CompMap* _comp;
+    ForkedMap* _forked;
+
+    int _terminal_num;
+
+    FilterMap *_filter;
+    TreeMap *_tree;
+
+  public:
+
+    /// \brief Constructor
+    
+    /// Constructor
+    ///
+    SteinerTree(const UGraph &graph, const CostMap &cost)
+      : _graph(graph), _cost(cost), _dijkstra(graph, _cost), 
+        _pred(0), _comp(0), _forked(0), _filter(0), _tree(0) {}
+
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures.
+    void init() {
+      if (!_pred) _pred = new PredMap(_graph);
+      if (!_comp) _comp = new CompMap(_graph);
+      if (!_forked) _forked = new ForkedMap(*_pred, *_comp);
+      if (!_filter) _filter = new FilterMap(_graph);
+      if (!_tree) _tree = new TreeMap(_graph);
+      _dijkstra.predMap(*_forked);
+      _dijkstra.init();
+      _terminal_num = 0;
+      for (NodeIt it(_graph); it != INVALID; ++it) {
+        _filter->set(it, false);
+      }
+    }
+
+    /// \brief Adds a new terminal node.
+    ///
+    /// Adds a new terminal node to the Steiner-tree problem.
+    void addTerminal(const Node& node) {
+      if (!_dijkstra.reached(node)) {
+        _dijkstra.addSource(node);
+        _comp->comp(node, _terminal_num);
+        ++_terminal_num;
+      }
+    }
+    
+    /// \brief Executes the algorithm.
+    ///
+    /// Executes the algorithm.
+    ///
+    /// \pre init() must be called and at least some nodes should be
+    /// added with addTerminal() before using this function.
+    ///
+    /// This method constructs an approximation of the Steiner-Tree.
+    void start() {
+      _dijkstra.start();
+      
+      std::vector<External> externals;
+      for (UEdgeIt it(_graph); it != INVALID; ++it) {
+        Node s = _graph.source(it);
+        Node t = _graph.target(it);
+        if (_comp->comp(s) == _comp->comp(t)) continue;
+
+        Value cost = _dijkstra.dist(s) + _dijkstra.dist(t) + _cost[it];
+
+        if (_comp->comp(s) < _comp->comp(t)) {
+          externals.push_back(External(_comp->comp(s), _comp->comp(t), 
+                                       it, cost));
+        } else {
+          externals.push_back(External(_comp->comp(t), _comp->comp(s), 
+                                       it, cost));
+        }
+      }
+      std::sort(externals.begin(), externals.end(), ExternalLess());
+
+      SmartUGraph aux_graph;
+      std::vector<SmartUGraph::Node> aux_nodes;
+
+      for (int i = 0; i < _terminal_num; ++i) {
+        aux_nodes.push_back(aux_graph.addNode());
+      }
+
+      SmartUGraph::UEdgeMap<Value> aux_cost(aux_graph);
+      SmartUGraph::UEdgeMap<UEdge> cross(aux_graph);
+      {
+        int i = 0;
+        while (i < (int)externals.size()) {
+          int sn = externals[i].source;
+          int tn = externals[i].target;
+          Value ev = externals[i].value;
+          UEdge ee = externals[i].uedge;
+          ++i;
+          while (i < (int)externals.size() && 
+                 sn == externals[i].source && tn == externals[i].target) {
+            if (externals[i].value < ev) {
+              ev = externals[i].value;
+              ee = externals[i].uedge;
+            }
+            ++i;
+          }
+          SmartUGraph::UEdge ne = 
+            aux_graph.addEdge(aux_nodes[sn], aux_nodes[tn]);
+          aux_cost.set(ne, ev);
+          cross.set(ne, ee);
+        }
+      }
+
+      std::vector<SmartUGraph::UEdge> aux_tree_edges;
+      BackInserterBoolMap<std::vector<SmartUGraph::UEdge> >
+        aux_tree_map(aux_tree_edges);
+      prim(aux_graph, aux_cost, aux_tree_map);
+
+      for (std::vector<SmartUGraph::UEdge>::iterator 
+             it = aux_tree_edges.begin(); it != aux_tree_edges.end(); ++it) {
+        Node node;
+        node = _graph.source(cross[*it]);
+        while (node != INVALID && !(*_filter)[node]) {
+          _filter->set(node, true);
+          node = (*_pred)[node] != INVALID ? 
+            _graph.source((*_pred)[node]) : INVALID;
+        }
+        node = _graph.target(cross[*it]);
+        while (node != INVALID && !(*_filter)[node]) {
+          _filter->set(node, true);
+          node = (*_pred)[node] != INVALID ? 
+            _graph.source((*_pred)[node]) : INVALID;
+        }
+      }
+
+      prim(nodeSubUGraphAdaptor(_graph, *_filter), _cost, *_tree);
+            
+    }
+
+    /// \brief Checks if an edge is in the Steiner-tree or not.
+    ///
+    /// Checks if an edge is in the Steiner-tree or not.
+    /// \param e is the edge that will be checked
+    /// \return \c true if e is in the Steiner-tree, \c false otherwise
+    bool tree(UEdge e){
+      return (*_tree)[e];
+    }
+
+    /// \brief Checks if the node is in the Steiner-tree or not.
+    ///
+    /// Checks if a node is in the Steiner-tree or not.
+    /// \param n is the node that will be checked
+    /// \return \c true if n is in the Steiner-tree, \c false otherwise
+    bool tree(Node n){
+      return (*_filter)[n];
+    }
+    
+
+  };
+
+} //END OF NAMESPACE LEMON
+
+#endif
