Ticket #9420: trac_9420.patch

sage/graphs/base/c_graph.pxd

@@ -28 +28 @@
     cpdef add_arc(self, int u, int v)
     cpdef bint has_arc(self, int u, int v)
     cpdef del_all_arcs(self, int u, int v)
-    cdef int *adjacency_sequence(self, int n, int *vertices, int v)
+    cdef int *adjacency_sequence_in(self, int n, int *vertices, int v)
+    cdef int *adjacency_sequence_out(self, int n, int *vertices, int v)
     cpdef list all_arcs(self, int u, int v)
     cpdef list in_neighbors(self, int v)
     cpdef list out_neighbors(self, int u)

sage/graphs/base/c_graph.pyx

@@ -840 +840 @@
         """
         raise NotImplementedError()
 
-    cdef int *adjacency_sequence(self, int n, int *vertices, int v):
+    cdef int * adjacency_sequence_in(self, int n, int *vertices, int v):
         r""" 
         Returns the adjacency sequence corresponding to a list of vertices
-        and a vertex. See the function ``_test_adjacency_sequence()`` of
+        and a vertex. 
+
+        See the OUTPUT section for a formal definition.
+
+        See the function ``_test_adjacency_sequence()`` of
         ``dense_graph.pyx`` and ``sparse_graph.pyx`` for unit tests.
         
         INPUT: 
 
-        - ``n`` -- nonnegative integer; the maximum index in ``vertices`` up
-          to which we want to consider. If ``n = 0``, we only want to know if
-          ``v`` is adjacent to ``vertices[0]``. If ``n = 1``, we want to know
-          if ``v`` is adjacent to both ``vertices[0]`` and ``vertices[1]``.
-          Let ``k`` be the length of ``vertices``. If ``0 <= n < k``, then we
-          want to know if ``v`` is adjacent to each of
-          ``vertices[0], vertices[1], ..., vertices[n]``. Where ``n = k - 1``,
-          then we consider all elements in the list ``vertices``.
+        - ``n`` -- nonnegative integer; the maximum index in
+          ``vertices`` up to which we want to consider. If ``n = 0``,
+          we only want to know if ``(vertices[0],v)`` is an edge. If
+          ``n = 1``, we want to know whether ``(vertices[0],v)`` and
+          ``(vertices[1],v)`` are edges.  Let ``k`` be the length of
+          ``vertices``. If ``0 <= n < k``, then we want to know if
+          ``v`` is adjacent to each of ``vertices[0], vertices[1],
+          ..., vertices[n]``. Where ``n = k - 1``, then we consider
+          all elements in the list ``vertices``.
 
         - ``vertices`` -- list of vertices.
 
         - ``v`` -- a vertex.
 
+        - ``reverse`` (bool) -- if set to ``True``, considers the
+          edges ``(v,vertices[i])`` instead of ``(v,vertices[i])``
+          (only useful for digraphs).
+
         OUTPUT: 
 
