Index: CMakeLists.txt
===================================================================
--- CMakeLists.txt	(revision 1185)
+++ CMakeLists.txt	(revision 1138)
@@ -87,34 +87,35 @@
 IF(LEMON_ENABLE_GLPK) 
   FIND_PACKAGE(GLPK 4.33)
-  IF(GLPK_FOUND)
-    SET(LEMON_HAVE_LP TRUE)
-    SET(LEMON_HAVE_MIP TRUE)
-    SET(LEMON_HAVE_GLPK TRUE)
-  ENDIF(GLPK_FOUND)
 ENDIF(LEMON_ENABLE_GLPK)
 IF(LEMON_ENABLE_ILOG)
   FIND_PACKAGE(ILOG)
-  IF(ILOG_FOUND)
-    SET(LEMON_HAVE_LP TRUE)
-    SET(LEMON_HAVE_MIP TRUE)
-    SET(LEMON_HAVE_CPLEX TRUE)
-  ENDIF(ILOG_FOUND)
 ENDIF(LEMON_ENABLE_ILOG)
 IF(LEMON_ENABLE_COIN)
   FIND_PACKAGE(COIN)
-  IF(COIN_FOUND)
-    SET(LEMON_HAVE_LP TRUE)
-    SET(LEMON_HAVE_MIP TRUE)
-    SET(LEMON_HAVE_CLP TRUE)
-    SET(LEMON_HAVE_CBC TRUE)
-  ENDIF(COIN_FOUND)
 ENDIF(LEMON_ENABLE_COIN)
 IF(LEMON_ENABLE_SOPLEX)
   FIND_PACKAGE(SOPLEX)
-  IF(SOPLEX_FOUND)
-    SET(LEMON_HAVE_LP TRUE)
-    SET(LEMON_HAVE_SOPLEX TRUE)
-  ENDIF(SOPLEX_FOUND)
 ENDIF(LEMON_ENABLE_SOPLEX)
+
+IF(GLPK_FOUND)
+  SET(LEMON_HAVE_LP TRUE)
+  SET(LEMON_HAVE_MIP TRUE)
+  SET(LEMON_HAVE_GLPK TRUE)
+ENDIF(GLPK_FOUND)
+IF(ILOG_FOUND)
+  SET(LEMON_HAVE_LP TRUE)
+  SET(LEMON_HAVE_MIP TRUE)
+  SET(LEMON_HAVE_CPLEX TRUE)
+ENDIF(ILOG_FOUND)
+IF(COIN_FOUND)
+  SET(LEMON_HAVE_LP TRUE)
+  SET(LEMON_HAVE_MIP TRUE)
+  SET(LEMON_HAVE_CLP TRUE)
+  SET(LEMON_HAVE_CBC TRUE)
+ENDIF(COIN_FOUND)
+IF(SOPLEX_FOUND)
+  SET(LEMON_HAVE_LP TRUE)
+  SET(LEMON_HAVE_SOPLEX TRUE)
+ENDIF(SOPLEX_FOUND)
 
 IF(ILOG_FOUND)
Index: cmake/FindCOIN.cmake
===================================================================
--- cmake/FindCOIN.cmake	(revision 1180)
+++ cmake/FindCOIN.cmake	(revision 1064)
@@ -66,10 +66,4 @@
 )
 
-FIND_LIBRARY(COIN_PTHREADS_LIBRARY
-  NAMES pthreads libpthreads
-  HINTS ${COIN_ROOT_DIR}/lib/coin
-  HINTS ${COIN_ROOT_DIR}/lib
-)
-
 INCLUDE(FindPackageHandleStandardArgs)
 FIND_PACKAGE_HANDLE_STANDARD_ARGS(COIN DEFAULT_MSG
@@ -89,5 +83,5 @@
 IF(COIN_FOUND)
   SET(COIN_INCLUDE_DIRS ${COIN_INCLUDE_DIR})
-  SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARY};${COIN_COIN_UTILS_LIBRARY}")
+  SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARY};${COIN_COIN_UTILS_LIBRARY};${COIN_ZLIB_LIBRARY};${COIN_BZ2_LIBRARY}")
   IF(COIN_ZLIB_LIBRARY)
     SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_ZLIB_LIBRARY}")
@@ -96,8 +90,5 @@
     SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_BZ2_LIBRARY}")
   ENDIF(COIN_BZ2_LIBRARY)
-   IF(COIN_PTHREADS_LIBRARY)
-    SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_PTHREADS_LIBRARY}")
-  ENDIF(COIN_PTHREADS_LIBRARY)
-  SET(COIN_CBC_LIBRARIES "${COIN_CBC_LIBRARY};${COIN_CBC_SOLVER_LIBRARY};${COIN_CGL_LIBRARY};${COIN_OSI_LIBRARY};${COIN_OSI_CBC_LIBRARY};${COIN_OSI_CLP_LIBRARY};${COIN_CLP_LIBRARIES}")
+  SET(COIN_CBC_LIBRARIES "${COIN_CBC_LIBRARY};${COIN_CBC_SOLVER_LIBRARY};${COIN_CGL_LIBRARY};${COIN_OSI_LIBRARY};${COIN_OSI_CBC_LIBRARY};${COIN_OSI_CLP_LIBRARY};${COIN_ZLIB_LIBRARY};${COIN_BZ2_LIBRARY};${COIN_CLP_LIBRARIES}")
   SET(COIN_LIBRARIES ${COIN_CBC_LIBRARIES})
 ENDIF(COIN_FOUND)
Index: lemon/adaptors.h
===================================================================
--- lemon/adaptors.h	(revision 1184)
+++ lemon/adaptors.h	(revision 1092)
@@ -3447,8 +3447,8 @@
     /// This map adaptor class adapts two node maps of the original digraph
     /// to get a node map of the split digraph.
-    /// Its value type is inherited from the first node map type (\c In).
-    /// \tparam In The type of the node map for the in-nodes.
-    /// \tparam Out The type of the node map for the out-nodes.
-    template <typename In, typename Out>
+    /// Its value type is inherited from the first node map type (\c IN).
+    /// \tparam IN The type of the node map for the in-nodes.
+    /// \tparam OUT The type of the node map for the out-nodes.
+    template <typename IN, typename OUT>
     class CombinedNodeMap {
     public:
@@ -3457,14 +3457,14 @@
       typedef Node Key;
       /// The value type of the map
-      typedef typename In::Value Value;
-
-      typedef typename MapTraits<In>::ReferenceMapTag ReferenceMapTag;
-      typedef typename MapTraits<In>::ReturnValue ReturnValue;
-      typedef typename MapTraits<In>::ConstReturnValue ConstReturnValue;
-      typedef typename MapTraits<In>::ReturnValue Reference;
-      typedef typename MapTraits<In>::ConstReturnValue ConstReference;
+      typedef typename IN::Value Value;
+
+      typedef typename MapTraits<IN>::ReferenceMapTag ReferenceMapTag;
+      typedef typename MapTraits<IN>::ReturnValue ReturnValue;
+      typedef typename MapTraits<IN>::ConstReturnValue ConstReturnValue;
+      typedef typename MapTraits<IN>::ReturnValue Reference;
+      typedef typename MapTraits<IN>::ConstReturnValue ConstReference;
 
       /// Constructor
-      CombinedNodeMap(In& in_map, Out& out_map)
+      CombinedNodeMap(IN& in_map, OUT& out_map)
         : _in_map(in_map), _out_map(out_map) {}
 
@@ -3498,6 +3498,6 @@
     private:
 
-      In& _in_map;
-      Out& _out_map;
+      IN& _in_map;
+      OUT& _out_map;
 
     };
@@ -3507,26 +3507,26 @@
     ///
     /// This function just returns a combined node map.
-    template <typename In, typename Out>
-    static CombinedNodeMap<In, Out>
-    combinedNodeMap(In& in_map, Out& out_map) {
-      return CombinedNodeMap<In, Out>(in_map, out_map);
-    }
-
-    template <typename In, typename Out>
-    static CombinedNodeMap<const In, Out>
-    combinedNodeMap(const In& in_map, Out& out_map) {
-      return CombinedNodeMap<const In, Out>(in_map, out_map);
-    }
-
-    template <typename In, typename Out>
-    static CombinedNodeMap<In, const Out>
-    combinedNodeMap(In& in_map, const Out& out_map) {
-      return CombinedNodeMap<In, const Out>(in_map, out_map);
-    }
-
-    template <typename In, typename Out>
-    static CombinedNodeMap<const In, const Out>
-    combinedNodeMap(const In& in_map, const Out& out_map) {
-      return CombinedNodeMap<const In, const Out>(in_map, out_map);
+    template <typename IN, typename OUT>
+    static CombinedNodeMap<IN, OUT>
+    combinedNodeMap(IN& in_map, OUT& out_map) {
+      return CombinedNodeMap<IN, OUT>(in_map, out_map);
+    }
+
+    template <typename IN, typename OUT>
+    static CombinedNodeMap<const IN, OUT>
+    combinedNodeMap(const IN& in_map, OUT& out_map) {
+      return CombinedNodeMap<const IN, OUT>(in_map, out_map);
+    }
+
+    template <typename IN, typename OUT>
+    static CombinedNodeMap<IN, const OUT>
+    combinedNodeMap(IN& in_map, const OUT& out_map) {
+      return CombinedNodeMap<IN, const OUT>(in_map, out_map);
+    }
+
+    template <typename IN, typename OUT>
+    static CombinedNodeMap<const IN, const OUT>
+    combinedNodeMap(const IN& in_map, const OUT& out_map) {
+      return CombinedNodeMap<const IN, const OUT>(in_map, out_map);
     }
 
Index: lemon/arg_parser.cc
===================================================================
--- lemon/arg_parser.cc	(revision 1179)
+++ lemon/arg_parser.cc	(revision 877)
@@ -222,7 +222,7 @@
   {
     Opts::iterator o = _opts.find(opt);
+    Opts::iterator s = _opts.find(syn);
     LEMON_ASSERT(o!=_opts.end(), "Unknown option: '"+opt+"'");
-    LEMON_ASSERT(_opts.find(syn)==_opts.end(),
-                 "Option already used: '"+syn+"'");
+    LEMON_ASSERT(s==_opts.end(), "Option already used: '"+syn+"'");
     ParData p;
     p.help=opt;
Index: lemon/bits/vf2_internals.h
===================================================================
--- lemon/bits/vf2_internals.h	(revision 1186)
+++ lemon/bits/vf2_internals.h	(revision 1186)
@@ -0,0 +1,48 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2015-2017
+ * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
+ *
+ * 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 VF2_INTERNALS_H
+#define VF2_INTERNALS_H
+
+
+///\ingroup graph_properties
+///\file
+///\brief Mapping types for graph matching algorithms.
+
+namespace lemon {
+  ///\ingroup graph_isomorphism
+  ///The \ref Vf2 "VF2" algorithm is capable of finding different kind of
+  ///embeddings, this enum specifies its type.
+  ///
+  ///See \ref graph_isomorphism for a more detailed description.
+  enum MappingType {
+    /// Subgraph isomorphism
+    SUBGRAPH = 0,
+    /// Induced subgraph isomorphism
+    INDUCED = 1,
+    /// Graph isomorphism
+
+    /// If the two graph has the same number of nodes, than it is
+    /// equivalent to \ref INDUCED, and if they also have the same
+    /// number of edges, then it is also equivalent to \ref SUBGRAPH.
+    ///
+    /// However, using this setting is faster than the other two
+    /// options.
+    ISOMORPH = 2
+  };
+}
+#endif
Index: lemon/lgf_reader.h
===================================================================
--- lemon/lgf_reader.h	(revision 1161)
+++ lemon/lgf_reader.h	(revision 1150)
@@ -452,11 +452,11 @@
   ///
   ///\code
-  /// DigraphReader<DGR>(digraph, std::cin)
-  ///   .nodeMap("coordinates", coord_map)
-  ///   .arcMap("capacity", cap_map)
-  ///   .node("source", src)
-  ///   .node("target", trg)
-  ///   .attribute("caption", caption)
-  ///   .run();
+  /// DigraphReader<DGR>(digraph, std::cin).
+  ///   nodeMap("coordinates", coord_map).
+  ///   arcMap("capacity", cap_map).
+  ///   node("source", src).
+  ///   node("target", trg).
+  ///   attribute("caption", caption).
+  ///   run();
   ///\endcode
   ///
@@ -1247,9 +1247,9 @@
   ///ListDigraph::ArcMap<int> cap(digraph);
   ///ListDigraph::Node src, trg;
-  ///digraphReader(digraph, std::cin)
-  ///  .arcMap("capacity", cap)
-  ///  .node("source", src)
-  ///  .node("target", trg)
-  ///  .run();
+  ///digraphReader(digraph, std::cin).
+  ///  arcMap("capacity", cap).
+  ///  node("source", src).
+  ///  node("target", trg).
+  ///  run();
   ///\endcode
   ///
@@ -2124,7 +2124,7 @@
   ///ListGraph graph;
   ///ListGraph::EdgeMap<int> weight(graph);
-  ///graphReader(graph, std::cin)
-  ///  .edgeMap("weight", weight)
-  ///  .run();
+  ///graphReader(graph, std::cin).
+  ///  edgeMap("weight", weight).
+  ///  run();
   ///\endcode
   ///
@@ -3192,7 +3192,7 @@
   ///ListBpGraph graph;
   ///ListBpGraph::EdgeMap<int> weight(graph);
-  ///bpGraphReader(graph, std::cin)
-  ///  .edgeMap("weight", weight)
-  ///  .run();
+  ///bpGraphReader(graph, std::cin).
+  ///  edgeMap("weight", weight).
+  ///  run();
   ///\endcode
   ///
Index: lemon/lgf_writer.h
===================================================================
--- lemon/lgf_writer.h	(revision 1161)
+++ lemon/lgf_writer.h	(revision 1092)
@@ -409,13 +409,13 @@
   ///
   ///\code
-  /// DigraphWriter<DGR>(digraph, std::cout)
-  ///   .nodeMap("coordinates", coord_map)
-  ///   .nodeMap("size", size)
-  ///   .nodeMap("title", title)
-  ///   .arcMap("capacity", cap_map)
-  ///   .node("source", src)
-  ///   .node("target", trg)
-  ///   .attribute("caption", caption)
-  ///   .run();
+  /// DigraphWriter<DGR>(digraph, std::cout).
+  ///   nodeMap("coordinates", coord_map).
+  ///   nodeMap("size", size).
+  ///   nodeMap("title", title).
+  ///   arcMap("capacity", cap_map).
+  ///   node("source", src).
+  ///   node("target", trg).
+  ///   attribute("caption", caption).
+  ///   run();
   ///\endcode
   ///
@@ -962,9 +962,9 @@
   ///ListDigraph::Node src, trg;
   ///  // Setting the capacity map and source and target nodes
-  ///digraphWriter(digraph, std::cout)
-  ///  .arcMap("capacity", cap)
-  ///  .node("source", src)
-  ///  .node("target", trg)
-  ///  .run();
+  ///digraphWriter(digraph, std::cout).
+  ///  arcMap("capacity", cap).
+  ///  node("source", src).
+  ///  node("target", trg).
+  ///  run();
   ///\endcode
   ///
@@ -1600,7 +1600,7 @@
   ///ListGraph::EdgeMap<int> weight(graph);
   ///  // Setting the weight map
-  ///graphWriter(graph, std::cout)
-  ///  .edgeMap("weight", weight)
-  ///  .run();
+  ///graphWriter(graph, std::cout).
+  ///  edgeMap("weight", weight).
+  ///  run();
   ///\endcode
   ///
@@ -2420,7 +2420,7 @@
   ///ListBpGraph::EdgeMap<int> weight(graph);
   ///  // Setting the weight map
-  ///bpGraphWriter(graph, std::cout)
-  ///  .edgeMap("weight", weight)
-  ///  .run();
+  ///bpGraphWriter(graph, std::cout).
+  ///  edgeMap("weight", weight).
+  ///  run();
   ///\endcode
   ///
Index: lemon/planarity.h
===================================================================
--- lemon/planarity.h	(revision 1182)
+++ lemon/planarity.h	(revision 1092)
@@ -2384,5 +2384,5 @@
       if (!pe.run()) return false;
 
-      run(pe.embeddingMap());
+      run(pe);
       return true;
     }
@@ -2398,13 +2398,4 @@
     void run(const EmbeddingMap& embedding) {
       typedef SmartEdgeSet<Graph> AuxGraph;
-
-      if (countNodes(_graph) < 3) {
-        int y = 0;
-        for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
-          _point_map[n].x = 0;
-          _point_map[n].y = y++;
-        }
-        return;
-      }
 
       if (3 * countNodes(_graph) - 6 == countEdges(_graph)) {
Index: lemon/preflow.h
===================================================================
--- lemon/preflow.h	(revision 1169)
+++ lemon/preflow.h	(revision 1092)
@@ -477,5 +477,5 @@
     /// flow to the given \c flowMap. The \c flowMap should contain a
     /// flow or at least a preflow, i.e. at each node excluding the
-    /// source node the incoming flow should be greater or equal to the
+    /// source node the incoming flow should greater or equal to the
     /// outgoing flow.
     /// \return \c false if the given \c flowMap is not a preflow.
@@ -496,5 +496,5 @@
           excess -= (*_flow)[e];
         }
-        if (_tolerance.negative(excess) && n != _source) return false;
+        if (excess < 0 && n != _source) return false;
         (*_excess)[n] = excess;
       }
@@ -640,5 +640,5 @@
           (*_excess)[n] = excess;
 
