lemon/hypercube_graph.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 03 Apr 2009 18:59:15 +0200
changeset 600 6ac5d9ae1d3d
parent 372 7b6466ed488a
child 550 c5fd2d996909
permissions -rw-r--r--
Support real types + numerical stability fix in NS (#254)

- Real types are supported by appropriate inicialization.
- A feature of the XTI spanning tree structure is removed to ensure
numerical stability (could cause problems using integer types).
The node potentials are updated always on the lower subtree,
in order to prevent overflow problems.
The former method isn't notably faster during to our tests.
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef HYPERCUBE_GRAPH_H
    20 #define HYPERCUBE_GRAPH_H
    21 
    22 #include <vector>
    23 #include <lemon/core.h>
    24 #include <lemon/assert.h>
    25 #include <lemon/bits/graph_extender.h>
    26 
    27 ///\ingroup graphs
    28 ///\file
    29 ///\brief HypercubeGraph class.
    30 
    31 namespace lemon {
    32 
    33   class HypercubeGraphBase {
    34 
    35   public:
    36 
    37     typedef HypercubeGraphBase Graph;
    38 
    39     class Node;
    40     class Edge;
    41     class Arc;
    42 
    43   public:
    44 
    45     HypercubeGraphBase() {}
    46 
    47   protected:
    48 
    49     void construct(int dim) {
    50       LEMON_ASSERT(dim >= 1, "The number of dimensions must be at least 1.");
    51       _dim = dim;
    52       _node_num = 1 << dim;
    53       _edge_num = dim * (1 << (dim-1));
    54     }
    55 
    56   public:
    57 
    58     typedef True NodeNumTag;
    59     typedef True EdgeNumTag;
    60     typedef True ArcNumTag;
    61 
    62     int nodeNum() const { return _node_num; }
    63     int edgeNum() const { return _edge_num; }
    64     int arcNum() const { return 2 * _edge_num; }
    65 
    66     int maxNodeId() const { return _node_num - 1; }
    67     int maxEdgeId() const { return _edge_num - 1; }
    68     int maxArcId() const { return 2 * _edge_num - 1; }
    69 
    70     static Node nodeFromId(int id) { return Node(id); }
    71     static Edge edgeFromId(int id) { return Edge(id); }
    72     static Arc arcFromId(int id) { return Arc(id); }
    73 
    74     static int id(Node node) { return node._id; }
    75     static int id(Edge edge) { return edge._id; }
    76     static int id(Arc arc) { return arc._id; }
    77 
    78     Node u(Edge edge) const {
    79       int base = edge._id & ((1 << (_dim-1)) - 1);
    80       int k = edge._id >> (_dim-1);
    81       return ((base >> k) << (k+1)) | (base & ((1 << k) - 1));
    82     }
    83 
    84     Node v(Edge edge) const {
    85       int base = edge._id & ((1 << (_dim-1)) - 1);
    86       int k = edge._id >> (_dim-1);
    87       return ((base >> k) << (k+1)) | (base & ((1 << k) - 1)) | (1 << k);
    88     }
    89 
    90     Node source(Arc arc) const {
    91       return (arc._id & 1) == 1 ? u(arc) : v(arc);
    92     }
    93 
    94     Node target(Arc arc) const {
    95       return (arc._id & 1) == 1 ? v(arc) : u(arc);
    96     }
    97 
    98     typedef True FindEdgeTag;
    99     typedef True FindArcTag;
   100 
   101     Edge findEdge(Node u, Node v, Edge prev = INVALID) const {
   102       if (prev != INVALID) return INVALID;
   103       int d = u._id ^ v._id;
   104       int k = 0;
   105       if (d == 0) return INVALID;
   106       for ( ; (d & 1) == 0; d >>= 1) ++k;
   107       if (d >> 1 != 0) return INVALID;
   108       return (k << (_dim-1)) | ((u._id >> (k+1)) << k) |
   109         (u._id & ((1 << k) - 1));
   110     }
   111 
   112     Arc findArc(Node u, Node v, Arc prev = INVALID) const {
   113       Edge edge = findEdge(u, v, prev);
   114       if (edge == INVALID) return INVALID;
   115       int k = edge._id >> (_dim-1);
   116       return ((u._id >> k) & 1) == 1 ? edge._id << 1 : (edge._id << 1) | 1;
   117     }
   118 
   119     class Node {
   120       friend class HypercubeGraphBase;
   121 
   122     protected:
   123       int _id;
   124       Node(int id) : _id(id) {}
   125     public:
   126       Node() {}
   127       Node (Invalid) : _id(-1) {}
   128       bool operator==(const Node node) const {return _id == node._id;}
   129       bool operator!=(const Node node) const {return _id != node._id;}
   130       bool operator<(const Node node) const {return _id < node._id;}
   131     };
   132 
   133     class Edge {
   134       friend class HypercubeGraphBase;
   135       friend class Arc;
   136 
   137     protected:
   138       int _id;
   139 
   140       Edge(int id) : _id(id) {}
   141 
   142     public:
   143       Edge() {}
   144       Edge (Invalid) : _id(-1) {}
   145       bool operator==(const Edge edge) const {return _id == edge._id;}
   146       bool operator!=(const Edge edge) const {return _id != edge._id;}
   147       bool operator<(const Edge edge) const {return _id < edge._id;}
   148     };
   149 
   150     class Arc {
   151       friend class HypercubeGraphBase;
   152 
   153     protected:
   154       int _id;
   155 
   156       Arc(int id) : _id(id) {}
   157 
   158     public:
   159       Arc() {}
   160       Arc (Invalid) : _id(-1) {}
   161       operator Edge() const { return _id != -1 ? Edge(_id >> 1) : INVALID; }
   162       bool operator==(const Arc arc) const {return _id == arc._id;}
   163       bool operator!=(const Arc arc) const {return _id != arc._id;}
   164       bool operator<(const Arc arc) const {return _id < arc._id;}
   165     };
   166 
   167     void first(Node& node) const {
   168       node._id = _node_num - 1;
   169     }
   170 
   171     static void next(Node& node) {
   172       --node._id;
   173     }
   174 
   175     void first(Edge& edge) const {
   176       edge._id = _edge_num - 1;
   177     }
   178 
   179     static void next(Edge& edge) {
   180       --edge._id;
   181     }
   182 
   183     void first(Arc& arc) const {
   184       arc._id = 2 * _edge_num - 1;
   185     }
   186 
   187     static void next(Arc& arc) {
   188       --arc._id;
   189     }
   190 
   191     void firstInc(Edge& edge, bool& dir, const Node& node) const {
   192       edge._id = node._id >> 1;
   193       dir = (node._id & 1) == 0;
   194     }
   195 
   196     void nextInc(Edge& edge, bool& dir) const {
   197       Node n = dir ? u(edge) : v(edge);
   198       int k = (edge._id >> (_dim-1)) + 1;
   199       if (k < _dim) {
   200         edge._id = (k << (_dim-1)) |
   201           ((n._id >> (k+1)) << k) | (n._id & ((1 << k) - 1));
   202         dir = ((n._id >> k) & 1) == 0;
   203       } else {
   204         edge._id = -1;
   205         dir = true;
   206       }
   207     }
   208 
   209     void firstOut(Arc& arc, const Node& node) const {
   210       arc._id = ((node._id >> 1) << 1) | (~node._id & 1);
   211     }
   212 
   213     void nextOut(Arc& arc) const {
   214       Node n = (arc._id & 1) == 1 ? u(arc) : v(arc);
   215       int k = (arc._id >> _dim) + 1;
   216       if (k < _dim) {
   217         arc._id = (k << (_dim-1)) |
   218           ((n._id >> (k+1)) << k) | (n._id & ((1 << k) - 1));
   219         arc._id = (arc._id << 1) | (~(n._id >> k) & 1);
   220       } else {
   221         arc._id = -1;
   222       }
   223     }
   224 
   225     void firstIn(Arc& arc, const Node& node) const {
   226       arc._id = ((node._id >> 1) << 1) | (node._id & 1);
   227     }
   228 
   229     void nextIn(Arc& arc) const {
   230       Node n = (arc._id & 1) == 1 ? v(arc) : u(arc);
   231       int k = (arc._id >> _dim) + 1;
   232       if (k < _dim) {
   233         arc._id = (k << (_dim-1)) |
   234           ((n._id >> (k+1)) << k) | (n._id & ((1 << k) - 1));
   235         arc._id = (arc._id << 1) | ((n._id >> k) & 1);
   236       } else {
   237         arc._id = -1;
   238       }
   239     }
   240 
   241     static bool direction(Arc arc) {
   242       return (arc._id & 1) == 1;
   243     }
   244 
   245     static Arc direct(Edge edge, bool dir) {
   246       return Arc((edge._id << 1) | (dir ? 1 : 0));
   247     }
   248 
   249     int dimension() const {
   250       return _dim;
   251     }
   252 
   253     bool projection(Node node, int n) const {
   254       return static_cast<bool>(node._id & (1 << n));
   255     }
   256 
   257     int dimension(Edge edge) const {
   258       return edge._id >> (_dim-1);
   259     }
   260 
   261     int dimension(Arc arc) const {
   262       return arc._id >> _dim;
   263     }
   264 
   265     int index(Node node) const {
   266       return node._id;
   267     }
   268 
   269     Node operator()(int ix) const {
   270       return Node(ix);
   271     }
   272 
   273   private:
   274     int _dim;
   275     int _node_num, _edge_num;
   276   };
   277 
   278 
   279   typedef GraphExtender<HypercubeGraphBase> ExtendedHypercubeGraphBase;
   280 
   281   /// \ingroup graphs
   282   ///
   283   /// \brief Hypercube graph class
   284   ///
   285   /// This class implements a special graph type. The nodes of the graph
   286   /// are indiced with integers with at most \c dim binary digits.
   287   /// Two nodes are connected in the graph if and only if their indices
   288   /// differ only on one position in the binary form.
   289   ///
   290   /// \note The type of the indices is chosen to \c int for efficiency
   291   /// reasons. Thus the maximum dimension of this implementation is 26
   292   /// (assuming that the size of \c int is 32 bit).
   293   ///
   294   /// This graph type is fully conform to the \ref concepts::Graph
   295   /// "Graph" concept, and it also has an important extra feature
   296   /// that its maps are real \ref concepts::ReferenceMap
   297   /// "reference map"s.
   298   class HypercubeGraph : public ExtendedHypercubeGraphBase {
   299   public:
   300 
   301     typedef ExtendedHypercubeGraphBase Parent;
   302 
   303     /// \brief Constructs a hypercube graph with \c dim dimensions.
   304     ///
   305     /// Constructs a hypercube graph with \c dim dimensions.
   306     HypercubeGraph(int dim) { construct(dim); }
   307 
   308     /// \brief The number of dimensions.
   309     ///
   310     /// Gives back the number of dimensions.
   311     int dimension() const {
   312       return Parent::dimension();
   313     }
   314 
   315     /// \brief Returns \c true if the n'th bit of the node is one.
   316     ///
   317     /// Returns \c true if the n'th bit of the node is one.
   318     bool projection(Node node, int n) const {
   319       return Parent::projection(node, n);
   320     }
   321 
   322     /// \brief The dimension id of an edge.
   323     ///
   324     /// Gives back the dimension id of the given edge.
   325     /// It is in the [0..dim-1] range.
   326     int dimension(Edge edge) const {
   327       return Parent::dimension(edge);
   328     }
   329 
   330     /// \brief The dimension id of an arc.
   331     ///
   332     /// Gives back the dimension id of the given arc.
   333     /// It is in the [0..dim-1] range.
   334     int dimension(Arc arc) const {
   335       return Parent::dimension(arc);
   336     }
   337 
   338     /// \brief The index of a node.
   339     ///
   340     /// Gives back the index of the given node.
   341     /// The lower bits of the integer describes the node.
   342     int index(Node node) const {
   343       return Parent::index(node);
   344     }
   345 
   346     /// \brief Gives back a node by its index.
   347     ///
   348     /// Gives back a node by its index.
   349     Node operator()(int ix) const {
   350       return Parent::operator()(ix);
   351     }
   352 
   353     /// \brief Number of nodes.
   354     int nodeNum() const { return Parent::nodeNum(); }
   355     /// \brief Number of edges.
   356     int edgeNum() const { return Parent::edgeNum(); }
   357     /// \brief Number of arcs.
   358     int arcNum() const { return Parent::arcNum(); }
   359 
   360     /// \brief Linear combination map.
   361     ///
   362     /// This map makes possible to give back a linear combination
   363     /// for each node. It works like the \c std::accumulate function,
   364     /// so it accumulates the \c bf binary function with the \c fv first
   365     /// value. The map accumulates only on that positions (dimensions)
   366     /// where the index of the node is one. The values that have to be
   367     /// accumulated should be given by the \c begin and \c end iterators
   368     /// and the length of this range should be equal to the dimension
   369     /// number of the graph.
   370     ///
   371     ///\code
   372     /// const int DIM = 3;
   373     /// HypercubeGraph graph(DIM);
   374     /// dim2::Point<double> base[DIM];
   375     /// for (int k = 0; k < DIM; ++k) {
   376     ///   base[k].x = rnd();
   377     ///   base[k].y = rnd();
   378     /// }
   379     /// HypercubeGraph::HyperMap<dim2::Point<double> >
   380     ///   pos(graph, base, base + DIM, dim2::Point<double>(0.0, 0.0));
   381     ///\endcode
   382     ///
   383     /// \see HypercubeGraph
   384     template <typename T, typename BF = std::plus<T> >
   385     class HyperMap {
   386     public:
   387 
   388       /// \brief The key type of the map
   389       typedef Node Key;
   390       /// \brief The value type of the map
   391       typedef T Value;
   392 
   393       /// \brief Constructor for HyperMap.
   394       ///
   395       /// Construct a HyperMap for the given graph. The values that have
   396       /// to be accumulated should be given by the \c begin and \c end
   397       /// iterators and the length of this range should be equal to the
   398       /// dimension number of the graph.
   399       ///
   400       /// This map accumulates the \c bf binary function with the \c fv
   401       /// first value on that positions (dimensions) where the index of
   402       /// the node is one.
   403       template <typename It>
   404       HyperMap(const Graph& graph, It begin, It end,
   405                T fv = 0, const BF& bf = BF())
   406         : _graph(graph), _values(begin, end), _first_value(fv), _bin_func(bf)
   407       {
   408         LEMON_ASSERT(_values.size() == graph.dimension(),
   409                      "Wrong size of range");
   410       }
   411 
   412       /// \brief The partial accumulated value.
   413       ///
   414       /// Gives back the partial accumulated value.
   415       Value operator[](const Key& k) const {
   416         Value val = _first_value;
   417         int id = _graph.index(k);
   418         int n = 0;
   419         while (id != 0) {
   420           if (id & 1) {
   421             val = _bin_func(val, _values[n]);
   422           }
   423           id >>= 1;
   424           ++n;
   425         }
   426         return val;
   427       }
   428 
   429     private:
   430       const Graph& _graph;
   431       std::vector<T> _values;
   432       T _first_value;
   433       BF _bin_func;
   434     };
   435 
   436   };
   437 
   438 }
   439 
   440 #endif