-        Returns a list of ``n`` integers, whose i-th element 
-        is set to 1 iff ``v`` is adjacent to ``vertices[i]``.
+        Returns a list of ``n`` integers, whose i-th element is set to
+        `1` iff ``(v,vertices[i])`` is an edge.
+
+        .. SEEALSO::
+
+            - :meth:`adjacency_sequence_out` -- Similar method for
+            ``(v, vertices[i])`` instead of ``(vertices[i], v)`` (the
+            difference only matters for digraphs).
+        """ 
+        cdef int i = 0 
+        cdef int *seq = <int *>sage_malloc(n * sizeof(int))
+        for 0 <= i < n: 
+            seq[i] = 1 if self.has_arc_unsafe(vertices[i], v) else 0
+
+        return seq
+
+    cdef int * adjacency_sequence_out(self, int n, int *vertices, int v):
+        r""" 
+        Returns the adjacency sequence corresponding to a list of vertices
+        and a vertex. 
+
+        See the OUTPUT section for a formal definition.
+
+        See the function ``_test_adjacency_sequence()`` of
+        ``dense_graph.pyx`` and ``sparse_graph.pyx`` for unit tests.
+        
+        INPUT: 
+
+        - ``n`` -- nonnegative integer; the maximum index in
+          ``vertices`` up to which we want to consider. If ``n = 0``,
+          we only want to know if ``(v, vertices[0])`` is an edge. If
+          ``n = 1``, we want to know whether ``(v, vertices[0])`` and
+          ``(v, vertices[1])`` are edges.  Let ``k`` be the length of
+          ``vertices``. If ``0 <= n < k``, then we want to know if
+          each of ``vertices[0], vertices[1], ..., vertices[n]`` is
+          adjacent to ``v``. Where ``n = k - 1``, then we consider all
+          elements in the list ``vertices``.
+
+        - ``vertices`` -- list of vertices.
+
+        - ``v`` -- a vertex.
+
+        OUTPUT: 
+
+        Returns a list of ``n`` integers, whose i-th element is set to
+        `1` iff ``(v,vertices[i])`` is an edge.
+
+        .. SEEALSO::
+
+            - :meth:`adjacency_sequence_in` -- Similar method for
+            ``(vertices[i],v)`` instead of ``(v,vertices[i])`` (the
+            difference only matters for digraphs).
+
         """ 
         cdef int i = 0 
         cdef int *seq = <int *>sage_malloc(n * sizeof(int))
         for 0 <= i < n: 
             seq[i] = 1 if self.has_arc_unsafe(v, vertices[i]) else 0
+
         return seq
 
+
     cpdef list all_arcs(self, int u, int v):
         """
         Return the labels of all arcs from ``u`` to ``v``.

sage/graphs/base/dense_graph.pyx

@@ -625 +625 @@
         if Gnew.in_degrees[i] != Gold.in_degree(i):
             raise RuntimeError( "NO" )
 
-def _test_adjacency_sequence():
+def _test_adjacency_sequence_out():
     """
-    Randomly test the method ``DenseGraph.adjacency_sequence()``. No output
+    Randomly test the method ``DenseGraph.adjacency_sequence_out()``. No output
     indicates that no errors were found.
 
     TESTS::
 
-        sage: from sage.graphs.base.dense_graph import _test_adjacency_sequence
-        sage: _test_adjacency_sequence()  # long time
+        sage: from sage.graphs.base.dense_graph import _test_adjacency_sequence_out
+        sage: _test_adjacency_sequence_out()  # long time
     """
     from sage.graphs.digraph import DiGraph
     from sage.graphs.graph_generators import GraphGenerators
@@ -658 +658 @@
     cdef int *seq    
     for 0 <= i < randint(50, 101):
         u = randint(low, n - 1)
-        seq = g.adjacency_sequence(n, V, u)
+        seq = g.adjacency_sequence_out(n, V, u)
         A = [seq[k] for k in range(n)]
         try:
             assert A == list(M[u])

sage/graphs/base/sparse_graph.pyx

@@ -1306 +1306 @@
             if Gnew.has_arc_unsafe(i,j) != Gold.has_edge(i,j):
                 raise RuntimeError( "NO" )
 
-def _test_adjacency_sequence():
+def _test_adjacency_sequence_out():
     """
-    Randomly test the method ``SparseGraph.adjacency_sequence()``. No output
+    Randomly test the method ``SparseGraph.adjacency_sequence_out()``. No output
     indicates that no errors were found.
 
     TESTS::
 
-        sage: from sage.graphs.base.sparse_graph import _test_adjacency_sequence
-        sage: _test_adjacency_sequence()  # long time
+        sage: from sage.graphs.base.sparse_graph import _test_adjacency_sequence_out
+        sage: _test_adjacency_sequence_out()  # long time
     """
     from sage.graphs.digraph import DiGraph
     from sage.graphs.graph_generators import GraphGenerators
