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