-          if (_tolerance.nonZero(excess)) {
+          if (excess != 0) {
             if (new_level + 1 < _level->maxLevel()) {
               _level->liftHighestActive(new_level + 1);
@@ -721,5 +721,5 @@
           (*_excess)[n] = excess;
 
-          if (_tolerance.nonZero(excess)) {
+          if (excess != 0) {
             if (new_level + 1 < _level->maxLevel()) {
               _level->liftActiveOn(level, new_level + 1);
@@ -792,5 +792,5 @@
         if (!reached[n]) {
           _level->dirtyTopButOne(n);
-        } else if (_tolerance.positive((*_excess)[n]) && _target != n) {
+        } else if ((*_excess)[n] > 0 && _target != n) {
           _level->activate(n);
         }
@@ -853,5 +853,5 @@
         (*_excess)[n] = excess;
 
-        if (_tolerance.nonZero(excess)) {
+        if (excess != 0) {
           if (new_level + 1 < _level->maxLevel()) {
             _level->liftHighestActive(new_level + 1);
Index: lemon/random.h
===================================================================
--- lemon/random.h	(revision 1185)
+++ lemon/random.h	(revision 1136)
@@ -111,4 +111,5 @@
       static const Word loMask = (1u << 31) - 1;
       static const Word hiMask = ~loMask;
+
 
       static Word tempering(Word rnd) {
@@ -242,4 +243,5 @@
 
     private:
+
 
       void fillState() {
@@ -270,4 +272,5 @@
       }
 
+
       Word *current;
       Word state[length];
@@ -294,5 +297,5 @@
     template <typename Result, typename Word,
               int rest = std::numeric_limits<Result>::digits, int shift = 0,
-              bool last = (rest <= std::numeric_limits<Word>::digits)>
+              bool last = rest <= std::numeric_limits<Word>::digits>
     struct IntConversion {
       static const int bits = std::numeric_limits<Word>::digits;
@@ -340,5 +343,5 @@
           num = rnd() & mask;
         } while (num > max);
-        return static_cast<Result>(num);
+        return num;
       }
     };
@@ -463,544 +466,5 @@
     };
 
-    /// \ingroup misc
-    ///
-    /// \brief Mersenne Twister random number generator
-    ///
-    /// The Mersenne Twister is a twisted generalized feedback
-    /// shift-register generator of Matsumoto and Nishimura. The period
-    /// of this generator is \f$ 2^{19937} - 1\f$ and it is
-    /// equi-distributed in 623 dimensions for 32-bit numbers. The time
-    /// performance of this generator is comparable to the commonly used
-    /// generators.
-    ///
-    /// This is a template implementation of both 32-bit and
-    /// 64-bit architecture optimized versions. The generators differ
-    /// sligthly in the initialization and generation phase so they
-    /// produce two completly different sequences.
-    ///
-    /// \alert Do not use this class directly, but instead one of \ref
-    /// Random, \ref Random32 or \ref Random64.
-    ///
-    /// The generator gives back random numbers of serveral types. To
-    /// get a random number from a range of a floating point type, you
-    /// can use one form of the \c operator() or the \c real() member
-    /// function. If you want to get random number from the {0, 1, ...,
-    /// n-1} integer range, use the \c operator[] or the \c integer()
-    /// method. And to get random number from the whole range of an
-    /// integer type, you can use the argumentless \c integer() or
-    /// \c uinteger() functions. Finally, you can get random bool with
-    /// equal chance of true and false or with given probability of true
-    /// result using the \c boolean() member functions.
-    ///
-    /// Various non-uniform distributions are also supported: normal (Gauss),
-    /// exponential, gamma, Poisson, etc.; and a few two-dimensional
-    /// distributions, too.
-    ///
-    ///\code
-    /// // The commented code is identical to the other
-    /// double a = rnd();                     // [0.0, 1.0)
-    /// // double a = rnd.real();             // [0.0, 1.0)
-    /// double b = rnd(100.0);                // [0.0, 100.0)
-    /// // double b = rnd.real(100.0);        // [0.0, 100.0)
-    /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
-    /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
-    /// int d = rnd[100000];                  // 0..99999
-    /// // int d = rnd.integer(100000);       // 0..99999
-    /// int e = rnd[6] + 1;                   // 1..6
-    /// // int e = rnd.integer(1, 1 + 6);     // 1..6
-    /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
-    /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
-    /// bool g = rnd.boolean();               // P(g = true) = 0.5
-    /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
-    ///\endcode
-    ///
-    /// LEMON provides a global instance of the random number generator:
-    /// \ref lemon::rnd "rnd". In most cases, it is a good practice
-    /// to use this global generator to get random numbers.
-    ///
-    /// \sa \ref Random, \ref Random32 or \ref Random64.
-    template<class Word>
-    class Random {
-    private:
-
-      _random_bits::RandomCore<Word> core;
-      _random_bits::BoolProducer<Word> bool_producer;
-
-
-    public:
-
-      ///\name Initialization
-      ///
-      /// @{
-
-      /// \brief Default constructor
-      ///
-      /// Constructor with constant seeding.
-      Random() { core.initState(); }
-
-      /// \brief Constructor with seed
-      ///
-      /// Constructor with seed. The current number type will be converted
-      /// to the architecture word type.
-      template <typename Number>
-      Random(Number seed) {
-        _random_bits::Initializer<Number, Word>::init(core, seed);
-      }
-
-      /// \brief Constructor with array seeding
-      ///
-      /// Constructor with array seeding. The given range should contain
-      /// any number type and the numbers will be converted to the
-      /// architecture word type.
-      template <typename Iterator>
-      Random(Iterator begin, Iterator end) {
-        typedef typename std::iterator_traits<Iterator>::value_type Number;
-        _random_bits::Initializer<Number, Word>::init(core, begin, end);
-      }
-
-      /// \brief Copy constructor
-      ///
-      /// Copy constructor. The generated sequence will be identical to
-      /// the other sequence. It can be used to save the current state
-      /// of the generator and later use it to generate the same
-      /// sequence.
-      Random(const Random& other) {
-        core.copyState(other.core);
-      }
-
-      /// \brief Assign operator
-      ///
-      /// Assign operator. The generated sequence will be identical to
-      /// the other sequence. It can be used to save the current state
-      /// of the generator and later use it to generate the same
-      /// sequence.
-      Random& operator=(const Random& other) {
-        if (&other != this) {
-          core.copyState(other.core);
-        }
-        return *this;
-      }
-
-      /// \brief Seeding random sequence
-      ///
-      /// Seeding the random sequence. The current number type will be
-      /// converted to the architecture word type.
-      template <typename Number>
-      void seed(Number seed) {
-        _random_bits::Initializer<Number, Word>::init(core, seed);
-      }
-
-      /// \brief Seeding random sequence
-      ///
-      /// Seeding the random sequence. The given range should contain
-      /// any number type and the numbers will be converted to the
-      /// architecture word type.
-      template <typename Iterator>
-      void seed(Iterator begin, Iterator end) {
-        typedef typename std::iterator_traits<Iterator>::value_type Number;
-        _random_bits::Initializer<Number, Word>::init(core, begin, end);
-      }
-
-      /// \brief Seeding from file or from process id and time
-      ///
-      /// By default, this function calls the \c seedFromFile() member
-      /// function with the <tt>/dev/urandom</tt> file. If it does not success,
-      /// it uses the \c seedFromTime().
-      /// \return Currently always \c true.
-      bool seed() {
-#ifndef LEMON_WIN32
-        if (seedFromFile("/dev/urandom", 0)) return true;
-#endif
-        if (seedFromTime()) return true;
-        return false;
-      }
-
-      /// \brief Seeding from file
-      ///
-      /// Seeding the random sequence from file. The linux kernel has two
-      /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
-      /// could give good seed values for pseudo random generators (The
-      /// difference between two devices is that the <tt>random</tt> may
-      /// block the reading operation while the kernel can give good
-      /// source of randomness, while the <tt>urandom</tt> does not
-      /// block the input, but it could give back bytes with worse
-      /// entropy).
-      /// \param file The source file
-      /// \param offset The offset, from the file read.
-      /// \return \c true when the seeding successes.
-#ifndef LEMON_WIN32
-      bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
-#else
-        bool seedFromFile(const std::string& file = "", int offset = 0)
-#endif
-      {
-        std::ifstream rs(file.c_str());
-        const int size = 4;
-        Word buf[size];
-        if (offset != 0 && !rs.seekg(offset)) return false;
-        if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
-        seed(buf, buf + size);
-        return true;
-      }
-
-      /// \brief Seeding from process id and time
-      ///
-      /// Seeding from process id and time. This function uses the
-      /// current process id and the current time for initialize the
-      /// random sequence.
-      /// \return Currently always \c true.
-      bool seedFromTime() {
-#ifndef LEMON_WIN32
-        timeval tv;
-        gettimeofday(&tv, 0);
-        seed(getpid() + tv.tv_sec + tv.tv_usec);
-#else
-        seed(bits::getWinRndSeed());
-#endif
-        return true;
-      }
-
-      /// @}
-
-      ///\name Uniform Distributions
-      ///
-      /// @{
-
-      /// \brief Returns a random real number from the range [0, 1)
-      ///
-      /// It returns a random real number from the range [0, 1). The
-      /// default Number type is \c double.
-      template <typename Number>
-      Number real() {
-        return _random_bits::RealConversion<Number, Word>::convert(core);
-      }
-
-      double real() {
-        return real<double>();
-      }
-
-      /// \brief Returns a random real number from the range [0, 1)
-      ///
-      /// It returns a random double from the range [0, 1).
-      double operator()() {
-        return real<double>();
-      }
-
-      /// \brief Returns a random real number from the range [0, b)
-      ///
-      /// It returns a random real number from the range [0, b).
-      double operator()(double b) {
-        return real<double>() * b;
-      }
-
-      /// \brief Returns a random real number from the range [a, b)
-      ///
-      /// It returns a random real number from the range [a, b).
-      double operator()(double a, double b) {
-        return real<double>() * (b - a) + a;
-      }
-
-      /// \brief Returns a random integer from a range
-      ///
-      /// It returns a random integer from the range {0, 1, ..., b - 1}.
-      template <typename Number>
-      Number integer(Number b) {
-        return _random_bits::Mapping<Number, Word>::map(core, b);
-      }
-
-      /// \brief Returns a random integer from a range
-      ///
-      /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
-      template <typename Number>
-      Number integer(Number a, Number b) {
-        return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
-      }
-
-      /// \brief Returns a random integer from a range
-      ///
-      /// It returns a random integer from the range {0, 1, ..., b - 1}.
-      template <typename Number>
-      Number operator[](Number b) {
-        return _random_bits::Mapping<Number, Word>::map(core, b);
-      }
-
-      /// \brief Returns a random non-negative integer
-      ///
-      /// It returns a random non-negative integer uniformly from the
-      /// whole range of the current \c Number type. The default result
-      /// type of this function is <tt>unsigned int</tt>.
-      template <typename Number>
-      Number uinteger() {
-        return _random_bits::IntConversion<Number, Word>::convert(core);
-      }
-
-      unsigned int uinteger() {
-        return uinteger<unsigned int>();
-      }
-
-      /// \brief Returns a random integer
-      ///
-      /// It returns a random integer uniformly from the whole range of
-      /// the current \c Number type. The default result type of this
-      /// function is \c int.
-      template <typename Number>
-      Number integer() {
-        static const int nb = std::numeric_limits<Number>::digits +
-          (std::numeric_limits<Number>::is_signed ? 1 : 0);
-        return _random_bits::IntConversion<Number, Word, nb>::convert(core);
-      }
-
-      int integer() {
-        return integer<int>();
-      }
-
-      /// \brief Returns a random bool
-      ///
-      /// It returns a random bool. The generator holds a buffer for
-      /// random bits. Every time when it become empty the generator makes
-      /// a new random word and fill the buffer up.
-      bool boolean() {
-        return bool_producer.convert(core);
-      }
-
-      /// @}
-
-      ///\name Non-uniform Distributions
-      ///
-      ///@{
-
-      /// \brief Returns a random bool with given probability of true result.
-      ///
-      /// It returns a random bool with given probability of true result.
-      bool boolean(double p) {
-        return operator()() < p;
-      }
-
-      /// Standard normal (Gauss) distribution
-
-      /// Standard normal (Gauss) distribution.
-      /// \note The Cartesian form of the Box-Muller
-      /// transformation is used to generate a random normal distribution.
-      double gauss()
-      {
-        double V1,V2,S;
-        do {
-          V1=2*real<double>()-1;
-          V2=2*real<double>()-1;
-          S=V1*V1+V2*V2;
-        } while(S>=1);
-        return std::sqrt(-2*std::log(S)/S)*V1;
-      }
-      /// Normal (Gauss) distribution with given mean and standard deviation
-
-      /// Normal (Gauss) distribution with given mean and standard deviation.
-      /// \sa gauss()
-      double gauss(double mean,double std_dev)
-      {
-        return gauss()*std_dev+mean;
-      }
-
-      /// Lognormal distribution
-
-      /// Lognormal distribution. The parameters are the mean and the standard
-      /// deviation of <tt>exp(X)</tt>.
-      ///
-      double lognormal(double n_mean,double n_std_dev)
-      {
-        return std::exp(gauss(n_mean,n_std_dev));
-      }
-      /// Lognormal distribution
-
-      /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
-      /// the mean and the standard deviation of <tt>exp(X)</tt>.
-      ///
-      double lognormal(const std::pair<double,double> &params)
-      {
-        return std::exp(gauss(params.first,params.second));
-      }
-      /// Compute the lognormal parameters from mean and standard deviation
-
-      /// This function computes the lognormal parameters from mean and
-      /// standard deviation. The return value can direcly be passed to
-      /// lognormal().
-      std::pair<double,double> lognormalParamsFromMD(double mean,
-                                                     double std_dev)
-      {
-        double fr=std_dev/mean;
-        fr*=fr;
-        double lg=std::log(1+fr);
-        return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
-      }
-      /// Lognormal distribution with given mean and standard deviation
-
-      /// Lognormal distribution with given mean and standard deviation.
-      ///
-      double lognormalMD(double mean,double std_dev)
-      {
-        return lognormal(lognormalParamsFromMD(mean,std_dev));
-      }
-
-      /// Exponential distribution with given mean
-
-      /// This function generates an exponential distribution random number
-      /// with mean <tt>1/lambda</tt>.
-      ///
-      double exponential(double lambda=1.0)
-      {
-        return -std::log(1.0-real<double>())/lambda;
-      }
-
-      /// Gamma distribution with given integer shape
-
-      /// This function generates a gamma distribution random number.
-      ///
-      ///\param k shape parameter (<tt>k>0</tt> integer)
-      double gamma(int k)
-      {
-        double s = 0;
-        for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
-        return s;
-      }
-
-      /// Gamma distribution with given shape and scale parameter
-
-      /// This function generates a gamma distribution random number.
-      ///
-      ///\param k shape parameter (<tt>k>0</tt>)
-      ///\param theta scale parameter
-      ///
-      double gamma(double k,double theta=1.0)
-      {
-        double xi,nu;
-        const double delta = k-std::floor(k);
-        const double v0=E/(E-delta);
-        do {
-          double V0=1.0-real<double>();
-          double V1=1.0-real<double>();
-          double V2=1.0-real<double>();
-          if(V2<=v0)
-            {
-              xi=std::pow(V1,1.0/delta);
-              nu=V0*std::pow(xi,delta-1.0);
-            }
-          else
-            {
-              xi=1.0-std::log(V1);
-              nu=V0*std::exp(-xi);
-            }
-        } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
-        return theta*(xi+gamma(int(std::floor(k))));
-      }
-
-      /// Weibull distribution
-
-      /// This function generates a Weibull distribution random number.
-      ///
-      ///\param k shape parameter (<tt>k>0</tt>)
-      ///\param lambda scale parameter (<tt>lambda>0</tt>)
-      ///
-      double weibull(double k,double lambda)
-      {
-        return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
-      }
-
-      /// Pareto distribution
-
-      /// This function generates a Pareto distribution random number.
-      ///
-      ///\param k shape parameter (<tt>k>0</tt>)
-      ///\param x_min location parameter (<tt>x_min>0</tt>)
-      ///
-      double pareto(double k,double x_min)
-      {
-        return exponential(gamma(k,1.0/x_min))+x_min;
-      }
-
-      /// Poisson distribution
-
-      /// This function generates a Poisson distribution random number with
-      /// parameter \c lambda.
-      ///
-      /// The probability mass function of this distribusion is
-      /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
-      /// \note The algorithm is taken from the book of Donald E. Knuth titled
-      /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
-      /// return value.
-
-      int poisson(double lambda)
-      {
-        const double l = std::exp(-lambda);
-        int k=0;
-        double p = 1.0;
-        do {
-          k++;
-          p*=real<double>();
-        } while (p>=l);
-        return k-1;
-      }
-
-      ///@}
-
-      ///\name Two-Dimensional Distributions
-      ///
-      ///@{
-
-      /// Uniform distribution on the full unit circle
-
-      /// Uniform distribution on the full unit circle.
-      ///
-      dim2::Point<double> disc()
-      {
-        double V1,V2;
-        do {
-          V1=2*real<double>()-1;
-          V2=2*real<double>()-1;
-
-        } while(V1*V1+V2*V2>=1);
-        return dim2::Point<double>(V1,V2);
-      }
-      /// A kind of two-dimensional normal (Gauss) distribution
-
-      /// This function provides a turning symmetric two-dimensional distribution.
-      /// Both coordinates are of standard normal distribution, but they are not
-      /// independent.
-      ///
-      /// \note The coordinates are the two random variables provided by
-      /// the Box-Muller method.
-      dim2::Point<double> gauss2()
-      {
-        double V1,V2,S;
-        do {
-          V1=2*real<double>()-1;
-          V2=2*real<double>()-1;
-          S=V1*V1+V2*V2;
-        } while(S>=1);
-        double W=std::sqrt(-2*std::log(S)/S);
-        return dim2::Point<double>(W*V1,W*V2);
-      }
-      /// A kind of two-dimensional exponential distribution
-
-      /// This function provides a turning symmetric two-dimensional distribution.
-      /// The x-coordinate is of conditionally exponential distribution
-      /// with the condition that x is positive and y=0. If x is negative and
-      /// y=0 then, -x is of exponential distribution. The same is true for the
-      /// y-coordinate.
-      dim2::Point<double> exponential2()
-      {
-        double V1,V2,S;
-        do {
-          V1=2*real<double>()-1;
-          V2=2*real<double>()-1;
-          S=V1*V1+V2*V2;
-        } while(S>=1);
-        double W=-std::log(S)/S;
-        return dim2::Point<double>(W*V1,W*V2);
-      }
-
-      ///@}
-    };
-
-
-  };
+  }
 
   /// \ingroup misc
@@ -1008,47 +472,535 @@
   /// \brief Mersenne Twister random number generator
   ///
-  /// This class implements either the 32-bit or the 64-bit version of
-  /// the Mersenne Twister random number generator algorithm
-  /// depending on the system architecture.
-  /// 
-  /// For the API description, see its base class
-  /// \ref _random_bits::Random.
+  /// The Mersenne Twister is a twisted generalized feedback
+  /// shift-register generator of Matsumoto and Nishimura. The period
+  /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
+  /// equi-distributed in 623 dimensions for 32-bit numbers. The time
+  /// performance of this generator is comparable to the commonly used
+  /// generators.
   ///
-  /// \sa \ref _random_bits::Random
-  typedef _random_bits::Random<unsigned long> Random;
-
-  /// \ingroup misc
+  /// This implementation is specialized for both 32-bit and 64-bit
+  /// architectures. The generators differ sligthly in the
+  /// initialization and generation phase so they produce two
+  /// completly different sequences.
   ///
-  /// \brief Mersenne Twister random number generator (32-bit version)
+  /// The generator gives back random numbers of serveral types. To
+  /// get a random number from a range of a floating point type you
+  /// can use one form of the \c operator() or the \c real() member
+  /// function. If you want to get random number from the {0, 1, ...,
+  /// n-1} integer range use the \c operator[] or the \c integer()
+  /// method. And to get random number from the whole range of an
+  /// integer type you can use the argumentless \c integer() or \c
+  /// uinteger() functions. After all you can get random bool with
+  /// equal chance of true and false or given probability of true
+  /// result with the \c boolean() member functions.
   ///
-  /// This class implements the 32-bit version of the Mersenne Twister
-  /// random number generator algorithm. It is recommended to be used
-  /// when someone wants to make sure that the \e same pseudo random
-  /// sequence will be generated on every platfrom.
+  ///\code
+  /// // The commented code is identical to the other
+  /// double a = rnd();                     // [0.0, 1.0)
+  /// // double a = rnd.real();             // [0.0, 1.0)
+  /// double b = rnd(100.0);                // [0.0, 100.0)
+  /// // double b = rnd.real(100.0);        // [0.0, 100.0)
+  /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
+  /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
+  /// int d = rnd[100000];                  // 0..99999
+  /// // int d = rnd.integer(100000);       // 0..99999
+  /// int e = rnd[6] + 1;                   // 1..6
+  /// // int e = rnd.integer(1, 1 + 6);     // 1..6
+  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
+  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
+  /// bool g = rnd.boolean();               // P(g = true) = 0.5
+  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
+  ///\endcode
   ///
-  /// For the API description, see its base class
-  /// \ref _random_bits::Random.
-  ///
-  /// \sa \ref _random_bits::Random
-  typedef _random_bits::Random<unsigned int> Random32;
-
-  /// \ingroup misc
-  ///
-  /// \brief Mersenne Twister random number generator (64-bit version)
-  ///
-  /// This class implements the 64-bit version of the Mersenne Twister
-  /// random number generator algorithm. (Even though it runs
-  /// on 32-bit architectures, too.) It is recommended to be used when
-  /// someone wants to make sure that the \e same pseudo random sequence
-  /// will be generated on every platfrom.
-  ///
-  /// For the API description, see its base class
-  /// \ref _random_bits::Random.
-  ///
-  /// \sa \ref _random_bits::Random
-  typedef _random_bits::Random<unsigned long long> Random64;
+  /// LEMON provides a global instance of the random number
+  /// generator which name is \ref lemon::rnd "rnd". Usually it is a
+  /// good programming convenience to use this global generator to get
+  /// random numbers.
+  class Random {
+  private:
+
+    // Architecture word
+    typedef unsigned long Word;
+
+    _random_bits::RandomCore<Word> core;
+    _random_bits::BoolProducer<Word> bool_producer;
+
+
+  public:
+
+    ///\name Initialization
+    ///
+    /// @{
+
+    /// \brief Default constructor
+    ///
+    /// Constructor with constant seeding.
+    Random() { core.initState(); }
+
+    /// \brief Constructor with seed
+    ///
+    /// Constructor with seed. The current number type will be converted
+    /// to the architecture word type.
+    template <typename Number>
+    Random(Number seed) {
+      _random_bits::Initializer<Number, Word>::init(core, seed);
+    }
+
+    /// \brief Constructor with array seeding
+    ///
+    /// Constructor with array seeding. The given range should contain
+    /// any number type and the numbers will be converted to the
+    /// architecture word type.
+    template <typename Iterator>
+    Random(Iterator begin, Iterator end) {
+      typedef typename std::iterator_traits<Iterator>::value_type Number;
+      _random_bits::Initializer<Number, Word>::init(core, begin, end);
+    }
+
+    /// \brief Copy constructor
+    ///
+    /// Copy constructor. The generated sequence will be identical to
+    /// the other sequence. It can be used to save the current state
+    /// of the generator and later use it to generate the same
+    /// sequence.
+    Random(const Random& other) {
+      core.copyState(other.core);
+    }
+
+    /// \brief Assign operator
+    ///
+    /// Assign operator. The generated sequence will be identical to
+    /// the other sequence. It can be used to save the current state
+    /// of the generator and later use it to generate the same
+    /// sequence.
+    Random& operator=(const Random& other) {
+      if (&other != this) {
+        core.copyState(other.core);
+      }
+      return *this;
+    }
+
+    /// \brief Seeding random sequence
+    ///
+    /// Seeding the random sequence. The current number type will be
+    /// converted to the architecture word type.
+    template <typename Number>
+    void seed(Number seed) {
+      _random_bits::Initializer<Number, Word>::init(core, seed);
+    }
+
+    /// \brief Seeding random sequence
+    ///
+    /// Seeding the random sequence. The given range should contain
+    /// any number type and the numbers will be converted to the
+    /// architecture word type.
+    template <typename Iterator>
+    void seed(Iterator begin, Iterator end) {
+      typedef typename std::iterator_traits<Iterator>::value_type Number;
+      _random_bits::Initializer<Number, Word>::init(core, begin, end);
+    }
+
+    /// \brief Seeding from file or from process id and time
+    ///
+    /// By default, this function calls the \c seedFromFile() member
+    /// function with the <tt>/dev/urandom</tt> file. If it does not success,
+    /// it uses the \c seedFromTime().
+    /// \return Currently always \c true.
+    bool seed() {
+#ifndef LEMON_WIN32
+      if (seedFromFile("/dev/urandom", 0)) return true;
+#endif
+      if (seedFromTime()) return true;
+      return false;
+    }
+
+    /// \brief Seeding from file
+    ///
+    /// Seeding the random sequence from file. The linux kernel has two
+    /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
+    /// could give good seed values for pseudo random generators (The
+    /// difference between two devices is that the <tt>random</tt> may
+    /// block the reading operation while the kernel can give good
+    /// source of randomness, while the <tt>urandom</tt> does not
+    /// block the input, but it could give back bytes with worse
+    /// entropy).
+    /// \param file The source file
+    /// \param offset The offset, from the file read.
+    /// \return \c true when the seeding successes.
+#ifndef LEMON_WIN32
+    bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
+#else
+    bool seedFromFile(const std::string& file = "", int offset = 0)
+#endif
+    {
+      std::ifstream rs(file.c_str());
+      const int size = 4;
+      Word buf[size];
+      if (offset != 0 && !rs.seekg(offset)) return false;
+      if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
+      seed(buf, buf + size);
+      return true;
+    }
+
+    /// \brief Seding from process id and time
+    ///
+    /// Seding from process id and time. This function uses the
+    /// current process id and the current time for initialize the
+    /// random sequence.
+    /// \return Currently always \c true.
+    bool seedFromTime() {
+#ifndef LEMON_WIN32
+      timeval tv;
+      gettimeofday(&tv, 0);
+      seed(getpid() + tv.tv_sec + tv.tv_usec);
+#else
+      seed(bits::getWinRndSeed());
+#endif
+      return true;
+    }
+
+    /// @}
+
+    ///\name Uniform Distributions
+    ///
+    /// @{
+
+    /// \brief Returns a random real number from the range [0, 1)
+    ///
+    /// It returns a random real number from the range [0, 1). The
+    /// default Number type is \c double.
+    template <typename Number>
+    Number real() {
+      return _random_bits::RealConversion<Number, Word>::convert(core);
+    }
+
+    double real() {
+      return real<double>();
+    }
+
+    /// \brief Returns a random real number from the range [0, 1)
+    ///
+    /// It returns a random double from the range [0, 1).
+    double operator()() {
+      return real<double>();
+    }
+
+    /// \brief Returns a random real number from the range [0, b)
+    ///
+    /// It returns a random real number from the range [0, b).
+    double operator()(double b) {
+      return real<double>() * b;
+    }
+
+    /// \brief Returns a random real number from the range [a, b)
+    ///
+    /// It returns a random real number from the range [a, b).
+    double operator()(double a, double b) {
+      return real<double>() * (b - a) + a;
+    }
+
+    /// \brief Returns a random integer from a range
+    ///
+    /// It returns a random integer from the range {0, 1, ..., b - 1}.
+    template <typename Number>
+    Number integer(Number b) {
+      return _random_bits::Mapping<Number, Word>::map(core, b);
+    }
+
+    /// \brief Returns a random integer from a range
+    ///
+    /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
+    template <typename Number>
+    Number integer(Number a, Number b) {
+      return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
+    }
+
+    /// \brief Returns a random integer from a range
+    ///
+    /// It returns a random integer from the range {0, 1, ..., b - 1}.
+    template <typename Number>
+    Number operator[](Number b) {
+      return _random_bits::Mapping<Number, Word>::map(core, b);
+    }
+
+    /// \brief Returns a random non-negative integer
+    ///
+    /// It returns a random non-negative integer uniformly from the
+    /// whole range of the current \c Number type. The default result
+    /// type of this function is <tt>unsigned int</tt>.
+    template <typename Number>
+    Number uinteger() {
+      return _random_bits::IntConversion<Number, Word>::convert(core);
+    }
+
+    unsigned int uinteger() {
+      return uinteger<unsigned int>();
+    }
+
+    /// \brief Returns a random integer
+    ///
+    /// It returns a random integer uniformly from the whole range of
+    /// the current \c Number type. The default result type of this
+    /// function is \c int.
+    template <typename Number>
+    Number integer() {
+      static const int nb = std::numeric_limits<Number>::digits +
+        (std::numeric_limits<Number>::is_signed ? 1 : 0);
+      return _random_bits::IntConversion<Number, Word, nb>::convert(core);
+    }
+
+    int integer() {
+      return integer<int>();
+    }
+
+    /// \brief Returns a random bool
+    ///
+    /// It returns a random bool. The generator holds a buffer for
+    /// random bits. Every time when it become empty the generator makes
+    /// a new random word and fill the buffer up.
+    bool boolean() {
+      return bool_producer.convert(core);
+    }
+
+    /// @}
+
+    ///\name Non-uniform Distributions
+    ///
+    ///@{
+
+    /// \brief Returns a random bool with given probability of true result.
+    ///
+    /// It returns a random bool with given probability of true result.
+    bool boolean(double p) {
+      return operator()() < p;
+    }
+
+    /// Standard normal (Gauss) distribution
+
+    /// Standard normal (Gauss) distribution.
+    /// \note The Cartesian form of the Box-Muller
+    /// transformation is used to generate a random normal distribution.
+    double gauss()
+    {
+      double V1,V2,S;
+      do {
+        V1=2*real<double>()-1;
+        V2=2*real<double>()-1;
+        S=V1*V1+V2*V2;
+      } while(S>=1);
+      return std::sqrt(-2*std::log(S)/S)*V1;
+    }
+    /// Normal (Gauss) distribution with given mean and standard deviation
+
+    /// Normal (Gauss) distribution with given mean and standard deviation.
+    /// \sa gauss()
+    double gauss(double mean,double std_dev)
+    {
+      return gauss()*std_dev+mean;
+    }
+
+    /// Lognormal distribution
+
+    /// Lognormal distribution. The parameters are the mean and the standard
+    /// deviation of <tt>exp(X)</tt>.
+    ///
+    double lognormal(double n_mean,double n_std_dev)
+    {
+      return std::exp(gauss(n_mean,n_std_dev));
+    }
+    /// Lognormal distribution
+
+    /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
+    /// the mean and the standard deviation of <tt>exp(X)</tt>.
+    ///
+    double lognormal(const std::pair<double,double> &params)
+    {
+      return std::exp(gauss(params.first,params.second));
+    }
+    /// Compute the lognormal parameters from mean and standard deviation
+
+    /// This function computes the lognormal parameters from mean and
+    /// standard deviation. The return value can direcly be passed to
+    /// lognormal().
+    std::pair<double,double> lognormalParamsFromMD(double mean,
+                                                   double std_dev)
+    {
+      double fr=std_dev/mean;
+      fr*=fr;
+      double lg=std::log(1+fr);
+      return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
+    }
+    /// Lognormal distribution with given mean and standard deviation
+
+    /// Lognormal distribution with given mean and standard deviation.
+    ///
+    double lognormalMD(double mean,double std_dev)
+    {
+      return lognormal(lognormalParamsFromMD(mean,std_dev));
+    }
+
+    /// Exponential distribution with given mean
+
+    /// This function generates an exponential distribution random number
+    /// with mean <tt>1/lambda</tt>.
+    ///
+    double exponential(double lambda=1.0)
+    {
+      return -std::log(1.0-real<double>())/lambda;
+    }
+
+    /// Gamma distribution with given integer shape
+
+    /// This function generates a gamma distribution random number.
+    ///
+    ///\param k shape parameter (<tt>k>0</tt> integer)
+    double gamma(int k)
+    {
+      double s = 0;
+      for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
+      return s;
+    }
+
+    /// Gamma distribution with given shape and scale parameter
+
+    /// This function generates a gamma distribution random number.
+    ///
+    ///\param k shape parameter (<tt>k>0</tt>)
+    ///\param theta scale parameter
+    ///
+    double gamma(double k,double theta=1.0)
+    {
+      double xi,nu;
+      const double delta = k-std::floor(k);
+      const double v0=E/(E-delta);
+      do {
+        double V0=1.0-real<double>();
+        double V1=1.0-real<double>();
+        double V2=1.0-real<double>();
+        if(V2<=v0)
+          {
+            xi=std::pow(V1,1.0/delta);
+            nu=V0*std::pow(xi,delta-1.0);
+          }
+        else
+          {
+            xi=1.0-std::log(V1);
+            nu=V0*std::exp(-xi);
+          }
+      } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
+      return theta*(xi+gamma(int(std::floor(k))));
+    }
+
+    /// Weibull distribution
+
+    /// This function generates a Weibull distribution random number.
+    ///
+    ///\param k shape parameter (<tt>k>0</tt>)
+    ///\param lambda scale parameter (<tt>lambda>0</tt>)
+    ///
+    double weibull(double k,double lambda)
+    {
+      return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
+    }
+
+    /// Pareto distribution
+
+    /// This function generates a Pareto distribution random number.
+    ///
+    ///\param k shape parameter (<tt>k>0</tt>)
+    ///\param x_min location parameter (<tt>x_min>0</tt>)
+    ///
+    double pareto(double k,double x_min)
+    {
+      return exponential(gamma(k,1.0/x_min))+x_min;
+    }
+
+    /// Poisson distribution
+
+    /// This function generates a Poisson distribution random number with
+    /// parameter \c lambda.
+    ///
+    /// The probability mass function of this distribusion is
+    /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
+    /// \note The algorithm is taken from the book of Donald E. Knuth titled
+    /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
+    /// return value.
+
+    int poisson(double lambda)
+    {
+      const double l = std::exp(-lambda);
+      int k=0;
+      double p = 1.0;
+      do {
+        k++;
+        p*=real<double>();
+      } while (p>=l);
+      return k-1;
+    }
+
+    ///@}
+
+    ///\name Two Dimensional Distributions
+    ///
+    ///@{
+
+    /// Uniform distribution on the full unit circle
+
+    /// Uniform distribution on the full unit circle.
+    ///
+    dim2::Point<double> disc()
+    {
+      double V1,V2;
+      do {
+        V1=2*real<double>()-1;
+        V2=2*real<double>()-1;
+
+      } while(V1*V1+V2*V2>=1);
+      return dim2::Point<double>(V1,V2);
+    }
+    /// A kind of two dimensional normal (Gauss) distribution
+
+    /// This function provides a turning symmetric two-dimensional distribution.
+    /// Both coordinates are of standard normal distribution, but they are not
+    /// independent.
+    ///
+    /// \note The coordinates are the two random variables provided by
+    /// the Box-Muller method.
+    dim2::Point<double> gauss2()
+    {
+      double V1,V2,S;
+      do {
+        V1=2*real<double>()-1;
+        V2=2*real<double>()-1;
+        S=V1*V1+V2*V2;
+      } while(S>=1);
+      double W=std::sqrt(-2*std::log(S)/S);
+      return dim2::Point<double>(W*V1,W*V2);
+    }
+    /// A kind of two dimensional exponential distribution
+
+    /// This function provides a turning symmetric two-dimensional distribution.
+    /// The x-coordinate is of conditionally exponential distribution
+    /// with the condition that x is positive and y=0. If x is negative and
+    /// y=0 then, -x is of exponential distribution. The same is true for the
+    /// y-coordinate.
+    dim2::Point<double> exponential2()
+    {
+      double V1,V2,S;
+      do {
+        V1=2*real<double>()-1;
+        V2=2*real<double>()-1;
+        S=V1*V1+V2*V2;
+      } while(S>=1);
+      double W=-std::log(S)/S;
+      return dim2::Point<double>(W*V1,W*V2);
+    }
+
+    ///@}
+  };
+
 
   extern Random rnd;
-  
+
 }
 
Index: lemon/vf2.h
===================================================================
--- lemon/vf2.h	(revision 1156)
+++ lemon/vf2.h	(revision 1186)
@@ -3,5 +3,5 @@
  * This file is a part of LEMON, a generic C++ optimization library.
  *
- * Copyright (C) 2015
+ * Copyright (C) 2015-2017
  * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
  *
@@ -27,20 +27,16 @@
 #include <lemon/dfs.h>
 #include <lemon/bfs.h>
+#include <lemon/bits/vf2_internals.h>
 
 #include <vector>
-#include <set>
-
-namespace lemon
-{
-  namespace bits
-  {
-    namespace vf2
-    {
-      class AlwaysEq
-      {
+
+namespace lemon {
+  namespace bits {
+    namespace vf2 {
+
+      class AlwaysEq {
       public:
         template<class T1, class T2>
-        bool operator()(T1, T2) const
-        {
+        bool operator()(T1&, T2&) const {
           return true;
         }
@@ -48,19 +44,18 @@
 
       template<class M1, class M2>
-      class MapEq
-      {
+      class MapEq{
         const M1 &_m1;
         const M2 &_m2;
       public:
-        MapEq(const M1 &m1, const M2 &m2) : _m1(m1), _m2(m2) {}
-        bool operator()(typename M1::Key k1, typename M2::Key k2) const
-        {
+        MapEq(const M1 &m1, const M2 &m2) : _m1(m1), _m2(m2) { }
+        bool operator()(typename M1::Key k1, typename M2::Key k2) const {
           return _m1[k1] == _m2[k2];
         }
       };
 
+
+
       template <class G>
-      class DfsLeaveOrder : public DfsVisitor<G>
-      {
+      class DfsLeaveOrder : public DfsVisitor<G> {
         const G &_g;
         std::vector<typename G::Node> &_order;
@@ -68,8 +63,6 @@
       public:
         DfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
-          : i(countNodes(g)), _g(g), _order(order)
-        {}
-        void leave(const typename G::Node &node)
-        {
+          : i(countNodes(g)), _g(g), _order(order) { }
+        void leave(const typename G::Node &node) {
           _order[--i]=node;
         }
@@ -77,6 +70,5 @@
 
       template <class G>
-      class BfsLeaveOrder : public BfsVisitor<G>
-      {
+      class BfsLeaveOrder : public BfsVisitor<G> {
         int i;
         const G &_g;
@@ -84,8 +76,7 @@
       public:
         BfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
-          : i(0), _g(g), _order(order)
-        {}
-        void process(const typename G::Node &node)
-        {
+          : i(0), _g(g), _order(order){
+        }
+        void process(const typename G::Node &node) {
           _order[i++]=node;
         }
@@ -94,26 +85,4 @@
   }
 
-  ///Graph mapping types.
-
-  ///\ingroup graph_isomorphism
-  ///The \ref Vf2 "VF2" algorithm is capable of finding different kind of
-  ///embeddings, this enum specifies its type.
-  ///
-  ///See \ref graph_isomorphism for a more detailed description.
-  enum Vf2MappingType {
-    /// Subgraph isomorphism
-    SUBGRAPH = 0,
-    /// Induced subgraph isomorphism
-    INDUCED = 1,
-    /// Graph isomorphism
-
-    /// If the two graph has the same number of nodes, than it is
-    /// equivalent to \ref INDUCED, and if they also have the same
-    /// number of edges, then it is also equivalent to \ref SUBGRAPH.
-    ///
-    /// However, using this setting is faster than the other two
-    /// options.
-    ISOMORPH = 2
-  };
 
   ///%VF2 algorithm class.
@@ -145,10 +114,9 @@
            class NEQ = bits::vf2::AlwaysEq >
 #endif
-  class Vf2
-  {
+  class Vf2 {
     //Current depth in the DFS tree.
     int _depth;
     //Functor with bool operator()(G1::Node,G2::Node), which returns 1
-    //if and only if the 2 nodes are equivalent.
+    //ifff the two nodes are equivalent.
     NEQ _nEq;
 
@@ -157,5 +125,5 @@
     //a node of g2.
     M &_mapping;
-    //order[i] is the node of g1, for which we find a pair if depth=i
+    //order[i] is the node of g1, for which we search a pair if depth=i
     std::vector<typename G1::Node> order;
     //currEdgeIts[i] is an edge iterator, witch is last used in the ith
@@ -164,109 +132,102 @@
     //The small graph.
     const G1 &_g1;
-    //The big graph.
+    //The large graph.
     const G2 &_g2;
-    //lookup tables for cut the searchtree
+    //lookup tables for cutting the searchtree
     typename G1::template NodeMap<int> rNew1t,rInOut1t;
 
-    Vf2MappingType _mapping_type;
+    MappingType _mapping_type;
+
+    bool _deallocMappingAfterUse;
 
     //cut the search tree
-    template<Vf2MappingType MT>
-    bool cut(const typename G1::Node n1,const typename G2::Node n2) const
-    {
+    template<MappingType MT>
+    bool cut(const typename G1::Node n1,const typename G2::Node n2) const {
       int rNew2=0,rInOut2=0;
-      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
-        {
-          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
-          if(_conn[currNode]>0)
-            ++rInOut2;
-          else if(MT!=SUBGRAPH&&_conn[currNode]==0)
-            ++rNew2;
-        }
-      switch(MT)
-        {
-        case INDUCED:
-          return rInOut1t[n1]<=rInOut2&&rNew1t[n1]<=rNew2;
-        case SUBGRAPH:
-          return rInOut1t[n1]<=rInOut2;
-        case ISOMORPH:
-          return rInOut1t[n1]==rInOut2&&rNew1t[n1]==rNew2;
-        default:
-          return false;
-        }
-    }
-
-    template<Vf2MappingType MT>
-    bool feas(const typename G1::Node n1,const typename G2::Node n2)
-    {
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
+        if(_conn[currNode]>0)
+          ++rInOut2;
+        else if(MT!=SUBGRAPH&&_conn[currNode]==0)
+          ++rNew2;
+      }
+      switch(MT) {
+      case INDUCED:
+        return rInOut1t[n1]<=rInOut2&&rNew1t[n1]<=rNew2;
+      case SUBGRAPH:
+        return rInOut1t[n1]<=rInOut2;
+      case ISOMORPH:
+        return rInOut1t[n1]==rInOut2&&rNew1t[n1]==rNew2;
+      default:
+        return false;
+      }
+    }
+
+    template<MappingType MT>
+    bool feas(const typename G1::Node n1,const typename G2::Node n2) {
       if(!_nEq(n1,n2))
         return 0;
 
-      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
-        {
-          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
-          if(_mapping[currNode]!=INVALID)
-            --_conn[_mapping[currNode]];
-        }
+      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+        const typename G1::Node& currNode=_g1.oppositeNode(n1,e1);
+        if(_mapping[currNode]!=INVALID)
+          --_conn[_mapping[currNode]];
+      }
       bool isIso=1;
-      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
-        {
-          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
-          if(_conn[currNode]<-1)
-            ++_conn[currNode];
-          else if(MT!=SUBGRAPH&&_conn[currNode]==-1)
-            {
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& connCurrNode = _conn[_g2.oppositeNode(n2,e2)];
+        if(connCurrNode<-1)
+          ++connCurrNode;
+        else if(MT!=SUBGRAPH&&connCurrNode==-1) {
+          isIso=0;
+          break;
+        }
+      }
+
+      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+        const typename G2::Node& currNodePair=_mapping[_g1.oppositeNode(n1,e1)];
+        int& connCurrNodePair=_conn[currNodePair];
+        if(currNodePair!=INVALID&&connCurrNodePair!=-1) {
+          switch(MT) {
+          case INDUCED:
+          case ISOMORPH:
+            isIso=0;
+            break;
+          case SUBGRAPH:
+            if(connCurrNodePair<-1)
               isIso=0;
-              break;
-            }
-        }
-
-      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
-        {
-          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
-          if(_mapping[currNode]!=INVALID&&_conn[_mapping[currNode]]!=-1)
-            {
-              switch(MT)
-                {
-                case INDUCED:
-                case ISOMORPH:
-                  isIso=0;
-                  break;
-                case SUBGRAPH:
-                  if(_conn[_mapping[currNode]]<-1)
-                    isIso=0;
-                  break;
-                }
-              _conn[_mapping[currNode]]=-1;
-            }
-        }
+            break;
+          }
+          connCurrNodePair=-1;
+        }
+      }
       return isIso&&cut<MT>(n1,n2);
     }
 
-    void addPair(const typename G1::Node n1,const typename G2::Node n2)
-    {
+    void addPair(const typename G1::Node n1,const typename G2::Node n2) {
       _conn[n2]=-1;
       _mapping.set(n1,n2);
-      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
-        if(_conn[_g2.oppositeNode(n2,e2)]!=-1)
-          ++_conn[_g2.oppositeNode(n2,e2)];
-    }
-
-    void subPair(const typename G1::Node n1,const typename G2::Node n2)
-    {
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
+        if(currConn!=-1)
+          ++currConn;
+      }
+    }
+
+    void subPair(const typename G1::Node n1,const typename G2::Node n2) {
       _conn[n2]=0;
       _mapping.set(n1,INVALID);
-      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
-        {
-          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
-          if(_conn[currNode]>0)
-            --_conn[currNode];
-          else if(_conn[currNode]==-1)
-            ++_conn[n2];
-        }
-    }
-
-    void setOrder()//we will find pairs for the nodes of g1 in this order
-    {
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
+        if(currConn>0)
+          --currConn;
+        else if(currConn==-1)
+          ++_conn[n2];
+      }
+    }
+
+    void setOrder() {
+      // we will find pairs for the nodes of g1 in this order
+
       // bits::vf2::DfsLeaveOrder<G1> v(_g1,order);
       //   DfsVisit<G1,bits::vf2::DfsLeaveOrder<G1> >dfs(_g1, v);
@@ -279,115 +240,101 @@
     }
 
-    template<Vf2MappingType MT>
-    bool extMatch()
-    {
-      while(_depth>=0)
-        {
-          //there are not nodes in g1, which has not pair in g2.
-          if(_depth==static_cast<int>(order.size()))
-            {
+    template<MappingType MT>
+    bool extMatch() {
+      while(_depth>=0) {
+        //there are not nodes in g1, which has not pair in g2.
+        if(_depth==static_cast<int>(order.size())) {
+          --_depth;
+          return true;
+        }
+        typename G1::Node& nodeOfDepth = order[_depth];
+        const typename G2::Node& pairOfNodeOfDepth = _mapping[nodeOfDepth];
+        typename G2::IncEdgeIt &edgeItOfDepth = currEdgeIts[_depth];
+        //the node of g2, which neighbours are the candidates for
+        //the pair of nodeOfDepth
+        typename G2::Node currPNode;
+        if(edgeItOfDepth==INVALID) {
+          typename G1::IncEdgeIt fstMatchedE(_g1,nodeOfDepth);
+          //if pairOfNodeOfDepth!=INVALID, we dont use
+          //fstMatchedE
+          if(pairOfNodeOfDepth==INVALID)
+            for(; fstMatchedE!=INVALID &&
+                  _mapping[_g1.oppositeNode(nodeOfDepth,
+                                            fstMatchedE)]==INVALID;
+                ++fstMatchedE) ; //find fstMatchedE
+          if(fstMatchedE==INVALID||pairOfNodeOfDepth!=INVALID) {
+            //We found no covered neighbours, this means
+            //the graph is not connected(or _depth==0).  Each
+            //uncovered(and there are some other properties due
+            //to the spec. problem types) node of g2 is
+            //candidate.  We can read the iterator of the last
+            //tried node from the match if it is not the first
+            //try(match[nodeOfDepth]!=INVALID)
+            typename G2::NodeIt n2(_g2);
+            //if it's not the first try
+            if(pairOfNodeOfDepth!=INVALID) {
+              n2=++typename G2::NodeIt(_g2,pairOfNodeOfDepth);
+              subPair(nodeOfDepth,pairOfNodeOfDepth);
+            }
+            for(; n2!=INVALID; ++n2)
+              if(MT!=SUBGRAPH) {
+                if(_conn[n2]==0&&feas<MT>(nodeOfDepth,n2))
+                  break;
+              }
+              else if(_conn[n2]>=0&&feas<MT>(nodeOfDepth,n2))
+                break;
+            // n2 is the next candidate
+            if(n2!=INVALID){
+              addPair(nodeOfDepth,n2);
+              ++_depth;
+            }
+            else // there are no more candidates
               --_depth;
-              return true;
-            }
-          //the node of g2, which neighbours are the candidates for
-          //the pair of order[_depth]
-          typename G2::Node currPNode;
-          if(currEdgeIts[_depth]==INVALID)
-            {
-              typename G1::IncEdgeIt fstMatchedE(_g1,order[_depth]);
-              //if _mapping[order[_depth]]!=INVALID, we dont use
-              //fstMatchedE
-              if(_mapping[order[_depth]]==INVALID)
-                for(; fstMatchedE!=INVALID &&
-                      _mapping[_g1.oppositeNode(order[_depth],
-                                              fstMatchedE)]==INVALID;
-                    ++fstMatchedE) ; //find fstMatchedE
-              if(fstMatchedE==INVALID||_mapping[order[_depth]]!=INVALID)
-                {
-                  //We did not find an covered neighbour, this means
-                  //the graph is not connected(or _depth==0).  Every
-                  //uncovered(and there are some other properties due
-                  //to the spec. problem types) node of g2 is
-                  //candidate.  We can read the iterator of the last
-                  //tryed node from the match if it is not the first
-                  //try(match[order[_depth]]!=INVALID)
-                  typename G2::NodeIt n2(_g2);
-                  //if its not the first try
-                  if(_mapping[order[_depth]]!=INVALID)
-                    {
-                      n2=++typename G2::NodeIt(_g2,_mapping[order[_depth]]);
-                      subPair(order[_depth],_mapping[order[_depth]]);
-                    }
-                  for(; n2!=INVALID; ++n2)
-                    if(MT!=SUBGRAPH&&_conn[n2]==0)
-                      {
-                        if(feas<MT>(order[_depth],n2))
-                          break;
-                      }
-                    else if(MT==SUBGRAPH&&_conn[n2]>=0)
-                      if(feas<MT>(order[_depth],n2))
-                        break;
-                  // n2 is the next candidate
-                  if(n2!=INVALID)
-                    {
-                      addPair(order[_depth],n2);
-                      ++_depth;
-                    }
-                  else // there is no more candidate
-                    --_depth;
-                  continue;
-                }
-              else
-                {
-                  currPNode=_mapping[_g1.oppositeNode(order[_depth],
-                                                      fstMatchedE)];
-                  currEdgeIts[_depth]=typename G2::IncEdgeIt(_g2,currPNode);
-                }
-            }
-          else
-            {
-              currPNode=_g2.oppositeNode(_mapping[order[_depth]],
-                                         currEdgeIts[_depth]);
-              subPair(order[_depth],_mapping[order[_depth]]);
-              ++currEdgeIts[_depth];
-            }
-          for(; currEdgeIts[_depth]!=INVALID; ++currEdgeIts[_depth])
-            {
-              const typename G2::Node currNode =
-                _g2.oppositeNode(currPNode, currEdgeIts[_depth]);
-              //if currNode is uncovered
-              if(_conn[currNode]>0&&feas<MT>(order[_depth],currNode))
-                {
-                  addPair(order[_depth],currNode);
-                  break;
-                }
-            }
-          currEdgeIts[_depth]==INVALID?--_depth:++_depth;
-        }
+            continue;
+          }
+          else {
+            currPNode=_mapping[_g1.oppositeNode(nodeOfDepth,
+                                                fstMatchedE)];
+            edgeItOfDepth=typename G2::IncEdgeIt(_g2,currPNode);
+          }
+        }
+        else {
+          currPNode=_g2.oppositeNode(pairOfNodeOfDepth,
+                                     edgeItOfDepth);
+          subPair(nodeOfDepth,pairOfNodeOfDepth);
+          ++edgeItOfDepth;
+        }
+        for(; edgeItOfDepth!=INVALID; ++edgeItOfDepth) {
+          const typename G2::Node currNode =
+            _g2.oppositeNode(currPNode, edgeItOfDepth);
+          if(_conn[currNode]>0&&feas<MT>(nodeOfDepth,currNode)) {
+            addPair(nodeOfDepth,currNode);
+            break;
+          }
+        }
+        edgeItOfDepth==INVALID?--_depth:++_depth;
+      }
       return false;
     }
 
     //calc. the lookup table for cut the searchtree
-    void setRNew1tRInOut1t()
-    {
+    void setRNew1tRInOut1t() {
       typename G1::template NodeMap<int> tmp(_g1,0);
-      for(unsigned int i=0; i<order.size(); ++i)
-        {
-          tmp[order[i]]=-1;
-          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
-            {
-              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
-              if(tmp[currNode]>0)
-                ++rInOut1t[order[i]];
-              else if(tmp[currNode]==0)
-                ++rNew1t[order[i]];
-            }
-          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
-            {
-              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
-              if(tmp[currNode]!=-1)
-                ++tmp[currNode];
-            }
-        }
+      for(unsigned int i=0; i<order.size(); ++i) {
+        const typename G1::Node& orderI = order[i];
+        tmp[orderI]=-1;
+        for(typename G1::IncEdgeIt e1(_g1,orderI); e1!=INVALID; ++e1) {
+          const typename G1::Node currNode=_g1.oppositeNode(orderI,e1);
+          if(tmp[currNode]>0)
+            ++rInOut1t[orderI];
+          else if(tmp[currNode]==0)
+            ++rNew1t[orderI];
+        }
+        for(typename G1::IncEdgeIt e1(_g1,orderI); e1!=INVALID; ++e1) {
+          const typename G1::Node currNode=_g1.oppositeNode(orderI,e1);
+          if(tmp[currNode]!=-1)
+            ++tmp[currNode];
+        }
+      }
     }
   public:
@@ -400,11 +347,10 @@
     ///\param m \ref concepts::ReadWriteMap "read-write" NodeMap
     ///storing the found mapping.
-    ///\param neq A bool-valued binary functor determinining whether a node is
+    ///\param neq A bool-valued binary functor determining whether a node is
     ///mappable to another. By default it is an always true operator.
-    Vf2(const G1 &g1, const G2 &g2,M &m, const NEQ &neq = NEQ() ) :
+    Vf2(const G1 &g1, const G2 &g2, M &m, const NEQ &neq = NEQ() ) :
       _nEq(neq), _conn(g2,0), _mapping(m), order(countNodes(g1)),
       currEdgeIts(countNodes(g1),INVALID), _g1(g1), _g2(g2), rNew1t(g1,0),
-      rInOut1t(g1,0), _mapping_type(SUBGRAPH)
-    {
+      rInOut1t(g1,0), _mapping_type(SUBGRAPH), _deallocMappingAfterUse(0) {
       _depth=0;
       setOrder();
@@ -414,9 +360,21 @@
     }
 
+    ///Destructor
+
+    ///Destructor.
+    ///
+
+    ~Vf2(){
+      if(_deallocMappingAfterUse)
+        delete &_mapping;
+    }
+
     ///Returns the current mapping type
 
     ///Returns the current mapping type
     ///
-    Vf2MappingType mappingType() const { return _mapping_type; }
+    MappingType mappingType() const {
+      return _mapping_type;
+    }
     ///Sets mapping type
 
@@ -425,11 +383,13 @@
     ///The mapping type is set to \ref SUBGRAPH by default.
     ///
-    ///\sa See \ref Vf2MappingType for the possible values.
-    void mappingType(Vf2MappingType m_type) { _mapping_type = m_type; }
+    ///\sa See \ref MappingType for the possible values.
+    void mappingType(MappingType m_type) {
+      _mapping_type = m_type;
+    }
 
     ///Finds a mapping
 
-    ///It finds a mapping between from g1 into g2 according to the mapping
-    ///type set by \ref mappingType(Vf2MappingType) "mappingType()".
+    ///It finds a mapping from g1 into g2 according to the mapping
+    ///type set by \ref mappingType(MappingType) "mappingType()".
     ///
     ///By subsequent calls, it returns all possible mappings one-by-one.
@@ -437,23 +397,20 @@
     ///\retval true if a mapping is found.
     ///\retval false if there is no (more) mapping.
-    bool find()
-    {
-      switch(_mapping_type)
-        {
-        case SUBGRAPH:
-          return extMatch<SUBGRAPH>();
-        case INDUCED:
-          return extMatch<INDUCED>();
-        case ISOMORPH:
-          return extMatch<ISOMORPH>();
-        default:
-          return false;
-        }
+    bool find() {
+      switch(_mapping_type) {
+      case SUBGRAPH:
+        return extMatch<SUBGRAPH>();
+      case INDUCED:
+        return extMatch<INDUCED>();
+      case ISOMORPH:
+        return extMatch<ISOMORPH>();
+      default:
+        return false;
+      }
     }
   };
 
   template<class G1, class G2>
-  class Vf2WizardBase
-  {
+  class Vf2WizardBase {
   protected:
     typedef G1 Graph1;
@@ -463,13 +420,15 @@
     const G2 &_g2;
 
-    Vf2MappingType _mapping_type;
+    MappingType _mapping_type;
 
     typedef typename G1::template NodeMap<typename G2::Node> Mapping;
     bool _local_mapping;
     void *_mapping;
-    void createMapping()
-    {
+    void createMapping() {
       _mapping = new Mapping(_g1);
     }
+
+    void *myVf2; //used in Vf2Wizard::find
+
 
     typedef bits::vf2::AlwaysEq NodeEq;
@@ -477,6 +436,7 @@
 
     Vf2WizardBase(const G1 &g1,const G2 &g2)
-      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) {}
+      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) { }
   };
+
 
   /// Auxiliary class for the function-type interface of %VF2 algorithm.
@@ -490,6 +450,5 @@
   /// algorithm.
   template<class TR>
-  class Vf2Wizard : public TR
-  {
+  class Vf2Wizard : public TR {
     typedef TR Base;
     typedef typename TR::Graph1 Graph1;
@@ -507,12 +466,16 @@
   public:
     ///Constructor
-    Vf2Wizard(const Graph1 &g1,const Graph2 &g2) : Base(g1,g2) {}
+    Vf2Wizard(const Graph1 &g1,const Graph2 &g2) : Base(g1,g2) {
+    }
 
     ///Copy constructor
-    Vf2Wizard(const Base &b) : Base(b) {}
+    Vf2Wizard(const Base &b) : Base(b) { }
+
+    ///Copy constructor
+    Vf2Wizard(const Vf2Wizard &b) : Base(b) {}
 
 
     template<class T>
-    struct SetMappingBase : public Base {
+    struct SetMappingBase : public Base{
       typedef T Mapping;
       SetMappingBase(const Base &b) : Base(b) {}
@@ -525,6 +488,5 @@
     ///the map that stores the found embedding.
     template<class T>
-    Vf2Wizard< SetMappingBase<T> > mapping(const T &t)
-    {
+    Vf2Wizard< SetMappingBase<T> > mapping(const T &t) {
       Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
       Base::_local_mapping = false;
@@ -537,5 +499,6 @@
       NodeEq _node_eq;
       SetNodeEqBase(const Base &b, const NE &node_eq)
-        : Base(b), _node_eq(node_eq) {}
+        : Base(b), _node_eq(node_eq){
+      }
     };
 
@@ -550,6 +513,5 @@
     ///always true operator.
     template<class T>
-    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq)
-    {
+    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq) {
       return Vf2Wizard<SetNodeEqBase<T> >(SetNodeEqBase<T>(*this,node_eq));
     }
@@ -561,7 +523,7 @@
     ///the node labels defining equivalence relation between them.
     ///
-    ///\param m1 It is arbitrary \ref concepts::ReadMap "readable node map"
+    ///\param m1 An arbitrary \ref concepts::ReadMap "readable node map"
     ///of g1.
-    ///\param m2 It is arbitrary \ref concepts::ReadMap "readable node map"
+    ///\param m2 An arbitrary \ref concepts::ReadMap "readable node map"
     ///of g2.
     ///
@@ -569,6 +531,5 @@
     template<class M1, class M2>
     Vf2Wizard< SetNodeEqBase<bits::vf2::MapEq<M1,M2> > >
-    nodeLabels(const M1 &m1,const M2 &m2)
-    {
+    nodeLabels(const M1 &m1,const M2 &m2){
       return nodeEq(bits::vf2::MapEq<M1,M2>(m1,m2));
     }
@@ -582,7 +543,6 @@
     ///The mapping type is set to \ref SUBGRAPH by default.
     ///
-    ///\sa See \ref Vf2MappingType for the possible values.
-    Vf2Wizard<Base> &mappingType(Vf2MappingType m_type)
-    {
+    ///\sa See \ref MappingType for the possible values.
+    Vf2Wizard<Base> &mappingType(MappingType m_type) {
       _mapping_type = m_type;
       return *this;
@@ -594,6 +554,5 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///the mapping type to \ref INDUCED.
-    Vf2Wizard<Base> &induced()
-    {
+    Vf2Wizard<Base> &induced() {
       _mapping_type = INDUCED;
       return *this;
@@ -605,10 +564,10 @@
     ///\ref named-templ-param "Named parameter" for setting
     ///the mapping type to \ref ISOMORPH.
-    Vf2Wizard<Base> &iso()
-    {
+    Vf2Wizard<Base> &iso() {
       _mapping_type = ISOMORPH;
       return *this;
     }
 
+
     ///Runs VF2 algorithm.
 
@@ -616,7 +575,6 @@
     ///
     ///\retval true if a mapping is found.
-    ///\retval false if there is no (more) mapping.
-    bool run()
-    {
+    ///\retval false if there is no mapping.
+    bool run(){
       if(Base::_local_mapping)
         Base::createMapping();
@@ -631,4 +589,44 @@
       if(Base::_local_mapping)
         delete reinterpret_cast<Mapping*>(_mapping);
+
+      return ret;
+    }
+
+    ///Get a pointer to the generated Vf2 object.
+
+    ///Gives a pointer to the generated Vf2 object.
+    ///
+    ///\return Pointer to the generated Vf2 object.
+    ///\warning Don't forget to delete the referred Vf2 object after use.
+    Vf2<Graph1, Graph2, Mapping, NodeEq >* getPtrToVf2Object() {
+      if(Base::_local_mapping)
+        Base::createMapping();
+      Vf2<Graph1, Graph2, Mapping, NodeEq >* ptr =
+        new Vf2<Graph1, Graph2, Mapping, NodeEq>
+        (_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);
+      ptr->mappingType(_mapping_type);
+      if(Base::_local_mapping)
+        ptr->_deallocMappingAfterUse = true;
+      return ptr;
+    }
+
+    ///Counts the number of mappings.
+
+    ///This method counts the number of mappings.
+    ///
+    /// \return The number of mappings.
+    int count() {
+      if(Base::_local_mapping)
+        Base::createMapping();
+
+      Vf2<Graph1, Graph2, Mapping, NodeEq>
+        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);
+      if(Base::_local_mapping)
+        alg._deallocMappingAfterUse = true;
+      alg.mappingType(_mapping_type);
+
+      int ret = 0;
+      while(alg.find())
+        ++ret;
 
       return ret;
@@ -645,18 +643,30 @@
   ///The following examples show how to use these parameters.
   ///\code
-  ///  // Find an embedding of graph g into graph h
+  ///  // Find an embedding of graph g1 into graph g2
   ///  ListGraph::NodeMap<ListGraph::Node> m(g);
-  ///  vf2(g,h).mapping(m).run();
-  ///
-  ///  // Check whether graphs g and h are isomorphic
-  ///  bool is_iso = vf2(g,h).iso().run();
+  ///  vf2(g1,g2).mapping(m).run();
+  ///
+  ///  // Check whether graphs g1 and g2 are isomorphic
+  ///  bool is_iso = vf2(g1,g2).iso().run();
+  ///
+  ///  // Count the number of isomorphisms
+  ///  int num_isos = vf2(g1,g2).iso().count();
+  ///
+  ///  // Iterate through all the induced subgraph mappings of graph g1 into g2
+  ///  auto* myVf2 = vf2(g1,g2).mapping(m).nodeLabels(c1,c2)
+  ///  .induced().getPtrToVf2Object();
+  ///  while(myVf2->find()){
+  ///    //process the current mapping m
+  ///  }
+  ///  delete myVf22;
   ///\endcode
-  ///\warning Don't forget to put the \ref Vf2Wizard::run() "run()"
+  ///\warning Don't forget to put the \ref Vf2Wizard::run() "run()",
+  ///\ref Vf2Wizard::count() "count()" or
+  ///the \ref Vf2Wizard::getPtrToVf2Object() "getPtrToVf2Object()"
   ///to the end of the expression.
   ///\sa Vf2Wizard
   ///\sa Vf2
   template<class G1, class G2>
-  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2)
-  {
+  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2) {
     return Vf2Wizard<Vf2WizardBase<G1,G2> >(g1,g2);
   }
Index: lemon/vf2pp.h
===================================================================
--- lemon/vf2pp.h	(revision 1186)
+++ lemon/vf2pp.h	(revision 1186)
@@ -0,0 +1,902 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2015-2017
+ * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
+ *
+ * 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_VF2PP_H
+#define LEMON_VF2PP_H
+
+///\ingroup graph_properties
+///\file
+///\brief VF2 Plus Plus algorithm.
+
+#include <lemon/core.h>
+#include <lemon/concepts/graph.h>
+#include <lemon/dfs.h>
+#include <lemon/bfs.h>
+#include <lemon/bits/vf2_internals.h>
+
+
+#include <vector>
+#include <algorithm>
+#include <utility>
+
+
+namespace lemon {
+  namespace bits {
+    namespace vf2pp {
+
+      template <class G>
+      class DfsLeaveOrder : public DfsVisitor<G> {
+        int i;
+        const G &_g;
+        std::vector<typename G::Node> &_order;
+      public:
+        DfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
+          : i(countNodes(g)), _g(g), _order(order) {
+        }
+        void leave(const typename G::Node &node) {
+          _order[--i]=node;
+        }
+      };
+
+      template <class G>
+      class BfsLeaveOrder : public BfsVisitor<G> {
+        int i;
+        const G &_g;
+        std::vector<typename G::Node> &_order;
+      public:
+        BfsLeaveOrder(const G &g, std::vector<typename G::Node> &order) { }
+        void process(const typename G::Node &node) {
+          _order[i++]=node;
+        }
+      };
+    }
+  }
+
+
+  ///%VF2 Plus Plus algorithm class.
+
+  ///\ingroup graph_isomorphism This class provides an efficient
+  ///implementation of the %VF2 Plus Plus algorithm
+  ///for variants of the (Sub)graph Isomorphism problem.
+  ///
+  ///There is also a \ref vf2pp() "function-type interface" called
+  ///\ref vf2pp() for the %VF2 Plus Plus algorithm, which is probably
+  ///more convenient in most use-cases.
+  ///
+  ///\tparam G1 The type of the graph to be embedded.
+  ///The default type is \ref ListDigraph.
+  ///\tparam G2 The type of the graph g1 will be embedded into.
+  ///The default type is \ref ListDigraph.
+  ///\tparam M The type of the NodeMap storing the mapping.
+  ///By default, it is G1::NodeMap<G2::Node>
+  ///\tparam M1 The type of the NodeMap storing the integer node labels of G1.
+  ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
+  ///different labels. By default, it is G1::NodeMap<int>.
+  ///\tparam M2 The type of the NodeMap storing the integer node labels of G2.
+  ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
+  ///different labels. By default, it is G2::NodeMap<int>.
+  ///
+  ///\sa vf2pp()
+#ifdef DOXYGEN
+  template<class G1, class G2, class M, class M1, class M2 >
+#else
+  template<class G1=ListDigraph,
+           class G2=ListDigraph,
+           class M = typename G1::template NodeMap<G2::Node>,//the mapping
+           //labels of G1,the labels are the numbers {0,1,2,..,K-1},
+           //where K is the number of different labels
+           class M1 = typename G1::template NodeMap<int>,
+           //labels of G2, ...
+           class M2 = typename G2::template NodeMap<int> >
+#endif
+  class Vf2pp {
+    //Current depth in the search tree.
+    int _depth;
+
+    //_conn[v2] = number of covered neighbours of v2
+    typename G2::template NodeMap<int> _conn;
+
+    //The current mapping. _mapping[v1]=v2 iff v1 has been mapped to v2,
+    //where v1 is a node of G1 and v2 is a node of G2
+    M &_mapping;
+
+    //order[i] is a node of g1, for which a pair is searched if depth=i
+    std::vector<typename G1::Node> order;
+
+    //currEdgeIts[i] is the last used edge iterator in the ith
+    //depth to find a pair for node order[i]
+    std::vector<typename G2::IncEdgeIt> currEdgeIts;
+
+    //The small graph.
+    const G1 &_g1;
+
+    //The large graph.
+    const G2 &_g2;
+
+    //rNewLabels1[v] is a pair of form
+    //(label; num. of uncov. nodes with such label and no covered neighbours)
+    typename G1::template NodeMap<std::vector<std::pair<int,int> > >
+    rNewLabels1;
+
+    //rInOutLabels1[v] is the number of covered neighbours of v for each label
+    //in form (label,number of such labels)
+    typename G1::template NodeMap<std::vector<std::pair<int,int> > >
+    rInOutLabels1;
+
+    //_intLabels1[v]==i means that vertex v has the i label in
+    //_g1 (i is in {0,1,2,..,K-1}, where K is the number of diff. labels)
+    M1 &_intLabels1;
+
+    //_intLabels2[v]==i means that vertex v has the i label in
+    //_g2 (i is in {0,1,2,..,K-1}, where K is the number of diff. labels)
+    M2 &_intLabels2;
+
+    //largest label
+    const int maxLabel;
+
+    //lookup tables for manipulating with label class cardinalities
+    //after use they have to be reset to 0..0
+    std::vector<int> labelTmp1,labelTmp2;
+
+    MappingType _mapping_type;
+
+    //indicates whether the mapping or the labels must be deleted in the ctor
+    bool _deallocMappingAfterUse,_deallocLabelsAfterUse;
+
+
+    //improved cutting function
+    template<MappingType MT>
+    bool cutByLabels(const typename G1::Node n1,const typename G2::Node n2) {
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
+        if(_conn[currNode]>0)
+          --labelTmp1[_intLabels2[currNode]];
+        else if(MT!=SUBGRAPH&&_conn[currNode]==0)
+          --labelTmp2[_intLabels2[currNode]];
+      }
+
+      bool ret=1;
+      if(ret) {
+        for(unsigned int i = 0; i < rInOutLabels1[n1].size(); ++i)
+          labelTmp1[rInOutLabels1[n1][i].first]+=rInOutLabels1[n1][i].second;
+
+        if(MT!=SUBGRAPH)
+          for(unsigned int i = 0; i < rNewLabels1[n1].size(); ++i)
+            labelTmp2[rNewLabels1[n1][i].first]+=rNewLabels1[n1][i].second;
+
+        switch(MT) {
+        case INDUCED:
+          for(unsigned int i = 0; i < rInOutLabels1[n1].size(); ++i)
+            if(labelTmp1[rInOutLabels1[n1][i].first]>0) {
+              ret=0;
+              break;
+            }
+          if(ret)
+            for(unsigned int i = 0; i < rNewLabels1[n1].size(); ++i)
+              if(labelTmp2[rNewLabels1[n1][i].first]>0) {
+                ret=0;
+                break;
+              }
+          break;
+        case SUBGRAPH:
+          for(unsigned int i = 0; i < rInOutLabels1[n1].size(); ++i)
+            if(labelTmp1[rInOutLabels1[n1][i].first]>0) {
+              ret=0;
+              break;
+            }
+          break;
+        case ISOMORPH:
+          for(unsigned int i = 0; i < rInOutLabels1[n1].size(); ++i)
+            if(labelTmp1[rInOutLabels1[n1][i].first]!=0) {
+              ret=0;
+              break;
+            }
+          if(ret)
+            for(unsigned int i = 0; i < rNewLabels1[n1].size(); ++i)
+              if(labelTmp2[rNewLabels1[n1][i].first]!=0) {
+                ret=0;
+                break;
+              }
+          break;
+        default:
+          return false;
+        }
+        for(unsigned int i = 0; i < rInOutLabels1[n1].size(); ++i)
+          labelTmp1[rInOutLabels1[n1][i].first]=0;
+
+        if(MT!=SUBGRAPH)
+          for(unsigned int i = 0; i < rNewLabels1[n1].size(); ++i)
+            labelTmp2[rNewLabels1[n1][i].first]=0;
+      }
+
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
+        labelTmp1[_intLabels2[currNode]]=0;
+        if(MT!=SUBGRAPH&&_conn[currNode]==0)
+          labelTmp2[_intLabels2[currNode]]=0;
+      }
+
+      return ret;
+    }
+
+
+    //try to exclude the matching of n1 and n2
+    template<MappingType MT>
+    bool feas(const typename G1::Node n1,const typename G2::Node n2) {
+      if(_intLabels1[n1]!=_intLabels2[n2])
+        return 0;
+
+      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+        const typename G1::Node& currNode=_g1.oppositeNode(n1,e1);
+        if(_mapping[currNode]!=INVALID)
+          --_conn[_mapping[currNode]];
+      }
+
+      bool isIso=1;
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& connCurrNode = _conn[_g2.oppositeNode(n2,e2)];
+        if(connCurrNode<-1)
+          ++connCurrNode;
+        else if(MT!=SUBGRAPH&&connCurrNode==-1) {
+          isIso=0;
+          break;
+        }
+      }
+
+      if(isIso)
+        for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+          const typename G2::Node& currNodePair =
+            _mapping[_g1.oppositeNode(n1,e1)];
+          int& connCurrNodePair=_conn[currNodePair];
+          if(currNodePair!=INVALID&&connCurrNodePair!=-1) {
+            switch(MT){
+            case INDUCED:
+            case ISOMORPH:
+              isIso=0;
+              break;
+            case SUBGRAPH:
+              if(connCurrNodePair<-1)
+                isIso=0;
+              break;
+            }
+            connCurrNodePair=-1;
+          }
+        }
+      else
+        for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+          const typename G2::Node currNode=_mapping[_g1.oppositeNode(n1,e1)];
+          if(currNode!=INVALID/*&&_conn[currNode]!=-1*/)
+            _conn[currNode]=-1;
+        }
+
+      return isIso&&cutByLabels<MT>(n1,n2);
+    }
+
+
+    //matches n1 and n2
+    void addPair(const typename G1::Node n1,const typename G2::Node n2) {
+      _conn[n2]=-1;
+      _mapping.set(n1,n2);
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
+        if(currConn!=-1)
+          ++currConn;
+      }
+    }
+
+
+    //dematches n1 and n2
+    void subPair(const typename G1::Node n1,const typename G2::Node n2) {
+      _conn[n2]=0;
+      _mapping.set(n1,INVALID);
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2){
+        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
+        if(currConn>0)
+          --currConn;
+        else if(currConn==-1)
+          ++_conn[n2];
+      }
+    }
+
+
+    void processBFSLevel(typename G1::Node source,unsigned int& orderIndex,
+                         typename G1::template NodeMap<int>& dm1,
+                         typename G1::template NodeMap<bool>& added) {
+      order[orderIndex]=source;
+      added[source]=1;
+
+      unsigned int endPosOfLevel=orderIndex,
+        startPosOfLevel=orderIndex,
+        lastAdded=orderIndex;
+
+      typename G1::template NodeMap<int> currConn(_g1,0);
+
+      while(orderIndex<=lastAdded){
+        typename G1::Node currNode = order[orderIndex];
+        for(typename G1::IncEdgeIt e(_g1,currNode); e!=INVALID; ++e) {
+          typename G1::Node n = _g1.oppositeNode(currNode,e);
+          if(!added[n]) {
+            order[++lastAdded]=n;
+            added[n]=1;
+          }
+        }
+        if(orderIndex>endPosOfLevel){
+          for(unsigned int j = startPosOfLevel; j <= endPosOfLevel; ++j) {
+            int minInd=j;
+            for(unsigned int i = j+1; i <= endPosOfLevel; ++i)
+              if(currConn[order[i]]>currConn[order[minInd]]||
+                 (currConn[order[i]]==currConn[order[minInd]]&&
+                  (dm1[order[i]]>dm1[order[minInd]]||
+                   (dm1[order[i]]==dm1[order[minInd]]&&
+                    labelTmp1[_intLabels1[order[minInd]]]>
+                    labelTmp1[_intLabels1[order[i]]]))))
+                minInd=i;
+
+            --labelTmp1[_intLabels1[order[minInd]]];
+            for(typename G1::IncEdgeIt e(_g1,order[minInd]); e!=INVALID; ++e)
+              ++currConn[_g1.oppositeNode(order[minInd],e)];
+            std::swap(order[j],order[minInd]);
+          }
+          startPosOfLevel=endPosOfLevel+1;
+          endPosOfLevel=lastAdded;
+        }
+        ++orderIndex;
+      }
+    }
+
+
+    //we will find pairs for the nodes of g1 in this order
+    void setOrder(){
+      for(typename G2::NodeIt n2(_g2); n2!=INVALID; ++n2)
+        ++labelTmp1[_intLabels2[n2]];
+
+      //       OutDegMap<G1> dm1(_g1);
+      typename G1::template NodeMap<int> dm1(_g1,0);
+      for(typename G1::EdgeIt e(_g1); e!=INVALID; ++e) {
+        ++dm1[_g1.u(e)];
+        ++dm1[_g1.v(e)];
+      }
+
+      typename G1::template NodeMap<bool> added(_g1,0);
+      unsigned int orderIndex=0;
+
+      for(typename G1::NodeIt n(_g1); n!=INVALID;) {
+        if(!added[n]){
+          typename G1::Node minNode = n;
+          for(typename G1::NodeIt n1(_g1,minNode); n1!=INVALID; ++n1)
+            if(!added[n1] &&
+               (labelTmp1[_intLabels1[minNode]]>
+                labelTmp1[_intLabels1[n1]]||(dm1[minNode]<dm1[n1]&&
+                                             labelTmp1[_intLabels1[minNode]]==
+                                             labelTmp1[_intLabels1[n1]])))
+              minNode=n1;
+          processBFSLevel(minNode,orderIndex,dm1,added);
+        }
+        else
+          ++n;
+      }
+      for(unsigned int i = 0; i < labelTmp1.size(); ++i)
+        labelTmp1[i]=0;
+    }
+
+
+    template<MappingType MT>
+    bool extMatch(){
+      while(_depth>=0) {
+        //there is no node in g1, which has not pair in g2.
+        if(_depth==static_cast<int>(order.size())) {
+          --_depth;
+          return true;
+        }
+        typename G1::Node& nodeOfDepth = order[_depth];
+        const typename G2::Node& pairOfNodeOfDepth = _mapping[nodeOfDepth];
+        typename G2::IncEdgeIt &edgeItOfDepth = currEdgeIts[_depth];
+        //the node of g2, which neighbours are the candidates for
+        //the pair of order[_depth]
+        typename G2::Node currPNode;
+        if(edgeItOfDepth==INVALID){
+          typename G1::IncEdgeIt fstMatchedE(_g1,nodeOfDepth);
+          //if _mapping[order[_depth]]!=INVALID, we dont need
+          //fstMatchedE
+          if(pairOfNodeOfDepth==INVALID)
+            for(; fstMatchedE!=INVALID &&
+                  _mapping[_g1.oppositeNode(nodeOfDepth,
+                                            fstMatchedE)]==INVALID;
+                ++fstMatchedE); //find fstMatchedE, it could be preprocessed
+          if(fstMatchedE==INVALID||pairOfNodeOfDepth!=INVALID) {
+            //We found no covered neighbours, this means
+            //the graph is not connected(or _depth==0).  Each
+            //uncovered(and there are some other properties due
+            //to the spec. problem types) node of g2 is
+            //candidate.  We can read the iterator of the last
+            //tried node from the match if it is not the first
+            //try(match[nodeOfDepth]!=INVALID)
+            typename G2::NodeIt n2(_g2);
+            //if it's not the first try
+            if(pairOfNodeOfDepth!=INVALID) {
+              n2=++typename G2::NodeIt(_g2,pairOfNodeOfDepth);
+              subPair(nodeOfDepth,pairOfNodeOfDepth);
+            }
+            for(; n2!=INVALID; ++n2)
+              if(MT!=SUBGRAPH) {
+                if(_conn[n2]==0&&feas<MT>(nodeOfDepth,n2))
+                  break;
+              }
+              else if(_conn[n2]>=0&&feas<MT>(nodeOfDepth,n2))
+                break;
+            // n2 is the next candidate
+            if(n2!=INVALID) {
+              addPair(nodeOfDepth,n2);
+              ++_depth;
+            }
+            else // there are no more candidates
+              --_depth;
+            continue;
+          }
+          else{
+            currPNode=_mapping[_g1.oppositeNode(nodeOfDepth,
+                                                fstMatchedE)];
+            edgeItOfDepth=typename G2::IncEdgeIt(_g2,currPNode);
+          }
+        }
+        else{
+          currPNode=_g2.oppositeNode(pairOfNodeOfDepth,
+                                     edgeItOfDepth);
+          subPair(nodeOfDepth,pairOfNodeOfDepth);
+          ++edgeItOfDepth;
+        }
+        for(; edgeItOfDepth!=INVALID; ++edgeItOfDepth) {
+          const typename G2::Node currNode =
+            _g2.oppositeNode(currPNode, edgeItOfDepth);
+          if(_conn[currNode]>0&&feas<MT>(nodeOfDepth,currNode)) {
+            addPair(nodeOfDepth,currNode);
+            break;
+          }
+        }
+        edgeItOfDepth==INVALID?--_depth:++_depth;
+      }
+      return false;
+    }
+
+    //calc. the lookup table for cutting the searchtree
+    void setRNew1tRInOut1t(){
+      typename G1::template NodeMap<int> tmp(_g1,0);
+      for(unsigned int i=0; i<order.size(); ++i) {
+        tmp[order[i]]=-1;
+        for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1) {
+          const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
+          if(tmp[currNode]>0)
+            ++labelTmp1[_intLabels1[currNode]];
+          else if(tmp[currNode]==0)
+            ++labelTmp2[_intLabels1[currNode]];
+        }
+        //labelTmp1[i]=number of neightbours with label i in set rInOut
+        //labelTmp2[i]=number of neightbours with label i in set rNew
+        for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1) {
+          const int& currIntLabel = _intLabels1[_g1.oppositeNode(order[i],e1)];
+          if(labelTmp1[currIntLabel]>0) {
+            rInOutLabels1[order[i]]
+              .push_back(std::make_pair(currIntLabel,
+                                        labelTmp1[currIntLabel]));
+            labelTmp1[currIntLabel]=0;
+          }
+          else if(labelTmp2[currIntLabel]>0) {
+            rNewLabels1[order[i]].
+              push_back(std::make_pair(currIntLabel,labelTmp2[currIntLabel]));
+            labelTmp2[currIntLabel]=0;
+          }
+        }
+
+        for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1) {
+          int& tmpCurrNode=tmp[_g1.oppositeNode(order[i],e1)];
+          if(tmpCurrNode!=-1)
+            ++tmpCurrNode;
+        }
+      }
+    }
+
+    int getMaxLabel() const{
+      int m=-1;
+      for(typename G1::NodeIt n1(_g1); n1!=INVALID; ++n1) {
+        const int& currIntLabel = _intLabels1[n1];
+        if(currIntLabel>m)
+          m=currIntLabel;
+      }
+      for(typename G2::NodeIt n2(_g2); n2!=INVALID; ++n2) {
+        const int& currIntLabel = _intLabels2[n2];
+        if(currIntLabel>m)
+          m=currIntLabel;
+      }
+      return m;
+    }
+
+  public:
+    ///Constructor
+
+    ///Constructor.
+    ///\param g1 The graph to be embedded.
+    ///\param g2 The graph \e g1 will be embedded into.
+    ///\param m The type of the NodeMap storing the mapping.
+    ///By default, it is G1::NodeMap<G2::Node>
+    ///\param intLabel1 The NodeMap storing the integer node labels of G1.
+    ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
+    ///different labels.
+    ///\param intLabel1 The NodeMap storing the integer node labels of G2.
+    ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
+    ///different labels.
+    Vf2pp(const G1 &g1, const G2 &g2,M &m, M1 &intLabels1, M2 &intLabels2) :
+      _depth(0), _conn(g2,0), _mapping(m), order(countNodes(g1),INVALID),
+      currEdgeIts(countNodes(g1),INVALID), _g1(g1), _g2(g2), rNewLabels1(_g1),
+      rInOutLabels1(_g1), _intLabels1(intLabels1) ,_intLabels2(intLabels2),
+      maxLabel(getMaxLabel()), labelTmp1(maxLabel+1),labelTmp2(maxLabel+1),
+      _mapping_type(SUBGRAPH), _deallocMappingAfterUse(0),
+      _deallocLabelsAfterUse(0)
+    {
+      setOrder();
+      setRNew1tRInOut1t();
+
+      //reset mapping
+      for(typename G1::NodeIt n(g1);n!=INVALID;++n)
+        m[n]=INVALID;
+    }
+
+    ///Destructor
+
+    ///Destructor.
+    ///
+    ~Vf2pp()
+    {
+      if(_deallocMappingAfterUse)
+        delete &_mapping;
+      if(_deallocLabelsAfterUse) {
+        delete &_intLabels1;
+        delete &_intLabels2;
+      }
+    }
+
+    ///Returns the current mapping type.
+
+    ///Returns the current mapping type.
+    ///
+    MappingType mappingType() const
+    {
+      return _mapping_type;
+    }
+
+    ///Sets the mapping type
+
+    ///Sets the mapping type.
+    ///
+    ///The mapping type is set to \ref SUBGRAPH by default.
+    ///
+    ///\sa See \ref MappingType for the possible values.
+    void mappingType(MappingType m_type)
+    {
+      _mapping_type = m_type;
+    }
+
+    ///Finds a mapping.
+
+    ///This method finds a mapping from g1 into g2 according to the mapping
+    ///type set by \ref mappingType(MappingType) "mappingType()".
+    ///
+    ///By subsequent calls, it returns all possible mappings one-by-one.
+    ///
+    ///\retval true if a mapping is found.
+    ///\retval false if there is no (more) mapping.
+    bool find()
+    {
+      switch(_mapping_type)
+        {
+        case SUBGRAPH:
+          return extMatch<SUBGRAPH>();
+        case INDUCED:
+          return extMatch<INDUCED>();
+        case ISOMORPH:
+          return extMatch<ISOMORPH>();
+        default:
+          return false;
+        }
+    }
+  };
+
+  template<typename G1, typename G2>
+  class Vf2ppWizardBase {
+  protected:
+    typedef G1 Graph1;
+    typedef G2 Graph2;
+
+    const G1 &_g1;
+    const G2 &_g2;
+
+    MappingType _mapping_type;
+
+    typedef typename G1::template NodeMap<typename G2::Node> Mapping;
+    bool _local_mapping;
+    void *_mapping;
+    void createMapping() {
+      _mapping = new Mapping(_g1);
+    }
+
+    bool _local_nodeLabels;
+    typedef typename G1::template NodeMap<int> NodeLabels1;
+    typedef typename G2::template NodeMap<int> NodeLabels2;
+    void *_nodeLabels1, *_nodeLabels2;
+    void createNodeLabels() {
+      _nodeLabels1 = new NodeLabels1(_g1,0);
+      _nodeLabels2 = new NodeLabels2(_g2,0);
+    }
+
+    Vf2ppWizardBase(const G1 &g1,const G2 &g2)
+      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH),
+        _local_mapping(1), _local_nodeLabels(1) { }
+  };
+
+
+  /// \brief Auxiliary class for the function-type interface of %VF2
+  /// Plus Plus algorithm.
+  ///
+  /// This auxiliary class implements the named parameters of
+  /// \ref vf2pp() "function-type interface" of \ref Vf2pp algorithm.
+  ///
+  /// \warning This class is not to be used directly.
+  ///
+  /// \tparam TR The traits class that defines various types used by the
+  /// algorithm.
+  template<typename TR>
+  class Vf2ppWizard : public TR {
+    typedef TR Base;
+    typedef typename TR::Graph1 Graph1;
+    typedef typename TR::Graph2 Graph2;
+    typedef typename TR::Mapping Mapping;
+    typedef typename TR::NodeLabels1 NodeLabels1;
+    typedef typename TR::NodeLabels2 NodeLabels2;
+
+    using TR::_g1;
+    using TR::_g2;
+    using TR::_mapping_type;
+    using TR::_mapping;
+    using TR::_nodeLabels1;
+    using TR::_nodeLabels2;
+
+  public:
+    ///Constructor
+    Vf2ppWizard(const Graph1 &g1,const Graph2 &g2) : Base(g1,g2) { }
+
+    ///Copy constructor
+    Vf2ppWizard(const Base &b) : Base(b) {}
+
+
+    template<typename T>
+    struct SetMappingBase : public Base {
+      typedef T Mapping;
+      SetMappingBase(const Base &b) : Base(b) { }
+    };
+
+    ///\brief \ref named-templ-param "Named parameter" for setting
+    ///the mapping.
+    ///
+    ///\ref named-templ-param "Named parameter" function for setting
+    ///the map that stores the found embedding.
+    template<typename T>
+    Vf2ppWizard< SetMappingBase<T> > mapping(const T &t) {
+      Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
+      Base::_local_mapping = 0;
+      return Vf2ppWizard<SetMappingBase<T> >(*this);
+    }
+
+    template<typename NL1, typename NL2>
+    struct SetNodeLabelsBase : public Base {
+      typedef NL1 NodeLabels1;
+      typedef NL2 NodeLabels2;
+      SetNodeLabelsBase(const Base &b) : Base(b) { }
+    };
+
+    ///\brief \ref named-templ-param "Named parameter" for setting the
+    ///node labels.
+    ///
+    ///\ref named-templ-param "Named parameter" function for setting
+    ///the node labels.
+    ///
+    ///\param nodeLabels1 A \ref concepts::ReadMap "readable node map"
+    ///of g1 with integer value. In case of K different labels, the labels
+    ///must be the {0,1,..,K-1} numbers.
+    ///\param nodeLabels2 A \ref concepts::ReadMap "readable node map"
+    ///of g2 with integer value. In case of K different labels, the labels
+    ///must be the {0,1,..,K-1} numbers.
+    template<typename NL1, typename NL2>
+    Vf2ppWizard< SetNodeLabelsBase<NL1,NL2> >
+    nodeLabels(const NL1 &nodeLabels1, const NL2 &nodeLabels2) {
+      Base::_local_nodeLabels = 0;
+      Base::_nodeLabels1=
+        reinterpret_cast<void*>(const_cast<NL1*>(&nodeLabels1));
+      Base::_nodeLabels2=
+        reinterpret_cast<void*>(const_cast<NL2*>(&nodeLabels2));
+      return Vf2ppWizard<SetNodeLabelsBase<NL1,NL2> >
+        (SetNodeLabelsBase<NL1,NL2>(*this));
+    }
+
+
+    ///\brief \ref named-templ-param "Named parameter" for setting
+    ///the mapping type.
+    ///
+    ///\ref named-templ-param "Named parameter" for setting
+    ///the mapping type.
+    ///
+    ///The mapping type is set to \ref SUBGRAPH by default.
+    ///
+    ///\sa See \ref MappingType for the possible values.
+    Vf2ppWizard<Base> &mappingType(MappingType m_type) {
+      _mapping_type = m_type;
+      return *this;
+    }
+
+    ///\brief \ref named-templ-param "Named parameter" for setting
+    ///the mapping type to \ref INDUCED.
+    ///
+    ///\ref named-templ-param "Named parameter" for setting
+    ///the mapping type to \ref INDUCED.
+    Vf2ppWizard<Base> &induced() {
+      _mapping_type = INDUCED;
+      return *this;
+    }
+
+    ///\brief \ref named-templ-param "Named parameter" for setting
+    ///the mapping type to \ref ISOMORPH.
+    ///
+    ///\ref named-templ-param "Named parameter" for setting
+    ///the mapping type to \ref ISOMORPH.
+    Vf2ppWizard<Base> &iso() {
+      _mapping_type = ISOMORPH;
+      return *this;
+    }
+
+    ///Runs the %VF2 Plus Plus algorithm.
+
+    ///This method runs the VF2 Plus Plus algorithm.
+    ///
+    ///\retval true if a mapping is found.
+    ///\retval false if there is no mapping.
+    bool run() {
+      if(Base::_local_mapping)
+        Base::createMapping();
+      if(Base::_local_nodeLabels)
+        Base::createNodeLabels();
+
+      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >
+        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
+            *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
+            *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
+
+      alg.mappingType(_mapping_type);
+
+      const bool ret = alg.find();
+
+      if(Base::_local_nodeLabels) {
+        delete reinterpret_cast<NodeLabels1*>(_nodeLabels1);
+        delete reinterpret_cast<NodeLabels2*>(_nodeLabels2);
+      }
+      if(Base::_local_mapping)
+        delete reinterpret_cast<Mapping*>(_mapping);
+
+      return ret;
+    }
+
+    ///Get a pointer to the generated Vf2pp object.
+
+    ///Gives a pointer to the generated Vf2pp object.
+    ///
+    ///\return Pointer to the generated Vf2pp object.
+    ///\warning Don't forget to delete the referred Vf2pp object after use.
+    Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >*
+    getPtrToVf2ppObject(){
+      if(Base::_local_mapping)
+        Base::createMapping();
+      if(Base::_local_nodeLabels)
+        Base::createNodeLabels();
+
+      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >* ptr =
+        new Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2>
+        (_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
+         *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
+         *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
+      ptr->mappingType(_mapping_type);
+      if(Base::_local_mapping)
+        ptr->_deallocMappingAfterUse=true;
+      if(Base::_local_nodeLabels)
+        ptr->_deallocLabelMapsAfterUse=true;
+
+      return ptr;
+    }
+
+    ///Counts the number of mappings.
+
+    ///This method counts the number of mappings.
+    ///
+    /// \return The number of mappings.
+    int count() {
+      if(Base::_local_mapping)
+        Base::createMapping();
+      if(Base::_local_nodeLabels)
+        Base::createNodeLabels();
+
+      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2>
+        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
+            *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
+            *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
+
+      alg.mappingType(_mapping_type);
+
+      int ret = 0;
+      while(alg.find())
+        ++ret;
+
+      if(Base::_local_nodeLabels) {
+        delete reinterpret_cast<NodeLabels1*>(_nodeLabels1);
+        delete reinterpret_cast<NodeLabels2*>(_nodeLabels2);
+      }
+      if(Base::_local_mapping)
+        delete reinterpret_cast<Mapping*>(_mapping);
+
+      return ret;
+    }
+  };
+
+
+  ///Function-type interface for VF2 Plus Plus algorithm.
+
+  /// \ingroup graph_isomorphism
+  ///Function-type interface for VF2 Plus Plus algorithm.
+  ///
+  ///This function has several \ref named-func-param "named parameters"
+  ///declared as the members of class \ref Vf2ppWizard.
+  ///The following examples show how to use these parameters.
+  ///\code
+  ///  ListGraph::NodeMap<ListGraph::Node> m(g);
+  ///  // Find an embedding of graph g1 into graph g2
+  ///  vf2pp(g1,g2).mapping(m).run();
+  ///
+  ///  // Check whether graphs g1 and g2 are isomorphic
+  ///  bool is_iso = vf2pp(g1,g2).iso().run();
+  ///
+  ///  // Count the number of isomorphisms
+  ///  int num_isos = vf2pp(g1,g2).iso().count();
+  ///
+  ///  // Iterate through all the induced subgraph mappings
+  ///  //of graph g1 into g2 using the labels c1 and c2
+  ///  auto* myVf2pp = vf2pp(g1,g2).mapping(m).nodeLabels(c1,c2)
+  ///  .induced().getPtrToVf2Object();
+  ///  while(myVf2pp->find()){
+  ///    //process the current mapping m
+  ///  }
+  ///  delete myVf22pp;
+  ///\endcode
+  ///\warning Don't forget to put the \ref Vf2ppWizard::run() "run()",
+  ///\ref Vf2ppWizard::count() "count()" or
+  ///the \ref Vf2ppWizard::getPtrToVf2ppObject() "getPtrToVf2ppObject()"
+  ///to the end of the expression.
+  ///\sa Vf2ppWizard
+  ///\sa Vf2pp
+  template<class G1, class G2>
+  Vf2ppWizard<Vf2ppWizardBase<G1,G2> > vf2pp(const G1 &g1, const G2 &g2) {
+    return Vf2ppWizard<Vf2ppWizardBase<G1,G2> >(g1,g2);
+  }
+
+}
+
+#endif
+
Index: test/CMakeLists.txt
===================================================================
--- test/CMakeLists.txt	(revision 1141)
+++ test/CMakeLists.txt	(revision 1186)
@@ -56,4 +56,5 @@
   unionfind_test
   vf2_test
+  vf2pp_test
 )
 
Index: test/bellman_ford_test.cc
===================================================================
--- test/bellman_ford_test.cc	(revision 1162)
+++ test/bellman_ford_test.cc	(revision 1131)
@@ -101,5 +101,4 @@
     pp = const_bf_test.negativeCycle();
 
-#ifdef LEMON_CXX11
     for (BF::ActiveIt it(const_bf_test); it != INVALID; ++it) {}
     for (auto n: const_bf_test.activeNodes()) { ::lemon::ignore_unused_variable_warning(n); }
@@ -107,5 +106,4 @@
       ::lemon::ignore_unused_variable_warning(n);
     }
-#endif
   }
   {
Index: test/max_flow_test.cc
===================================================================
--- test/max_flow_test.cc	(revision 1178)
+++ test/max_flow_test.cc	(revision 1092)
@@ -27,5 +27,4 @@
 #include <lemon/lgf_reader.h>
 #include <lemon/elevator.h>
-#include <lemon/tolerance.h>
 
 using namespace lemon;
@@ -67,35 +66,4 @@
   "target 8\n";
 
-char test_lgf_float[] =
-  "@nodes\n"
-  "label\n"
-  "0\n"
-  "1\n"
-  "2\n"
-  "3\n"
-  "4\n"
-  "5\n"
-  "6\n"
-  "7\n"
-  "8\n"
-  "9\n"
-  "@arcs\n"
-  "      capacity\n"
-  "0 1 0.1\n"
-  "0 2 0.1\n"
-  "0 3 0.1\n"
-  "1 4 0.1\n"
-  "2 4 0.1\n"
-  "3 4 0.1\n"
-  "4 5 0.3\n"
-  "5 6 0.1\n"
-  "5 7 0.1\n"
-  "5 8 0.1\n"
-  "6 9 0.1\n"
-  "7 9 0.1\n"
-  "8 9 0.1\n"
-  "@attributes\n"
-  "source 0\n"
-  "target 9\n";
 
 // Checks the general interface of a max flow algorithm
@@ -198,4 +166,6 @@
   typedef concepts::Digraph Digraph;
   typedef concepts::ReadMap<Digraph::Arc, Value> CapMap;
+  typedef Elevator<Digraph, Digraph::Node> Elev;
+  typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev;
 
   Digraph g;
@@ -215,11 +185,11 @@
 
 template <typename T>
-T cutValue(const SmartDigraph& g,
-           const SmartDigraph::NodeMap<bool>& cut,
-           const SmartDigraph::ArcMap<T>& cap) {
-
-  T c = 0;
-  for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {
-    if (cut[g.source(e)] && !cut[g.target(e)]) c += cap[e];
+T cutValue (const SmartDigraph& g,
+              const SmartDigraph::NodeMap<bool>& cut,
+              const SmartDigraph::ArcMap<T>& cap) {
+
+  T c=0;
+  for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) {
+    if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e];
   }
   return c;
@@ -230,9 +200,8 @@
                const SmartDigraph::ArcMap<T>& flow,
                const SmartDigraph::ArcMap<T>& cap,
-               SmartDigraph::Node s, SmartDigraph::Node t,
-               const Tolerance<T>& tol) {
+               SmartDigraph::Node s, SmartDigraph::Node t) {
 
   for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {
-    if (tol.negative(flow[e]) || tol.less(cap[e], flow[e])) return false;
+    if (flow[e] < 0 || flow[e] > cap[e]) return false;
   }
 
@@ -246,27 +215,26 @@
       sum -= flow[e];
     }
-    if (tol.nonZero(sum)) return false;
+    if (sum != 0) return false;
   }
   return true;
 }
 
-void checkInitPreflow()
+void initFlowTest()
 {
   DIGRAPH_TYPEDEFS(SmartDigraph);
 
   SmartDigraph g;
-  SmartDigraph::ArcMap<int> cap(g), iflow(g);
-  Node s = g.addNode(); Node t = g.addNode();
-  Node n1 = g.addNode(); Node n2 = g.addNode();
+  SmartDigraph::ArcMap<int> cap(g),iflow(g);
+  Node s=g.addNode(); Node t=g.addNode();
+  Node n1=g.addNode(); Node n2=g.addNode();
   Arc a;
-  a = g.addArc(s, n1); cap[a] = 20; iflow[a] = 20;
-  a = g.addArc(n1, n2); cap[a] = 10; iflow[a] = 0;
-  a = g.addArc(n2, t); cap[a] = 20; iflow[a] = 0;
-
-  Preflow<SmartDigraph> pre(g, cap, s, t);
+  a=g.addArc(s,n1); cap[a]=20; iflow[a]=20;
+  a=g.addArc(n1,n2); cap[a]=10; iflow[a]=0;
+  a=g.addArc(n2,t); cap[a]=20; iflow[a]=0;
+
+  Preflow<SmartDigraph> pre(g,cap,s,t);
   pre.init(iflow);
   pre.startFirstPhase();
-
-  check(pre.flowValue() == 10, "Incorrect max flow value.");
+  check(pre.flowValue() == 10, "The incorrect max flow value.");
   check(pre.minCut(s), "Wrong min cut (Node s).");
   check(pre.minCut(n1), "Wrong min cut (Node n1).");
@@ -276,5 +244,5 @@
 
 template <typename MF, typename SF>
-void checkMaxFlowAlg(const char *input_lgf,  typename MF::Value expected) {
+void checkMaxFlowAlg() {
   typedef SmartDigraph Digraph;
   DIGRAPH_TYPEDEFS(Digraph);
@@ -285,14 +253,12 @@
   typedef BoolNodeMap CutMap;
 
-  Tolerance<Value> tol;
-
   Digraph g;
   Node s, t;
   CapMap cap(g);
-  std::istringstream input(input_lgf);
-  DigraphReader<Digraph>(g, input)
+  std::istringstream input(test_lgf);
+  DigraphReader<Digraph>(g,input)
       .arcMap("capacity", cap)
-      .node("source", s)
-      .node("target", t)
+      .node("source",s)
+      .node("target",t)
       .run();
 
@@ -300,7 +266,5 @@
   max_flow.run();
 
-  check(!tol.different(expected, max_flow.flowValue()),
-        "Incorrect max flow value.");
-  check(checkFlow(g, max_flow.flowMap(), cap, s, t, tol),
+  check(checkFlow(g, max_flow.flowMap(), cap, s, t),
         "The flow is not feasible.");
 
@@ -309,10 +273,13 @@
   Value min_cut_value = cutValue(g, min_cut, cap);
 
-  check(!tol.different(expected, min_cut_value),
-        "Incorrect min cut value.");
+  check(max_flow.flowValue() == min_cut_value,
+        "The max flow value is not equal to the min cut value.");
 
   FlowMap flow(g);
-  for (ArcIt e(g); e != INVALID; ++e) flow[e] = 13 * max_flow.flowMap()[e];
-  for (ArcIt e(g); e != INVALID; ++e) cap[e] = 17 * cap[e];
+  for (ArcIt e(g); e != INVALID; ++e) flow[e] = max_flow.flowMap()[e];
+
+  Value flow_value = max_flow.flowValue();
+
+  for (ArcIt e(g); e != INVALID; ++e) cap[e] = 2 * cap[e];
   max_flow.init(flow);
 
@@ -323,12 +290,11 @@
   min_cut_value = cutValue(g, min_cut1, cap);
 
-  check(!tol.different(17 * expected, max_flow.flowValue()),
-        "Incorrect max flow value.");
-  check(!tol.different(17 * expected, min_cut_value),
-        "Incorrect min cut value.");
+  check(max_flow.flowValue() == min_cut_value &&
+        min_cut_value == 2 * flow_value,
+        "The max flow value or the min cut value is wrong.");
 
   SF::startSecondPhase(max_flow);       // start second phase of the algorithm
 
-  check(checkFlow(g, max_flow.flowMap(), cap, s, t, tol),
+  check(checkFlow(g, max_flow.flowMap(), cap, s, t),
         "The flow is not feasible.");
 
@@ -337,8 +303,8 @@
   min_cut_value = cutValue(g, min_cut2, cap);
 
-  check(!tol.different(17 * expected, max_flow.flowValue()),
-        "Incorrect max flow value.");
-  check(!tol.different(17 * expected, min_cut_value),
-        "Incorrect min cut value.");
+  check(max_flow.flowValue() == min_cut_value &&
+        min_cut_value == 2 * flow_value,
+        "The max flow value or the min cut value was not doubled");
+
 
   max_flow.flowMap(flow);
@@ -359,7 +325,7 @@
   CutMap min_cut3(g);
   max_flow.minCutMap(min_cut3);
-  min_cut_value = cutValue(g, min_cut3, cap);
-
-  check(!tol.different(max_flow.flowValue(), min_cut_value),
+  min_cut_value=cutValue(g, min_cut3, cap);
+
+  check(max_flow.flowValue() == min_cut_value,
         "The max flow value or the min cut value is wrong.");
 }
@@ -414,26 +380,15 @@
   typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<int> > PType1;
   typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<float> > PType2;
-  typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<double> > PType3;
-
-  checkMaxFlowAlg<PType1, PreflowStartFunctions<PType1> >(test_lgf, 13);
-  checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >(test_lgf, 13);
-  checkMaxFlowAlg<PType3, PreflowStartFunctions<PType3> >(test_lgf, 13);
-
-  checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >(test_lgf_float, 0.3f);
-  checkMaxFlowAlg<PType3, PreflowStartFunctions<PType3> >(test_lgf_float, 0.3);
-
-  checkInitPreflow();
+  checkMaxFlowAlg<PType1, PreflowStartFunctions<PType1> >();
+  checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >();
+  initFlowTest();
 
   // Check EdmondsKarp
   typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<int> > EKType1;
   typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<float> > EKType2;
-  typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<double> > EKType3;
-
-  checkMaxFlowAlg<EKType1, GeneralStartFunctions<EKType1> >(test_lgf, 13);
-  checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >(test_lgf, 13);
-  checkMaxFlowAlg<EKType3, GeneralStartFunctions<EKType3> >(test_lgf, 13);
-
-  checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >(test_lgf_float, 0.3f);
-  checkMaxFlowAlg<EKType3, GeneralStartFunctions<EKType3> >(test_lgf_float, 0.3);
+  checkMaxFlowAlg<EKType1, GeneralStartFunctions<EKType1> >();
+  checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >();
+
+  initFlowTest();
 
   return 0;
Index: test/planarity_test.cc
===================================================================
--- test/planarity_test.cc	(revision 1182)
+++ test/planarity_test.cc	(revision 798)
@@ -31,36 +31,6 @@
 using namespace lemon::dim2;
 
-const int lgfn = 8;
+const int lgfn = 4;
 const std::string lgf[lgfn] = {
-  "@nodes\n"
-  "label\n"
-  "@edges\n"
-  "     label\n",
-
-  "@nodes\n"
-  "label\n"
-  "0\n"
-  "@edges\n"
-  "     label\n",
-
-  "@nodes\n"
-  "label\n"
-  "0\n"
-  "1\n"
-  "@edges\n"
-  "     label\n"
-  "0 1  0\n",
-
-  "@nodes\n"
-  "label\n"
-  "0\n"
-  "1\n"
-  "2\n"
-  "@edges\n"
-  "     label\n"
-  "0 1  0\n"
-  "1 2  1\n"
-  "2 0  2\n",
-
   "@nodes\n"
   "label\n"
@@ -167,9 +137,6 @@
     }
   }
-
-  if (face_num != 0) {
-    check(face_num + countNodes(graph) - countConnectedComponents(graph) ==
-          countEdges(graph) + 1, "Euler test does not passed");
-  }
+  check(face_num + countNodes(graph) - countConnectedComponents(graph) ==
+        countEdges(graph) + 1, "Euler test does not passed");
 }
 
@@ -279,27 +246,11 @@
       checkEmbedding(graph, pe);
 
-      {
-        PlanarDrawing<Graph> pd(graph);
-        pd.run(pe.embeddingMap());
-        checkDrawing(graph, pd);
-      }
-
-      {
-        PlanarDrawing<Graph> pd(graph);
-        pd.run();
-        checkDrawing(graph, pd);
-      }
-
-      {
-        PlanarColoring<Graph> pc(graph);
-        pc.runFiveColoring(pe.embeddingMap());
-        checkColoring(graph, pc, 5);
-      }
-
-      {
-        PlanarColoring<Graph> pc(graph);
-        pc.runFiveColoring();
-        checkColoring(graph, pc, 5);
-      }
+      PlanarDrawing<Graph> pd(graph);
+      pd.run(pe.embeddingMap());
+      checkDrawing(graph, pd);
+
+      PlanarColoring<Graph> pc(graph);
+      pc.runFiveColoring(pe.embeddingMap());
+      checkColoring(graph, pc, 5);
 
     } else {
Index: test/random_test.cc
===================================================================
--- test/random_test.cc	(revision 1163)
+++ test/random_test.cc	(revision 440)
@@ -22,31 +22,4 @@
 int seed_array[] = {1, 2};
 
-int rnd_seq32[] = {
-2732, 43567, 42613, 52416, 45891, 21243, 30403, 32103, 
-62501, 33003, 12172, 5192, 32511, 50057, 43723, 7813, 
-23720, 35343, 6637, 30280, 44566, 31019, 18898, 33867, 
-5994, 1688, 11513, 59011, 48056, 25544, 39168, 25365, 
-17530, 8366, 27063, 49861, 55169, 63848, 11863, 49608
-};
-int rnd_seq64[] = {
-56382, 63883, 59577, 64750, 9644, 59886, 57647, 18152, 
-28520, 64078, 17818, 49294, 26424, 26697, 53684, 19209, 
-35404, 12121, 12837, 11827, 32156, 58333, 62553, 7907, 
-64427, 39399, 21971, 48789, 46981, 15716, 53335, 65256, 
-12999, 15308, 10906, 42162, 47587, 43006, 53921, 18716
-};
-
-void seq_test() {
-  for(int i=0;i<5;i++) {
-    lemon::Random32 r32(i);
-    lemon::Random64 r64(i);
-    for(int j=0;j<8;j++) {
-      check(r32[65536]==rnd_seq32[i*8+j], "Wrong random sequence");
-      check(r64[65536]==rnd_seq64[i*8+j], "Wrong random sequence");
-    }
-  }
-}
-
-
 int main()
 {
@@ -64,5 +37,4 @@
                   (sizeof(seed_array) / sizeof(seed_array[0])));
 
-  seq_test();
   return 0;
 }
Index: test/vf2_test.cc
===================================================================
--- test/vf2_test.cc	(revision 1143)
+++ test/vf2_test.cc	(revision 1186)
@@ -3,5 +3,5 @@
  * This file is a part of LEMON, a generic C++ optimization library.
  *
- * Copyright (C) 2015
+ * Copyright (C) 2015-2017
  * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
  *
@@ -184,5 +184,7 @@
 class EqComparable {
 public:
-  bool operator==(const EqComparable&) { return false; }
+  bool operator==(const EqComparable&) {
+    return false;
+  }
 };
 
@@ -190,10 +192,11 @@
 class EqClass {
 public:
-  bool operator()(A, B) { return false; }
+  bool operator()(A, B){
+    return false;
+  }
 };
 
 template<class G1,class G2>
-void checkVf2Compile()
-{
+void checkVf2Compile() {
   G1 g;
   G2 h;
@@ -206,6 +209,13 @@
   succ = vf2(g,h).iso().run();
   succ = vf2(g,h).mapping(r).run();
+
+  Vf2<G1,G2,concepts::ReadWriteMap<typename G1::Node, typename G2::Node>,
+      EqClass<typename G1::Node,typename G2::Node> >
+    myVf2(g,h,r,EqClass<typename G1::Node,typename G2::Node>());
+  myVf2.find();
+
   succ = vf2(g,h).induced().mapping(r).run();
   succ = vf2(g,h).iso().mapping(r).run();
+
   concepts::ReadMap<typename G1::Node, EqComparable> l1;
   concepts::ReadMap<typename G2::Node, EqComparable> l2;
@@ -215,6 +225,5 @@
 }
 
-void justCompile()
-{
+void justCompile() {
   checkVf2Compile<concepts::Graph,concepts::Graph>();
   checkVf2Compile<concepts::Graph,SmartGraph>();
@@ -223,14 +232,12 @@
 
 template<class G1, class G2, class I>
-void checkSub(const G1 &g1, const G2 &g2, const I &i)
-{
+void checkSub(const G1 &g1, const G2 &g2, const I &i) {
   {
     std::set<typename G2::Node> image;
-    for(typename G1::NodeIt n(g1);n!=INVALID;++n)
-      {
-          check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
-          check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
-          image.insert(i[n]);
-      }
+    for(typename G1::NodeIt n(g1);n!=INVALID;++n){
+      check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
+      check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
+      image.insert(i[n]);
+    }
   }
   for(typename G1::EdgeIt e(g1);e!=INVALID;++e)
@@ -240,85 +247,76 @@
 
 template<class G1, class G2, class I>
-void checkInd(const G1 &g1, const G2 &g2, const I &i)
-  {
+void checkInd(const G1 &g1, const G2 &g2, const I &i) {
   std::set<typename G2::Node> image;
-  for(typename G1::NodeIt n(g1);n!=INVALID;++n)
-    {
+  for(typename G1::NodeIt n(g1);n!=INVALID;++n) {
     check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
     check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
     image.insert(i[n]);
-    }
+  }
   for(typename G1::NodeIt n(g1); n!=INVALID; ++n)
     for(typename G1::NodeIt m(g1); m!=INVALID; ++m)
-      if((findEdge(g1,n,m)==INVALID) != (findEdge(g2,i[n],i[m])==INVALID))
-        {
+      if((findEdge(g1,n,m)==INVALID) != (findEdge(g2,i[n],i[m])==INVALID)) {
         std::cout << "Wrong isomorphism: edge mismatch";
         exit(1);
-        }
-  }
+      }
+}
 
 template<class G1,class G2>
-int checkSub(const G1 &g1, const G2 &g2)
-{
+int checkSub(const G1 &g1, const G2 &g2) {
   typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
-  if(vf2(g1,g2).mapping(iso).run())
-    {
-      checkSub(g1,g2,iso);
-      return true;
-    }
-  else return false;
+  if(vf2(g1,g2).mapping(iso).run()) {
+    checkSub(g1,g2,iso);
+    return true;
+  }
+  else
+    return false;
 }
 
 template<class G1,class G2>
-int checkInd(const G1 &g1, const G2 &g2)
-{
+int checkInd(const G1 &g1, const G2 &g2) {
   typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
-  if(vf2(g1,g2).induced().mapping(iso).run())
-    {
-      checkInd(g1,g2,iso);
-      return true;
-    }
-  else return false;
+  if(vf2(g1,g2).induced().mapping(iso).run()) {
+    checkInd(g1,g2,iso);
+    return true;
+  }
+  else
+    return false;
 }
 
 template<class G1,class G2>
-int checkIso(const G1 &g1, const G2 &g2)
-{
+int checkIso(const G1 &g1, const G2 &g2) {
   typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
-  if(vf2(g1,g2).iso().mapping(iso).run())
-    {
-      check(countNodes(g1)==countNodes(g2),
-            "Wrong iso alg.: they are not isomophic.");
-      checkInd(g1,g2,iso);
-      return true;
-    }
-  else return false;
+  if(vf2(g1,g2).iso().mapping(iso).run()) {
+    check(countNodes(g1)==countNodes(g2),
+          "Wrong iso alg.: they are not isomophic.");
+    checkInd(g1,g2,iso);
+    return true;
+  }
+  else
+    return false;
 }
 
 template<class G1, class G2, class L1, class L2, class I>
 void checkLabel(const G1 &g1, const G2 &,
-                const L1 &l1, const L2 &l2,const I &i)
-{
+                const L1 &l1, const L2 &l2,const I &i) {
   for(typename G1::NodeIt n(g1);n!=INVALID;++n)
-    {
-      check(l1[n]==l2[i[n]],"Wrong isomorphism: label mismatch.");
-    }
+    check(l1[n]==l2[i[n]],"Wrong isomorphism: label mismatch.");
 }
 
 template<class G1,class G2,class L1,class L2>
-int checkSub(const G1 &g1, const G2 &g2, const L1 &l1, const L2 &l2)
-{
+int checkSub(const G1 &g1, const G2 &g2, const L1 &l1, const L2 &l2) {
   typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
-  if(vf2(g1,g2).nodeLabels(l1,l2).mapping(iso).run())
-    {
-      checkSub(g1,g2,iso);
-      checkLabel(g1,g2,l1,l2,iso);
-      return true;
-    }
-  else return false;
+  if(vf2(g1,g2).nodeLabels(l1,l2).mapping(iso).run()){
+    checkSub(g1,g2,iso);
+    checkLabel(g1,g2,l1,l2,iso);
+    return true;
+  }
+  else
+    return false;
 }
 
 int main() {
   make_graphs();
+  //   justCompile();
   check(checkSub(c5,petersen), "There should exist a C5->Petersen mapping.");
   check(!checkSub(c7,petersen),
Index: test/vf2pp_test.cc
===================================================================
--- test/vf2pp_test.cc	(revision 1186)
+++ test/vf2pp_test.cc	(revision 1186)
@@ -0,0 +1,392 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2015-2017
+ * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
+ *
+ * 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.
+ *
+ */
+
+#include <lemon/vf2pp.h>
+#include <lemon/concepts/digraph.h>
+#include <lemon/smart_graph.h>
+#include <lemon/lgf_reader.h>
+#include <lemon/concepts/maps.h>
+#include <lemon/maps.h>
+#include <lemon/list_graph.h>
+
+#include <test/test_tools.h>
+#include <sstream>
+
+using namespace lemon;
+
+char petersen_lgf[] =
+  "@nodes\n"
+  "label col1 col2\n"
+  "0 1 1\n"
+  "1 1 2\n"
+  "2 1 3\n"
+  "3 1 4\n"
+  "4 2 5\n"
+  "5 2 1\n"
+  "6 2 2\n"
+  "7 2 3\n"
+  "8 2 4\n"
+  "9 2 5\n"
+  "@arcs\n"
+  "     -\n"
+  "0 1\n"
+  "1 2\n"
+  "2 3\n"
+  "3 4\n"
+  "4 0\n"
+  "0 5\n"
+  "1 6\n"
+  "2 7\n"
+  "3 8\n"
+  "4 9\n"
+  "5 8\n"
+  "5 7\n"
+  "9 6\n"
+  "9 7\n"
+  "6 8\n";
+
+char c5_lgf[] =
+  "@nodes\n"
+  "label col\n"
+  "0 1\n"
+  "1 2\n"
+  "2 3\n"
+  "3 4\n"
+  "4 5\n"
+  "@arcs\n"
+  "     -\n"
+  "0 1\n"
+  "1 2\n"
+  "2 3\n"
+  "3 4\n"
+  "4 0\n";
+
+char c7_lgf[] =
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "@arcs\n"
+  "     -\n"
+  "0 1\n"
+  "1 2\n"
+  "2 3\n"
+  "3 4\n"
+  "4 5\n"
+  "5 6\n"
+  "6 0\n";
+
+char c10_lgf[] =
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "8\n"
+  "9\n"
+  "@arcs\n"
+  "     -\n"
+  "0 1\n"
+  "1 2\n"
+  "2 3\n"
+  "3 4\n"
+  "4 5\n"
+  "5 6\n"
+  "6 7\n"
+  "7 8\n"
+  "8 9\n"
+  "9 0\n";
+
+char p10_lgf[] =
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "8\n"
+  "9\n"
+  "@arcs\n"
+  "     -\n"
+  "0 1\n"
+  "1 2\n"
+  "2 3\n"
+  "3 4\n"
+  "4 5\n"
+  "5 6\n"
+  "6 7\n"
+  "7 8\n"
+  "8 9\n";
+
+SmartGraph petersen, c5, c7, c10, p10;
+SmartGraph::NodeMap<int> petersen_col1(petersen);
+SmartGraph::NodeMap<int> petersen_col2(petersen);
+SmartGraph::NodeMap<int> c5_col(c5);
+
+void make_graphs(){
+  std::stringstream ss(petersen_lgf);
+  graphReader(petersen, ss)
+    .nodeMap("col1",petersen_col1)
+    .nodeMap("col2",petersen_col2)
+    .run();
+
+  ss.clear();
+  ss.str("");
+  ss<<c5_lgf;
+
+  graphReader(c5, ss)
+    .nodeMap("col",c5_col)
+    .run();
+
+  ss.clear();
+  ss.str("");
+  ss<<c7_lgf;
+  graphReader(c7, ss).run();
+
+  ss.clear();
+  ss.str("");
+  ss<<c10_lgf;
+  graphReader(c10, ss).run();
+
+  ss.clear();
+  ss.str("");
+  ss<<p10_lgf;
+  graphReader(p10, ss).run();
+
+}
+
+class IntConvertible1{
+public:
+  operator int(){
+    return 0;
+  }
+};
+
+class IntConvertible2{
+public:
+  operator int(){
+    return 0;
+  }
+};
+
+template<class G1,class G2>
+void checkVf2Compile() {
+  G1 g;
+  G2 h;
+  concepts::ReadWriteMap<typename G1::Node, typename G2::Node> r;
+  bool succ;
+  ::lemon::ignore_unused_variable_warning(succ);
+
+  succ = vf2pp(g,h).run();
+  succ = vf2pp(g,h).induced().run();
+  succ = vf2pp(g,h).iso().run();
+  succ = vf2pp(g,h).mapping(r).run();
+  succ = vf2pp(g,h).induced().mapping(r).run();
+  succ = vf2pp(g,h).iso().mapping(r).run();
+
+
+  concepts::ReadMap<typename G1::Node, int> c1;
+  concepts::ReadMap<typename G2::Node, int> c2;
+  Vf2pp<G1,G2,concepts::ReadWriteMap<typename G1::Node, typename G2::Node>,
+        concepts::ReadMap<typename G1::Node, int>,
+        concepts::ReadMap<typename G2::Node, int> >
+    myVf2pp(g,h,r,c1,c2);
+  myVf2pp.find();
+
+  succ = vf2pp(g,h).nodeLabels(c1,c2).mapping(r).run();
+  succ = vf2pp(g,h).nodeLabels(c1,c2)
+    .mapping(r).run();
+
+
+  concepts::ReadMap<typename G1::Node, char> c1_c;
+  concepts::ReadMap<typename G2::Node, char> c2_c;
+  Vf2pp<G1,G2,concepts::ReadWriteMap<typename G1::Node, typename G2::Node>,
+        concepts::ReadMap<typename G1::Node, char>,
+        concepts::ReadMap<typename G2::Node, char> >
+    myVf2pp_c(g,h,r,c1_c,c2_c);
+  myVf2pp_c.find();
+
+  succ = vf2pp(g,h).nodeLabels(c1_c,c2_c).mapping(r).run();
+  succ = vf2pp(g,h).nodeLabels(c1_c,c2_c)
+    .mapping(r).run();
+
+
+  concepts::ReadMap<typename G1::Node, IntConvertible1> c1_IntConv;
+  concepts::ReadMap<typename G2::Node, IntConvertible2> c2_IntConv;
+  Vf2pp<G1,G2,concepts::ReadWriteMap<typename G1::Node, typename G2::Node>,
+        concepts::ReadMap<typename G1::Node, IntConvertible1>,
+        concepts::ReadMap<typename G2::Node, IntConvertible2> >
+    myVf2pp_IntConv(g,h,r,c1_IntConv,c2_IntConv);
+  myVf2pp_IntConv.find();
+
+  succ = vf2pp(g,h).nodeLabels(c1_IntConv,c2_IntConv).mapping(r).run();
+  succ = vf2pp(g,h).nodeLabels(c1_IntConv,c2_IntConv)
+    .mapping(r).run();
+}
+
+void justCompile() {
+  checkVf2Compile<concepts::Graph,concepts::Graph>();
+  checkVf2Compile<concepts::Graph,SmartGraph>();
+  checkVf2Compile<SmartGraph,concepts::Graph>();
+}
+
+template<class G1, class G2, class I>
+void checkSub(const G1 &g1, const G2 &g2, const I &i) {
+  {
+    std::set<typename G2::Node> image;
+    for(typename G1::NodeIt n(g1);n!=INVALID;++n){
+      check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
+      check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
+      image.insert(i[n]);
+    }
+  }
+  for(typename G1::EdgeIt e(g1);e!=INVALID;++e)
+    check(findEdge(g2,i[g1.u(e)],i[g1.v(e)])!=INVALID,
+          "Wrong isomorphism: missing edge(checkSub).");
+}
+
+template<class G1, class G2, class I>
+void checkInd(const G1 &g1, const G2 &g2, const I &i) {
+  std::set<typename G2::Node> image;
+  for(typename G1::NodeIt n(g1);n!=INVALID;++n) {
+    check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
+    check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
+    image.insert(i[n]);
+  }
+  for(typename G1::NodeIt n(g1); n!=INVALID; ++n)
+    for(typename G1::NodeIt m(g1); m!=INVALID; ++m)
+      if((findEdge(g1,n,m)==INVALID) != (findEdge(g2,i[n],i[m])==INVALID)) {
+        std::cout << "Wrong isomorphism: edge mismatch";
+        exit(1);
+      }
+}
+
+template<class G1,class G2>
+int checkSub(const G1 &g1, const G2 &g2) {
+  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
+  if(vf2pp(g1,g2).mapping(iso).run()){
+    checkSub(g1,g2,iso);
+    return true;
+  }
+  else
+    return false;
+}
+
+template<class G1,class G2>
+int checkInd(const G1 &g1, const G2 &g2) {
+  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
+  if(vf2pp(g1,g2).induced().mapping(iso).run()) {
+    checkInd(g1,g2,iso);
+    return true;
+  }
+  else
+    return false;
+}
+
+template<class G1,class G2>
+int checkIso(const G1 &g1, const G2 &g2) {
+  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
+  if(vf2pp(g1,g2).iso().mapping(iso).run()) {
+    check(countNodes(g1)==countNodes(g2),
+          "Wrong iso alg.: they are not isomophic.");
+    checkInd(g1,g2,iso);
+    return true;
+  }
+  else
+    return false;
+}
+
+template<class G1, class G2, class L1, class L2, class I>
+void checkLabel(const G1 &g1, const G2 &,
+                const L1 &l1, const L2 &l2,const I &i) {
+  for(typename G1::NodeIt n(g1);n!=INVALID;++n)
+    check(l1[n]==l2[i[n]],"Wrong isomorphism: label mismatch.");
+}
+
+template<class G1,class G2,class L1,class L2>
+int checkSub(const G1 &g1, const G2 &g2, const L1 &l1, const L2 &l2) {
+  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
+  if(vf2pp(g1,g2).nodeLabels(l1,l2).mapping(iso).run()) {
+    checkSub(g1,g2,iso);
+    checkLabel(g1,g2,l1,l2,iso);
+    return true;
+  }
+  else
+    return false;
+}
+
+int main() {
+  make_graphs();
+//   justCompile();
+  check(checkSub(c5,petersen), "There should exist a C5->Petersen mapping.");
+  check(!checkSub(c7,petersen),
+        "There should not exist a C7->Petersen mapping.");
+  check(checkSub(p10,petersen), "There should exist a P10->Petersen mapping.");
+  check(!checkSub(c10,petersen),
+        "There should not exist a C10->Petersen mapping.");
+  check(checkSub(petersen,petersen),
+        "There should exist a Petersen->Petersen mapping.");
+
+  check(checkInd(c5,petersen),
+        "There should exist a C5->Petersen spanned mapping.");
+  check(!checkInd(c7,petersen),
+        "There should exist a C7->Petersen spanned mapping.");
+  check(!checkInd(p10,petersen),
+        "There should not exist a P10->Petersen spanned mapping.");
+  check(!checkInd(c10,petersen),
+        "There should not exist a C10->Petersen spanned mapping.");
+  check(checkInd(petersen,petersen),
+        "There should exist a Petersen->Petersen spanned mapping.");
+
+  check(!checkSub(petersen,c10),
+        "There should not exist a Petersen->C10 mapping.");
+  check(checkSub(p10,c10),
+        "There should exist a P10->C10 mapping.");
+  check(!checkInd(p10,c10),
+        "There should not exist a P10->C10 spanned mapping.");
+  check(!checkSub(c10,p10),
+        "There should not exist a C10->P10 mapping.");
+
+  check(!checkIso(p10,c10),
+        "P10 and C10 are not isomorphic.");
+  check(checkIso(c10,c10),
+        "C10 and C10 are isomorphic.");
+
+  check(!vf2pp(p10,c10).iso().run(),
+        "P10 and C10 are not isomorphic.");
+  check(vf2pp(c10,c10).iso().run(),
+        "C10 and C10 are isomorphic.");
+
+  check(!checkSub(c5,petersen,c5_col,petersen_col1),
+        "There should exist a C5->Petersen mapping.");
+  check(checkSub(c5,petersen,c5_col,petersen_col2),
+        "There should exist a C5->Petersen mapping.");
+}
Index: tools/dimacs-solver.cc
===================================================================
--- tools/dimacs-solver.cc	(revision 1170)
+++ tools/dimacs-solver.cc	(revision 1093)
@@ -223,5 +223,4 @@
         throw IoError("Cannot open the file for writing", ap.files()[1]);
       }
-      // fall through
     case 1:
       input.open(ap.files()[0].c_str());
@@ -229,5 +228,4 @@
         throw IoError("File cannot be found", ap.files()[0]);
       }
-      // fall through
     case 0:
       break;
@@ -254,5 +252,4 @@
         case DimacsDescriptor::SP:
           std::cout << "sp";
-          break;
         case DimacsDescriptor::MAT:
           std::cout << "mat";
Index: tools/dimacs-to-lgf.cc
===================================================================
--- tools/dimacs-to-lgf.cc	(revision 1179)
+++ tools/dimacs-to-lgf.cc	(revision 584)
@@ -74,5 +74,4 @@
         throw IoError("Cannot open the file for writing", ap.files()[1]);
       }
-      // fall through
     case 1:
       input.open(ap.files()[0].c_str());
@@ -80,5 +79,4 @@
         throw IoError("File cannot be found", ap.files()[0]);
       }
-      // fall through
     case 0:
       break;