@@ -1341 +1341 @@
     cdef int *seq    
     for 0 <= i < randint(50, 101):
         u = randint(low, n - 1)
-        seq = g.adjacency_sequence(n, V, u)
+        seq = g.adjacency_sequence_out(n, V, u)
         A = [seq[k] for k in range(n)]
         try:
             assert A == list(M[u])

sage/graphs/generic_graph.py

@@ -8008 +8008 @@
           is no copy (induced or otherwise) of ``G`` in ``self``, we return
           ``None``.
 
+        .. NOTE::
+
+            This method also works on digraphs.
+
         ALGORITHM:
 
         Brute-force search.
@@ -8058 +8062 @@
 
         The empty graph is a subgraph of every graph::
 
-            sage: g.subgraph_search(graphs.EmptyGraph())
-            Graph on 0 vertices
-            sage: g.subgraph_search(graphs.EmptyGraph(), induced=True)
-            Graph on 0 vertices
-        """
-        from sage.graphs.generic_graph_pyx import subgraph_search
+             sage: g.subgraph_search(graphs.EmptyGraph())
+             Graph on 0 vertices
+             sage: g.subgraph_search(graphs.EmptyGraph(), induced=True)
+             Graph on 0 vertices
+        """
+        from sage.graphs.generic_graph_pyx import SubgraphSearch
         from sage.graphs.graph_generators import GraphGenerators
         if G.order() == 0:
             return GraphGenerators().EmptyGraph()
-        H = subgraph_search(self, G, induced=induced)
-        if H == []:
-            return None
-        return self.subgraph(H)
-    
+
+        S = SubgraphSearch(self, G, induced = induced)
+
+        for g in S:
+            return self.subgraph(g)
+
+        return None
+
+    def subgraph_search_count(self, G, induced=False):
+        r"""
+        Returns the number of labelled occurences of ``G`` in ``self``.
+
+        INPUT:
+
+        - ``G`` -- the graph whose copies we are looking for in
+          ``self``.
+
+        - ``induced`` -- boolean (default: ``False``). Whether or not
+          to count induced copies of ``G`` in ``self``.
+
+        ALGORITHM:
+
+        Brute-force search.
+
+        .. NOTE::
+
+            This method also works on digraphs.
+
+        EXAMPLES:
+
+        Counting the number of paths `P_5` in a PetersenGraph::
+
+            sage: g = graphs.PetersenGraph()
+            sage: g.subgraph_search_count(graphs.PathGraph(5))
+            240
+
+        Requiring these subgraphs be induced::
+
+            sage: g.subgraph_search_count(graphs.PathGraph(5), induced = True)
+            120
+
+        If we define the graph `T_k` (the transitive tournament on `k`
+        vertices) as the graph on `\{0, ..., k-1\}` such that `ij \in
+        T_k` iif `i<j`, how many directed triangles can be found in
+        `T_5` ? The answer is of course `0` ::
+
+             sage: T5 = DiGraph()
+             sage: T5.add_edges([(i,j) for i in xrange(5) for j in xrange(i+1, 5)])
+             sage: T5.subgraph_search_count(digraphs.Circuit(3))
+             0
+
+        If we count instead the number of `T_3` in `T_5`, we expect
+        the answer to be `{5 \choose 3}`::
+
+             sage: T3 = T5.subgraph([0,1,2])
+             sage: T5.subgraph_search_count(T3)
+             10
+             sage: binomial(5,3)
+             10
+
+        The empty graph is a subgraph of every graph::
+
+            sage: g.subgraph_search_count(graphs.EmptyGraph())
+            1
+        """
+        from sage.graphs.generic_graph_pyx import SubgraphSearch
+
+        if G.order() == 0:
+            return 1
+
+        if self.order() == 0:
+            return 0
+
+        S = SubgraphSearch(self, G, induced = induced)
+
+        return S.cardinality()
+
+    def subgraph_search_iterator(self, G, induced=False):
+        r"""
+        Returns an iterator over the labelled copies of ``G`` in ``self``.
+
+        INPUT:
+
+        - ``G`` -- the graph whose copies we are looking for in
+          ``self``.
+
+        - ``induced`` -- boolean (default: ``False``). Whether or not
+          to iterate over the induced copies of ``G`` in ``self``.
+
+        ALGORITHM:
+
+        Brute-force search.
+
+        OUTPUT:
+
+            Iterator over the labelled copies of ``G`` in ``self``, as
+            *lists*. For each value `(v_1, v_2, ..., v_k)` returned,
+            the first vertex of `G` is associated with `v_1`, the
+            second with `v_2`, etc ...
+
+        .. NOTE::
+
+            This method also works on digraphs.
+
+        EXAMPLE:
+
+        Iterating through all the labelled `P_3` of `P_5`::
+        
+            sage: g = graphs.PathGraph(5)
+            sage: for p in g.subgraph_search_iterator(graphs.PathGraph(3)):
+            ...      print p
+            [0, 1, 2]
+            [1, 2, 3]
+            [2, 1, 0]
+            [2, 3, 4]
+            [3, 2, 1]
+            [4, 3, 2]
+        """
+
+        if G.order() == 0:
+            from sage.graphs.graph_generators import GraphGenerators
+            return [GraphGenerators().EmptyGraph()]
+
+        elif self.order() == 0:
+            return []
+
+        else:
+            from sage.graphs.generic_graph_pyx import SubgraphSearch
+            return SubgraphSearch(self, G, induced = induced)
+
     def random_subgraph(self, p, inplace=False):
         """
         Return a random subgraph that contains each vertex with prob. p.

sage/graphs/generic_graph_pyx.pxd

@@ -1 +1 @@
 
 from sage.structure.sage_object cimport SageObject
+from sage.graphs.base.dense_graph cimport DenseGraph
 
 cdef run_spring(int, int, double*, int*, int, bint)
 
 cdef class GenericGraph_pyx(SageObject):
     pass
 
+
+
+cdef class SubgraphSearch:
+    cdef int ng
+    cdef int nh
+    cdef (bint) (*is_admissible) (int, int *, int *)
+    cdef DenseGraph g
+    cdef DenseGraph h
+    cdef int *busy
+    cdef int *stack
+    cdef int active
+    cdef int *vertices
+    cdef int **line_h_out
+    cdef int **line_h_in
+    cdef list g_vertices
+    cdef int i
+    cdef bool directed
+
+
 cdef inline bint vectors_equal(int n, int *a, int *b)
 cdef inline bint vectors_inferior(int n, int *a, int *b)

sage/graphs/generic_graph_pyx.pyx

@@ -22 +22 @@
 # import from Python standard library
 from random import random
 
-# import from third-party library
-from sage.graphs.base.dense_graph cimport DenseGraph
 
 cdef extern from *:
     double sqrt(double)
@@ -458 +456 @@
 
 # Exhaustive search in graphs
 
-cpdef list subgraph_search(G, H, bint induced=False):
-    r"""
-    Returns a set of vertices in ``G`` representing a copy of ``H``.
+cdef class SubgraphSearch:
+    r""" 
+    This class implements methods to exhaustively search for labelled
+    copies of a graph `H` in a larger graph `G`.
+
+    It is possible to look for induced subgraphs instead, and to
+    iterate or count the number of their occurrences.
 
     ALGORITHM:
 
-    This algorithm is a brute-force search. 
-    Let `V(H) = \{h_1,\dots,h_k\}`.  It first tries
-    to find in `G` a possible representant of `h_1`, then a 
-    representant of `h_2` compatible with `h_1`, then
-    a representant of `h_3` compatible with the first
+    The algorithm is a brute-force search.  Let `V(H) =
+    \{h_1,\dots,h_k\}`.  It first tries to find in `G` a possible
+    representant of `h_1`, then a representant of `h_2` compatible
+    with `h_1`, then a representant of `h_3` compatible with the first
     two, etc.
 
-    This way, most of the time we need to test far less than
-    `k! \binom{|V(G)|}{k}` subsets, and hope this brute-force
-    technique can sometimes be useful.
+    This way, most of the time we need to test far less than `k!
+    \binom{|V(G)|}{k}` subsets, and hope this brute-force technique
+    can sometimes be useful.
+    """
+    def __init__(self, G, H, induced = False):
+        r"""
+        Constructor
 
-    INPUT:
+        This constructor only checks there is no inconsistency in the
+        input : `G` and `H` are both graphs or both digraphs.
+        """
+        if sum([G.is_directed(), H.is_directed()]) == 1:
+            raise ValueError("One graph can not be directed while the other is not.")
 
-    - ``G``, ``H`` -- two graphs such that ``H`` is a subgraph of ``G``.
+    def __iter__(self):
+        r""" 
+        Returns an iterator over all the labeleld subgraphs of `G`
+        isomorphic to `H`.
 
-    - ``induced`` -- boolean (default: ``False``); whether to require that
-      the subgraph is an induced subgraph.
+        EXAMPLE:
 
-    OUTPUT:
+        Iterating through all the `P_3` of `P_5`::
+        
+            sage: from sage.graphs.generic_graph_pyx import SubgraphSearch
+            sage: g = graphs.PathGraph(5)
+            sage: h = graphs.PathGraph(3)
+            sage: S = SubgraphSearch(g, h)
+            sage: for p in S:
+            ...      print p
+            [0, 1, 2]
+            [1, 2, 3]
+            [2, 1, 0]
+            [2, 3, 4]
+            [3, 2, 1]
+            [4, 3, 2]
+        """
+        self._initialization()
+        return self
 
-    A list of vertices inducing a copy of ``H`` in ``G``. If none is found,
-    an empty list is returned.
+    def cardinality(self):
+        r"""
+        Returns the number of labelled subgraphs of `G` isomorphic to
+        `H`.
 
-    EXAMPLES:
+        .. NOTE::
 
-    A Petersen graph contains an induced path graph `P_5`::
+           This method counts the subgraphs by enumerating them all !
+           Hence it probably is not a good idea to count their number
+           before enumerating them :-)
 
-        sage: from sage.graphs.generic_graph_pyx import subgraph_search
-        sage: g = graphs.PetersenGraph()
-        sage: subgraph_search(g, graphs.PathGraph(5), induced=True)
-        [0, 1, 2, 3, 8]
+        EXAMPLE:
 
-    It also contains a the claw `K_{1,3}`::
+        Counting the number of labelled `P_3` in `P_5`::
+        
+            sage: from sage.graphs.generic_graph_pyx import SubgraphSearch
+            sage: g = graphs.PathGraph(5)
+            sage: h = graphs.PathGraph(3)
+            sage: S = SubgraphSearch(g, h)
+            sage: S.cardinality()
+            6
+        """
+        if self.nh > self.ng:
+            return 0
 
-        sage: subgraph_search(g, graphs.ClawGraph())
-        [0, 1, 4, 5]
+        self._initialization()
+        cdef int i
 
-    Though it contains no induced `P_6`::
+        i=0
+        for _ in self:
+            i+=1
 
-        sage: subgraph_search(g, graphs.PathGraph(6), induced=True)
-        []
+        return i
 
-    TESTS:
+    def _initialization(self):
+        r"""
+        Initialization of the variables.
 
-    Let `G` and `H` be graphs having orders `m` and `n`, respectively. If
-    `m < n`, then there are no copies of `H` in `G`::
+        Once the memory allocation is done in :meth:`__cinit__`,
+        several variables need to be set to a default value. As this
+        operation needs to be performed before any call to
+        :meth:`__iter__` or to :meth:`cardinality`, it is cleaner to
+        create a dedicated method.
 
-        sage: from sage.graphs.generic_graph_pyx import subgraph_search
-        sage: m = randint(100, 200)
-        sage: n = randint(m + 1, 300)
-        sage: G = graphs.RandomGNP(m, random())
-        sage: H = graphs.RandomGNP(n, random())
-        sage: G.order() < H.order()
-        True
-        sage: subgraph_search(G, H)
-        []
-    """
-    # TODO: This is a brute-force search and can be very inefficient. Write
-    # a more efficient subgraph search implementation.
-    cdef int ng = G.order()
-    cdef int nh = H.order()
-    if ng < nh:
-        return []
-    cdef int i, j, k
-    cdef int *tmp_array
-    cdef (bint) (*is_admissible) (int, int *, int *)
-    if induced:
-        is_admissible = vectors_equal
-    else:
-        is_admissible = vectors_inferior
-    # static copies of the two graphs for more efficient operations
-    cdef DenseGraph g = DenseGraph(ng)
-    cdef DenseGraph h = DenseGraph(nh)
-    # copying the adjacency relations in both G and H
-    i = 0
-    for row in G.adjacency_matrix():
-        j = 0
-        for k in row:
-            if k:
-                g.add_arc(i, j)
-            j += 1
-        i += 1
-    i = 0
-    for row in H.adjacency_matrix():
-        j = 0
-        for k in row:
-            if k:
-                h.add_arc(i, j)
-            j += 1
-        i += 1
-    # A vertex is said to be busy if it is already part of the partial copy
-    # of H in G.
-    cdef int *busy = <int *>sage_malloc(ng * sizeof(int))
-    memset(busy, 0, ng * sizeof(int))
-    # 0 is the first vertex we use, so it is at first busy
-    busy[0] = 1
-    # stack -- list of the vertices which are part of the partial copy of H
-    # in G.
-    #
-    # stack[i] -- the integer corresponding to the vertex of G representing
-    # the i-th vertex of H.
-    #
-    # stack[i] = -1 means that i is not represented
-    # ... yet!
-    cdef int *stack = <int *>sage_malloc(nh * sizeof(int))
-    stack[0] = 0
-    stack[1] = -1
-    # Number of representants we have already found. Set to 1 as vertex 0
-    # is already part of the partial copy of H in G.
-    cdef int active = 1
-    # vertices is equal to range(nh), as an int *variable
-    cdef int *vertices = <int *>sage_malloc(nh * sizeof(int))
-    for 0 <= i < nh:
-        vertices[i] = i
-    # line_h[i] represents the adjacency sequence of vertex i
-    # in h relative to vertices 0, 1, ..., i-1
-    cdef int **line_h = <int **>sage_malloc(nh * sizeof(int *))
-    for 0 <= i < nh:
-        line_h[i] = <int *>h.adjacency_sequence(i, vertices, i)
-    # the sequence of vertices to be returned
-    cdef list value = []
-    
-    _sig_on
+        EXAMPLE:
 
-    # as long as there is a non-void partial copy of H in G
-    while active:
-        # If we are here and found nothing yet, we try the next possible
-        # vertex as a representant of the active i-th vertex of H.
-        i = stack[active] + 1
-        # Looking for a vertex that is not busy and compatible with the
-        # partial copy we have of H.
-        while i < ng:
-            if busy[i]:
-                i += 1
+        Finding two times the first occurrence through the
+        re-initialization of the instance ::
+
+            sage: from sage.graphs.generic_graph_pyx import SubgraphSearch
+            sage: g = graphs.PathGraph(5)
+            sage: h = graphs.PathGraph(3)
+            sage: S = SubgraphSearch(g, h)
+            sage: S.__next__()
+            [0, 1, 2]
+            sage: S._initialization()
+            sage: S.__next__()
+            [0, 1, 2]
+        """
+
+        memset(self.busy, 0, self.ng * sizeof(int))
+        # 0 is the first vertex we use, so it is at first busy
+        self.busy[0] = 1
+        # stack -- list of the vertices which are part of the partial copy of H
+        # in G.
+        #
+        # stack[i] -- the integer corresponding to the vertex of G representing
+        # the i-th vertex of H.
+        #
+        # stack[i] = -1 means that i is not represented
+        # ... yet!
+
+        self.stack[0] = 0
+        self.stack[1] = -1
+
+        # Number of representants we have already found. Set to 1 as vertex 0
+        # is already part of the partial copy of H in G.
+        self.active = 1
+
+    def __cinit__(self, G, H, induced = False):
+        r"""
+        Cython constructor
+
+        This method initializes all the C values.
+        """
+
+        # Storing the number of vertices
+        self.ng = G.order()
+        self.nh = H.order()
+
+
+
+        # Storing the list of vertices
+        self.g_vertices = G.vertices()
+
+        # Are the graphs directed (at the end of the code is checked
+        # whether both are of the same type)
+        self.directed = G.is_directed()
+
+        cdef int i, j, k
+        cdef int *tmp_array
+
+        # Should we look for induced subgraphs ?
+        if induced:
+            self.is_admissible = vectors_equal
+        else:
+            self.is_admissible = vectors_inferior
+
+        # static copies of the two graphs for more efficient operations
+        self.g = DenseGraph(self.ng)
+        self.h = DenseGraph(self.nh)
+
+        # copying the adjacency relations in both G and H
+        i = 0
+        for row in G.adjacency_matrix():
+            j = 0
+            for k in row:
+                if k:
+                    self.g.add_arc(i, j)
+                j += 1
+            i += 1
+        i = 0
+        for row in H.adjacency_matrix():
+            j = 0
+            for k in row:
+                if k:
+                    self.h.add_arc(i, j)
+                j += 1
+            i += 1
+
+        # A vertex is said to be busy if it is already part of the partial copy
+        # of H in G.
+        self.busy = <int *>sage_malloc(self.ng * sizeof(int))
+        self.stack = <int *>sage_malloc(self.nh * sizeof(int))
+
+        # vertices is equal to range(nh), as an int *variable
+        self.vertices = <int *>sage_malloc(self.nh * sizeof(int))
+        for 0 <= i < self.nh:
+            self.vertices[i] = i
+
+        # line_h_out[i] represents the adjacency sequence of vertex i
+        # in h relative to vertices 0, 1, ..., i-1
+        self.line_h_out = <int **>sage_malloc(self.nh * sizeof(int *))
+        for 0 <= i < self.nh:
+            self.line_h_out[i] = <int *>self.h.adjacency_sequence_out(i, self.vertices, i)
+
+        # Similarly in the opposite direction (only useful if the
+        # graphs are directed)
+        if self.directed:
+            self.line_h_in = <int **>sage_malloc(self.nh * sizeof(int *))
+            for 0 <= i < self.nh:
+                self.line_h_in[i] = <int *>self.h.adjacency_sequence_in(i, self.vertices, i)
+
+        self._initialization()
+
+    def __next__(self):
+        r"""
+        Returns the next isomorphic subgraph if any, and raises a
+        ``StopIteration`` otherwise.
+
+        EXAMPLE::
+
+            sage: from sage.graphs.generic_graph_pyx import SubgraphSearch
+            sage: g = graphs.PathGraph(5)
+            sage: h = graphs.PathGraph(3)
+            sage: S = SubgraphSearch(g, h)
+            sage: S.__next__()
+            [0, 1, 2]
+        """
+        _sig_on
+        cdef int *tmp_array_out
+        cdef int *tmp_array_in
+        cdef bool is_admissible
+
+        # as long as there is a non-void partial copy of H in G
+        while self.active >= 0:
+            # If we are here and found nothing yet, we try the next possible
+            # vertex as a representant of the active i-th vertex of H.
+            self.i = self.stack[self.active] + 1
+            # Looking for a vertex that is not busy and compatible with the
+            # partial copy we have of H.
+            while self.i < self.ng:
+                if self.busy[self.i]:
+                    self.i += 1
+                else:
+                    # Testing whether the vertex we picked is a
+                    # correct extension by checking the edges from the
+                    # vertices already selected to self.i satisfy the
+                    # constraints
+                    tmp_array_out = self.g.adjacency_sequence_out(self.active, self.stack, self.i)
+                    is_admissible = self.is_admissible(self.active, tmp_array_out, self.line_h_out[self.active])
+                    sage_free(tmp_array_out)
+
+                    # If G and H are digraphs, we also need to ensure
+                    # the edges going in the opposite direction
+                    # satisfy the constraints
+                    if is_admissible and self.directed:
+                        tmp_array_in = self.g.adjacency_sequence_in(self.active, self.stack, self.i)
+                        is_admissible = is_admissible and self.is_admissible(self.active, tmp_array_in, self.line_h_in[self.active])
+                        sage_free(tmp_array_in)
+
+                    if is_admissible:
+                        break
+                    else:
+                        self.i += 1                        
+
+            # If we have found a good representant of H's i-th vertex in G
+            if self.i < self.ng:
+
+                # updating the last vertex of the stack
+                if self.stack[self.active] != -1:
+                    self.busy[self.stack[self.active]] = 0
+                self.stack[self.active] = self.i
+
+                # We have found our copy !!!
+                if self.active == self.nh-1:
+                    return [self.g_vertices[self.stack[l]] for l in xrange(self.nh)]
+
+                # We are still missing several vertices ...
+                else:
+                    self.busy[self.stack[self.active]] = 1
+                    self.active += 1
+
+                    # we begin the search of the next vertex at 0
+                    self.stack[self.active] = -1
+
+            # If we found no representant for the i-th vertex, it
+            # means that we cannot extend the current copy of H so we
+            # update the status of stack[active] and prepare to change
+            # the previous vertex.
+
             else:
-                tmp_array = g.adjacency_sequence(active, stack, i)
-                if is_admissible(active, tmp_array, line_h[active]):
-                    sage_free(tmp_array)
-                    break
-                else:
-                    sage_free(tmp_array)
-                    i += 1
-        # If we found none, it means that we cannot extend the current copy
-        # of H so we update the status of stack[active] and prepare to change
-        # the previous vertex.
-        if i >= ng:
-            if stack[active] != -1:
-                busy[stack[active]] = 0
-            stack[active] = -1
-            active -= 1
-        # If we have found a good representant of H's i-th vertex in G
-        else:
-            if stack[active] != -1:
-                busy[stack[active]] = 0
-            stack[active] = i
-            busy[stack[active]] = 1
-            active += 1
-            # We have found our copy!!!
-            if active == nh:
-                g_vertices = G.vertices()
-                value = [g_vertices[stack[i]] for i in xrange(nh)]
-                break
-            else:
-                # we begin the search of the next vertex at 0
-                stack[active] = -1
+                if self.stack[self.active] != -1:
+                    self.busy[self.stack[self.active]] = 0
+                self.stack[self.active] = -1
+                self.active -= 1
 
-    _sig_off
+        _sig_off
+        raise StopIteration
 
-    # Free the memory
-    sage_free(busy)
-    sage_free(stack)
-    sage_free(vertices)
-    for 0 <= i < nh:
-        sage_free(line_h[i])
-    sage_free(line_h)
+    def __dealloc__(self):
+        r"""
+        Freeing the allocated memory.
+        """
 
-    return value
+        # Free the memory
+        sage_free(self.busy)
+        sage_free(self.stack)
+        sage_free(self.vertices)
+        for 0 <= i < self.nh:
+            sage_free(self.line_h_out[i])
+        sage_free(self.line_h_out)
+
+        if self.directed:
+            for 0 <= i < self.nh:
+                sage_free(self.line_h_in[i])
+            sage_free(self.line_h_in)
 
 cdef inline bint vectors_equal(int n, int *a, int *b):
     r"""