Ticket #9420: trac_9420.patch

File trac_9420.patch, 31.3 KB (added by ncohen, 11 years ago)
  • sage/graphs/base/c_graph.pxd

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1278157947 -7200
    # Node ID 20b8ec1a5d34e2633a8d6ce057f4f19b0165fd86
    # Parent  113df1e3b94772d3d5cf45f213f48a0c480a25ec
    trac 9420 -- SubgraphSearch class
    
    diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/base/c_graph.pxd
    a b  
    2828    cpdef add_arc(self, int u, int v)
    2929    cpdef bint has_arc(self, int u, int v)
    3030    cpdef del_all_arcs(self, int u, int v)
    31     cdef int *adjacency_sequence(self, int n, int *vertices, int v)
     31    cdef int *adjacency_sequence_in(self, int n, int *vertices, int v)
     32    cdef int *adjacency_sequence_out(self, int n, int *vertices, int v)
    3233    cpdef list all_arcs(self, int u, int v)
    3334    cpdef list in_neighbors(self, int v)
    3435    cpdef list out_neighbors(self, int u)
  • sage/graphs/base/c_graph.pyx

    diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/base/c_graph.pyx
    a b  
    840840        """
    841841        raise NotImplementedError()
    842842
    843     cdef int *adjacency_sequence(self, int n, int *vertices, int v):
     843    cdef int * adjacency_sequence_in(self, int n, int *vertices, int v):
    844844        r"""
    845845        Returns the adjacency sequence corresponding to a list of vertices
    846         and a vertex. See the function ``_test_adjacency_sequence()`` of
     846        and a vertex.
     847
     848        See the OUTPUT section for a formal definition.
     849
     850        See the function ``_test_adjacency_sequence()`` of
    847851        ``dense_graph.pyx`` and ``sparse_graph.pyx`` for unit tests.
    848852       
    849853        INPUT:
    850854
    851         - ``n`` -- nonnegative integer; the maximum index in ``vertices`` up
    852           to which we want to consider. If ``n = 0``, we only want to know if
    853           ``v`` is adjacent to ``vertices[0]``. If ``n = 1``, we want to know
    854           if ``v`` is adjacent to both ``vertices[0]`` and ``vertices[1]``.
    855           Let ``k`` be the length of ``vertices``. If ``0 <= n < k``, then we
    856           want to know if ``v`` is adjacent to each of
    857           ``vertices[0], vertices[1], ..., vertices[n]``. Where ``n = k - 1``,
    858           then we consider all elements in the list ``vertices``.
     855        - ``n`` -- nonnegative integer; the maximum index in
     856          ``vertices`` up to which we want to consider. If ``n = 0``,
     857          we only want to know if ``(vertices[0],v)`` is an edge. If
     858          ``n = 1``, we want to know whether ``(vertices[0],v)`` and
     859          ``(vertices[1],v)`` are edges.  Let ``k`` be the length of
     860          ``vertices``. If ``0 <= n < k``, then we want to know if
     861          ``v`` is adjacent to each of ``vertices[0], vertices[1],
     862          ..., vertices[n]``. Where ``n = k - 1``, then we consider
     863          all elements in the list ``vertices``.
    859864
    860865        - ``vertices`` -- list of vertices.
    861866
    862867        - ``v`` -- a vertex.
    863868
     869        - ``reverse`` (bool) -- if set to ``True``, considers the
     870          edges ``(v,vertices[i])`` instead of ``(v,vertices[i])``
     871          (only useful for digraphs).
     872
    864873        OUTPUT:
    865874
    866         Returns a list of ``n`` integers, whose i-th element
    867         is set to 1 iff ``v`` is adjacent to ``vertices[i]``.
     875        Returns a list of ``n`` integers, whose i-th element is set to
     876        `1` iff ``(v,vertices[i])`` is an edge.
     877
     878        .. SEEALSO::
     879
     880            - :meth:`adjacency_sequence_out` -- Similar method for
     881            ``(v, vertices[i])`` instead of ``(vertices[i], v)`` (the
     882            difference only matters for digraphs).
     883        """
     884        cdef int i = 0
     885        cdef int *seq = <int *>sage_malloc(n * sizeof(int))
     886        for 0 <= i < n:
     887            seq[i] = 1 if self.has_arc_unsafe(vertices[i], v) else 0
     888
     889        return seq
     890
     891    cdef int * adjacency_sequence_out(self, int n, int *vertices, int v):
     892        r"""
     893        Returns the adjacency sequence corresponding to a list of vertices
     894        and a vertex.
     895
     896        See the OUTPUT section for a formal definition.
     897
     898        See the function ``_test_adjacency_sequence()`` of
     899        ``dense_graph.pyx`` and ``sparse_graph.pyx`` for unit tests.
     900       
     901        INPUT:
     902
     903        - ``n`` -- nonnegative integer; the maximum index in
     904          ``vertices`` up to which we want to consider. If ``n = 0``,
     905          we only want to know if ``(v, vertices[0])`` is an edge. If
     906          ``n = 1``, we want to know whether ``(v, vertices[0])`` and
     907          ``(v, vertices[1])`` are edges.  Let ``k`` be the length of
     908          ``vertices``. If ``0 <= n < k``, then we want to know if
     909          each of ``vertices[0], vertices[1], ..., vertices[n]`` is
     910          adjacent to ``v``. Where ``n = k - 1``, then we consider all
     911          elements in the list ``vertices``.
     912
     913        - ``vertices`` -- list of vertices.
     914
     915        - ``v`` -- a vertex.
     916
     917        OUTPUT:
     918
     919        Returns a list of ``n`` integers, whose i-th element is set to
     920        `1` iff ``(v,vertices[i])`` is an edge.
     921
     922        .. SEEALSO::
     923
     924            - :meth:`adjacency_sequence_in` -- Similar method for
     925            ``(vertices[i],v)`` instead of ``(v,vertices[i])`` (the
     926            difference only matters for digraphs).
     927
    868928        """
    869929        cdef int i = 0
    870930        cdef int *seq = <int *>sage_malloc(n * sizeof(int))
    871931        for 0 <= i < n:
    872932            seq[i] = 1 if self.has_arc_unsafe(v, vertices[i]) else 0
     933
    873934        return seq
    874935
     936
    875937    cpdef list all_arcs(self, int u, int v):
    876938        """
    877939        Return the labels of all arcs from ``u`` to ``v``.
  • sage/graphs/base/dense_graph.pyx

    diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/base/dense_graph.pyx
    a b  
    625625        if Gnew.in_degrees[i] != Gold.in_degree(i):
    626626            raise RuntimeError( "NO" )
    627627
    628 def _test_adjacency_sequence():
     628def _test_adjacency_sequence_out():
    629629    """
    630     Randomly test the method ``DenseGraph.adjacency_sequence()``. No output
     630    Randomly test the method ``DenseGraph.adjacency_sequence_out()``. No output
    631631    indicates that no errors were found.
    632632
    633633    TESTS::
    634634
    635         sage: from sage.graphs.base.dense_graph import _test_adjacency_sequence
    636         sage: _test_adjacency_sequence()  # long time
     635        sage: from sage.graphs.base.dense_graph import _test_adjacency_sequence_out
     636        sage: _test_adjacency_sequence_out()  # long time
    637637    """
    638638    from sage.graphs.digraph import DiGraph
    639639    from sage.graphs.graph_generators import GraphGenerators
     
    658658    cdef int *seq   
    659659    for 0 <= i < randint(50, 101):
    660660        u = randint(low, n - 1)
    661         seq = g.adjacency_sequence(n, V, u)
     661        seq = g.adjacency_sequence_out(n, V, u)
    662662        A = [seq[k] for k in range(n)]
    663663        try:
    664664            assert A == list(M[u])
  • sage/graphs/base/sparse_graph.pyx

    diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/base/sparse_graph.pyx
    a b  
    13061306            if Gnew.has_arc_unsafe(i,j) != Gold.has_edge(i,j):
    13071307                raise RuntimeError( "NO" )
    13081308
    1309 def _test_adjacency_sequence():
     1309def _test_adjacency_sequence_out():
    13101310    """
    1311     Randomly test the method ``SparseGraph.adjacency_sequence()``. No output
     1311    Randomly test the method ``SparseGraph.adjacency_sequence_out()``. No output
    13121312    indicates that no errors were found.
    13131313
    13141314    TESTS::
    13151315
    1316         sage: from sage.graphs.base.sparse_graph import _test_adjacency_sequence
    1317         sage: _test_adjacency_sequence()  # long time
     1316        sage: from sage.graphs.base.sparse_graph import _test_adjacency_sequence_out
     1317        sage: _test_adjacency_sequence_out()  # long time
    13181318    """
    13191319    from sage.graphs.digraph import DiGraph
    13201320    from sage.graphs.graph_generators import GraphGenerators
     
    13411341    cdef int *seq   
    13421342    for 0 <= i < randint(50, 101):
    13431343        u = randint(low, n - 1)
    1344         seq = g.adjacency_sequence(n, V, u)
     1344        seq = g.adjacency_sequence_out(n, V, u)
    13451345        A = [seq[k] for k in range(n)]
    13461346        try:
    13471347            assert A == list(M[u])
  • sage/graphs/generic_graph.py

    diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/generic_graph.py
    a b  
    80088008          is no copy (induced or otherwise) of ``G`` in ``self``, we return
    80098009          ``None``.
    80108010
     8011        .. NOTE::
     8012
     8013            This method also works on digraphs.
     8014
    80118015        ALGORITHM:
    80128016
    80138017        Brute-force search.
     
    80588062
    80598063        The empty graph is a subgraph of every graph::
    80608064
    8061             sage: g.subgraph_search(graphs.EmptyGraph())
    8062             Graph on 0 vertices
    8063             sage: g.subgraph_search(graphs.EmptyGraph(), induced=True)
    8064             Graph on 0 vertices
    8065         """
    8066         from sage.graphs.generic_graph_pyx import subgraph_search
     8065             sage: g.subgraph_search(graphs.EmptyGraph())
     8066             Graph on 0 vertices
     8067             sage: g.subgraph_search(graphs.EmptyGraph(), induced=True)
     8068             Graph on 0 vertices
     8069        """
     8070        from sage.graphs.generic_graph_pyx import SubgraphSearch
    80678071        from sage.graphs.graph_generators import GraphGenerators
    80688072        if G.order() == 0:
    80698073            return GraphGenerators().EmptyGraph()
    8070         H = subgraph_search(self, G, induced=induced)
    8071         if H == []:
    8072             return None
    8073         return self.subgraph(H)
    8074    
     8074
     8075        S = SubgraphSearch(self, G, induced = induced)
     8076
     8077        for g in S:
     8078            return self.subgraph(g)
     8079
     8080        return None
     8081
     8082    def subgraph_search_count(self, G, induced=False):
     8083        r"""
     8084        Returns the number of labelled occurences of ``G`` in ``self``.
     8085
     8086        INPUT:
     8087
     8088        - ``G`` -- the graph whose copies we are looking for in
     8089          ``self``.
     8090
     8091        - ``induced`` -- boolean (default: ``False``). Whether or not
     8092          to count induced copies of ``G`` in ``self``.
     8093
     8094        ALGORITHM:
     8095
     8096        Brute-force search.
     8097
     8098        .. NOTE::
     8099
     8100            This method also works on digraphs.
     8101
     8102        EXAMPLES:
     8103
     8104        Counting the number of paths `P_5` in a PetersenGraph::
     8105
     8106            sage: g = graphs.PetersenGraph()
     8107            sage: g.subgraph_search_count(graphs.PathGraph(5))
     8108            240
     8109
     8110        Requiring these subgraphs be induced::
     8111
     8112            sage: g.subgraph_search_count(graphs.PathGraph(5), induced = True)
     8113            120
     8114
     8115        If we define the graph `T_k` (the transitive tournament on `k`
     8116        vertices) as the graph on `\{0, ..., k-1\}` such that `ij \in
     8117        T_k` iif `i<j`, how many directed triangles can be found in
     8118        `T_5` ? The answer is of course `0` ::
     8119
     8120             sage: T5 = DiGraph()
     8121             sage: T5.add_edges([(i,j) for i in xrange(5) for j in xrange(i+1, 5)])
     8122             sage: T5.subgraph_search_count(digraphs.Circuit(3))
     8123             0
     8124
     8125        If we count instead the number of `T_3` in `T_5`, we expect
     8126        the answer to be `{5 \choose 3}`::
     8127
     8128             sage: T3 = T5.subgraph([0,1,2])
     8129             sage: T5.subgraph_search_count(T3)
     8130             10
     8131             sage: binomial(5,3)
     8132             10
     8133
     8134        The empty graph is a subgraph of every graph::
     8135
     8136            sage: g.subgraph_search_count(graphs.EmptyGraph())
     8137            1
     8138        """
     8139        from sage.graphs.generic_graph_pyx import SubgraphSearch
     8140
     8141        if G.order() == 0:
     8142            return 1
     8143
     8144        if self.order() == 0:
     8145            return 0
     8146
     8147        S = SubgraphSearch(self, G, induced = induced)
     8148
     8149        return S.cardinality()
     8150
     8151    def subgraph_search_iterator(self, G, induced=False):
     8152        r"""
     8153        Returns an iterator over the labelled copies of ``G`` in ``self``.
     8154
     8155        INPUT:
     8156
     8157        - ``G`` -- the graph whose copies we are looking for in
     8158          ``self``.
     8159
     8160        - ``induced`` -- boolean (default: ``False``). Whether or not
     8161          to iterate over the induced copies of ``G`` in ``self``.
     8162
     8163        ALGORITHM:
     8164
     8165        Brute-force search.
     8166
     8167        OUTPUT:
     8168
     8169            Iterator over the labelled copies of ``G`` in ``self``, as
     8170            *lists*. For each value `(v_1, v_2, ..., v_k)` returned,
     8171            the first vertex of `G` is associated with `v_1`, the
     8172            second with `v_2`, etc ...
     8173
     8174        .. NOTE::
     8175
     8176            This method also works on digraphs.
     8177
     8178        EXAMPLE:
     8179
     8180        Iterating through all the labelled `P_3` of `P_5`::
     8181       
     8182            sage: g = graphs.PathGraph(5)
     8183            sage: for p in g.subgraph_search_iterator(graphs.PathGraph(3)):
     8184            ...      print p
     8185            [0, 1, 2]
     8186            [1, 2, 3]
     8187            [2, 1, 0]
     8188            [2, 3, 4]
     8189            [3, 2, 1]
     8190            [4, 3, 2]
     8191        """
     8192
     8193        if G.order() == 0:
     8194            from sage.graphs.graph_generators import GraphGenerators
     8195            return [GraphGenerators().EmptyGraph()]
     8196
     8197        elif self.order() == 0:
     8198            return []
     8199
     8200        else:
     8201            from sage.graphs.generic_graph_pyx import SubgraphSearch
     8202            return SubgraphSearch(self, G, induced = induced)
     8203
    80758204    def random_subgraph(self, p, inplace=False):
    80768205        """
    80778206        Return a random subgraph that contains each vertex with prob. p.
  • sage/graphs/generic_graph_pyx.pxd

    diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/generic_graph_pyx.pxd
    a b  
    11
    22from sage.structure.sage_object cimport SageObject
     3from sage.graphs.base.dense_graph cimport DenseGraph
    34
    45cdef run_spring(int, int, double*, int*, int, bint)
    56
    67cdef class GenericGraph_pyx(SageObject):
    78    pass
    89
     10
     11
     12cdef class SubgraphSearch:
     13    cdef int ng
     14    cdef int nh
     15    cdef (bint) (*is_admissible) (int, int *, int *)
     16    cdef DenseGraph g
     17    cdef DenseGraph h
     18    cdef int *busy
     19    cdef int *stack
     20    cdef int active
     21    cdef int *vertices
     22    cdef int **line_h_out
     23    cdef int **line_h_in
     24    cdef list g_vertices
     25    cdef int i
     26    cdef bool directed
     27
     28
    929cdef inline bint vectors_equal(int n, int *a, int *b)
    1030cdef inline bint vectors_inferior(int n, int *a, int *b)
  • sage/graphs/generic_graph_pyx.pyx

    diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/generic_graph_pyx.pyx
    a b  
    2222# import from Python standard library
    2323from random import random
    2424
    25 # import from third-party library
    26 from sage.graphs.base.dense_graph cimport DenseGraph
    2725
    2826cdef extern from *:
    2927    double sqrt(double)
     
    458456
    459457# Exhaustive search in graphs
    460458
    461 cpdef list subgraph_search(G, H, bint induced=False):
    462     r"""
    463     Returns a set of vertices in ``G`` representing a copy of ``H``.
     459cdef class SubgraphSearch:
     460    r"""
     461    This class implements methods to exhaustively search for labelled
     462    copies of a graph `H` in a larger graph `G`.
     463
     464    It is possible to look for induced subgraphs instead, and to
     465    iterate or count the number of their occurrences.
    464466
    465467    ALGORITHM:
    466468
    467     This algorithm is a brute-force search.
    468     Let `V(H) = \{h_1,\dots,h_k\}`.  It first tries
    469     to find in `G` a possible representant of `h_1`, then a
    470     representant of `h_2` compatible with `h_1`, then
    471     a representant of `h_3` compatible with the first
     469    The algorithm is a brute-force search.  Let `V(H) =
     470    \{h_1,\dots,h_k\}`.  It first tries to find in `G` a possible
     471    representant of `h_1`, then a representant of `h_2` compatible
     472    with `h_1`, then a representant of `h_3` compatible with the first
    472473    two, etc.
    473474
    474     This way, most of the time we need to test far less than
    475     `k! \binom{|V(G)|}{k}` subsets, and hope this brute-force
    476     technique can sometimes be useful.
     475    This way, most of the time we need to test far less than `k!
     476    \binom{|V(G)|}{k}` subsets, and hope this brute-force technique
     477    can sometimes be useful.
     478    """
     479    def __init__(self, G, H, induced = False):
     480        r"""
     481        Constructor
    477482
    478     INPUT:
     483        This constructor only checks there is no inconsistency in the
     484        input : `G` and `H` are both graphs or both digraphs.
     485        """
     486        if sum([G.is_directed(), H.is_directed()]) == 1:
     487            raise ValueError("One graph can not be directed while the other is not.")
    479488
    480     - ``G``, ``H`` -- two graphs such that ``H`` is a subgraph of ``G``.
     489    def __iter__(self):
     490        r"""
     491        Returns an iterator over all the labeleld subgraphs of `G`
     492        isomorphic to `H`.
    481493
    482     - ``induced`` -- boolean (default: ``False``); whether to require that
    483       the subgraph is an induced subgraph.
     494        EXAMPLE:
    484495
    485     OUTPUT:
     496        Iterating through all the `P_3` of `P_5`::
     497       
     498            sage: from sage.graphs.generic_graph_pyx import SubgraphSearch
     499            sage: g = graphs.PathGraph(5)
     500            sage: h = graphs.PathGraph(3)
     501            sage: S = SubgraphSearch(g, h)
     502            sage: for p in S:
     503            ...      print p
     504            [0, 1, 2]
     505            [1, 2, 3]
     506            [2, 1, 0]
     507            [2, 3, 4]
     508            [3, 2, 1]
     509            [4, 3, 2]
     510        """
     511        self._initialization()
     512        return self
    486513
    487     A list of vertices inducing a copy of ``H`` in ``G``. If none is found,
    488     an empty list is returned.
     514    def cardinality(self):
     515        r"""
     516        Returns the number of labelled subgraphs of `G` isomorphic to
     517        `H`.
    489518
    490     EXAMPLES:
     519        .. NOTE::
    491520
    492     A Petersen graph contains an induced path graph `P_5`::
     521           This method counts the subgraphs by enumerating them all !
     522           Hence it probably is not a good idea to count their number
     523           before enumerating them :-)
    493524
    494         sage: from sage.graphs.generic_graph_pyx import subgraph_search
    495         sage: g = graphs.PetersenGraph()
    496         sage: subgraph_search(g, graphs.PathGraph(5), induced=True)
    497         [0, 1, 2, 3, 8]
     525        EXAMPLE:
    498526
    499     It also contains a the claw `K_{1,3}`::
     527        Counting the number of labelled `P_3` in `P_5`::
     528       
     529            sage: from sage.graphs.generic_graph_pyx import SubgraphSearch
     530            sage: g = graphs.PathGraph(5)
     531            sage: h = graphs.PathGraph(3)
     532            sage: S = SubgraphSearch(g, h)
     533            sage: S.cardinality()
     534            6
     535        """
     536        if self.nh > self.ng:
     537            return 0
    500538
    501         sage: subgraph_search(g, graphs.ClawGraph())
    502         [0, 1, 4, 5]
     539        self._initialization()
     540        cdef int i
    503541
    504     Though it contains no induced `P_6`::
     542        i=0
     543        for _ in self:
     544            i+=1
    505545
    506         sage: subgraph_search(g, graphs.PathGraph(6), induced=True)
    507         []
     546        return i
    508547
    509     TESTS:
     548    def _initialization(self):
     549        r"""
     550        Initialization of the variables.
    510551
    511     Let `G` and `H` be graphs having orders `m` and `n`, respectively. If
    512     `m < n`, then there are no copies of `H` in `G`::
     552        Once the memory allocation is done in :meth:`__cinit__`,
     553        several variables need to be set to a default value. As this
     554        operation needs to be performed before any call to
     555        :meth:`__iter__` or to :meth:`cardinality`, it is cleaner to
     556        create a dedicated method.
    513557
    514         sage: from sage.graphs.generic_graph_pyx import subgraph_search
    515         sage: m = randint(100, 200)
    516         sage: n = randint(m + 1, 300)
    517         sage: G = graphs.RandomGNP(m, random())
    518         sage: H = graphs.RandomGNP(n, random())
    519         sage: G.order() < H.order()
    520         True
    521         sage: subgraph_search(G, H)
    522         []
    523     """
    524     # TODO: This is a brute-force search and can be very inefficient. Write
    525     # a more efficient subgraph search implementation.
    526     cdef int ng = G.order()
    527     cdef int nh = H.order()
    528     if ng < nh:
    529         return []
    530     cdef int i, j, k
    531     cdef int *tmp_array
    532     cdef (bint) (*is_admissible) (int, int *, int *)
    533     if induced:
    534         is_admissible = vectors_equal
    535     else:
    536         is_admissible = vectors_inferior
    537     # static copies of the two graphs for more efficient operations
    538     cdef DenseGraph g = DenseGraph(ng)
    539     cdef DenseGraph h = DenseGraph(nh)
    540     # copying the adjacency relations in both G and H
    541     i = 0
    542     for row in G.adjacency_matrix():
    543         j = 0
    544         for k in row:
    545             if k:
    546                 g.add_arc(i, j)
    547             j += 1
    548         i += 1
    549     i = 0
    550     for row in H.adjacency_matrix():
    551         j = 0
    552         for k in row:
    553             if k:
    554                 h.add_arc(i, j)
    555             j += 1
    556         i += 1
    557     # A vertex is said to be busy if it is already part of the partial copy
    558     # of H in G.
    559     cdef int *busy = <int *>sage_malloc(ng * sizeof(int))
    560     memset(busy, 0, ng * sizeof(int))
    561     # 0 is the first vertex we use, so it is at first busy
    562     busy[0] = 1
    563     # stack -- list of the vertices which are part of the partial copy of H
    564     # in G.
    565     #
    566     # stack[i] -- the integer corresponding to the vertex of G representing
    567     # the i-th vertex of H.
    568     #
    569     # stack[i] = -1 means that i is not represented
    570     # ... yet!
    571     cdef int *stack = <int *>sage_malloc(nh * sizeof(int))
    572     stack[0] = 0
    573     stack[1] = -1
    574     # Number of representants we have already found. Set to 1 as vertex 0
    575     # is already part of the partial copy of H in G.
    576     cdef int active = 1
    577     # vertices is equal to range(nh), as an int *variable
    578     cdef int *vertices = <int *>sage_malloc(nh * sizeof(int))
    579     for 0 <= i < nh:
    580         vertices[i] = i
    581     # line_h[i] represents the adjacency sequence of vertex i
    582     # in h relative to vertices 0, 1, ..., i-1
    583     cdef int **line_h = <int **>sage_malloc(nh * sizeof(int *))
    584     for 0 <= i < nh:
    585         line_h[i] = <int *>h.adjacency_sequence(i, vertices, i)
    586     # the sequence of vertices to be returned
    587     cdef list value = []
    588    
    589     _sig_on
     558        EXAMPLE:
    590559
    591     # as long as there is a non-void partial copy of H in G
    592     while active:
    593         # If we are here and found nothing yet, we try the next possible
    594         # vertex as a representant of the active i-th vertex of H.
    595         i = stack[active] + 1
    596         # Looking for a vertex that is not busy and compatible with the
    597         # partial copy we have of H.
    598         while i < ng:
    599             if busy[i]:
    600                 i += 1
     560        Finding two times the first occurrence through the
     561        re-initialization of the instance ::
     562
     563            sage: from sage.graphs.generic_graph_pyx import SubgraphSearch
     564            sage: g = graphs.PathGraph(5)
     565            sage: h = graphs.PathGraph(3)
     566            sage: S = SubgraphSearch(g, h)
     567            sage: S.__next__()
     568            [0, 1, 2]
     569            sage: S._initialization()
     570            sage: S.__next__()
     571            [0, 1, 2]
     572        """
     573
     574        memset(self.busy, 0, self.ng * sizeof(int))
     575        # 0 is the first vertex we use, so it is at first busy
     576        self.busy[0] = 1
     577        # stack -- list of the vertices which are part of the partial copy of H
     578        # in G.
     579        #
     580        # stack[i] -- the integer corresponding to the vertex of G representing
     581        # the i-th vertex of H.
     582        #
     583        # stack[i] = -1 means that i is not represented
     584        # ... yet!
     585
     586        self.stack[0] = 0
     587        self.stack[1] = -1
     588
     589        # Number of representants we have already found. Set to 1 as vertex 0
     590        # is already part of the partial copy of H in G.
     591        self.active = 1
     592
     593    def __cinit__(self, G, H, induced = False):
     594        r"""
     595        Cython constructor
     596
     597        This method initializes all the C values.
     598        """
     599
     600        # Storing the number of vertices
     601        self.ng = G.order()
     602        self.nh = H.order()
     603
     604
     605
     606        # Storing the list of vertices
     607        self.g_vertices = G.vertices()
     608
     609        # Are the graphs directed (at the end of the code is checked
     610        # whether both are of the same type)
     611        self.directed = G.is_directed()
     612
     613        cdef int i, j, k
     614        cdef int *tmp_array
     615
     616        # Should we look for induced subgraphs ?
     617        if induced:
     618            self.is_admissible = vectors_equal
     619        else:
     620            self.is_admissible = vectors_inferior
     621
     622        # static copies of the two graphs for more efficient operations
     623        self.g = DenseGraph(self.ng)
     624        self.h = DenseGraph(self.nh)
     625
     626        # copying the adjacency relations in both G and H
     627        i = 0
     628        for row in G.adjacency_matrix():
     629            j = 0
     630            for k in row:
     631                if k:
     632                    self.g.add_arc(i, j)
     633                j += 1
     634            i += 1
     635        i = 0
     636        for row in H.adjacency_matrix():
     637            j = 0
     638            for k in row:
     639                if k:
     640                    self.h.add_arc(i, j)
     641                j += 1
     642            i += 1
     643
     644        # A vertex is said to be busy if it is already part of the partial copy
     645        # of H in G.
     646        self.busy = <int *>sage_malloc(self.ng * sizeof(int))
     647        self.stack = <int *>sage_malloc(self.nh * sizeof(int))
     648
     649        # vertices is equal to range(nh), as an int *variable
     650        self.vertices = <int *>sage_malloc(self.nh * sizeof(int))
     651        for 0 <= i < self.nh:
     652            self.vertices[i] = i
     653
     654        # line_h_out[i] represents the adjacency sequence of vertex i
     655        # in h relative to vertices 0, 1, ..., i-1
     656        self.line_h_out = <int **>sage_malloc(self.nh * sizeof(int *))
     657        for 0 <= i < self.nh:
     658            self.line_h_out[i] = <int *>self.h.adjacency_sequence_out(i, self.vertices, i)
     659
     660        # Similarly in the opposite direction (only useful if the
     661        # graphs are directed)
     662        if self.directed:
     663            self.line_h_in = <int **>sage_malloc(self.nh * sizeof(int *))
     664            for 0 <= i < self.nh:
     665                self.line_h_in[i] = <int *>self.h.adjacency_sequence_in(i, self.vertices, i)
     666
     667        self._initialization()
     668
     669    def __next__(self):
     670        r"""
     671        Returns the next isomorphic subgraph if any, and raises a
     672        ``StopIteration`` otherwise.
     673
     674        EXAMPLE::
     675
     676            sage: from sage.graphs.generic_graph_pyx import SubgraphSearch
     677            sage: g = graphs.PathGraph(5)
     678            sage: h = graphs.PathGraph(3)
     679            sage: S = SubgraphSearch(g, h)
     680            sage: S.__next__()
     681            [0, 1, 2]
     682        """
     683        _sig_on
     684        cdef int *tmp_array_out
     685        cdef int *tmp_array_in
     686        cdef bool is_admissible
     687
     688        # as long as there is a non-void partial copy of H in G
     689        while self.active >= 0:
     690            # If we are here and found nothing yet, we try the next possible
     691            # vertex as a representant of the active i-th vertex of H.
     692            self.i = self.stack[self.active] + 1
     693            # Looking for a vertex that is not busy and compatible with the
     694            # partial copy we have of H.
     695            while self.i < self.ng:
     696                if self.busy[self.i]:
     697                    self.i += 1
     698                else:
     699                    # Testing whether the vertex we picked is a
     700                    # correct extension by checking the edges from the
     701                    # vertices already selected to self.i satisfy the
     702                    # constraints
     703                    tmp_array_out = self.g.adjacency_sequence_out(self.active, self.stack, self.i)
     704                    is_admissible = self.is_admissible(self.active, tmp_array_out, self.line_h_out[self.active])
     705                    sage_free(tmp_array_out)
     706
     707                    # If G and H are digraphs, we also need to ensure
     708                    # the edges going in the opposite direction
     709                    # satisfy the constraints
     710                    if is_admissible and self.directed:
     711                        tmp_array_in = self.g.adjacency_sequence_in(self.active, self.stack, self.i)
     712                        is_admissible = is_admissible and self.is_admissible(self.active, tmp_array_in, self.line_h_in[self.active])
     713                        sage_free(tmp_array_in)
     714
     715                    if is_admissible:
     716                        break
     717                    else:
     718                        self.i += 1                       
     719
     720            # If we have found a good representant of H's i-th vertex in G
     721            if self.i < self.ng:
     722
     723                # updating the last vertex of the stack
     724                if self.stack[self.active] != -1:
     725                    self.busy[self.stack[self.active]] = 0
     726                self.stack[self.active] = self.i
     727
     728                # We have found our copy !!!
     729                if self.active == self.nh-1:
     730                    return [self.g_vertices[self.stack[l]] for l in xrange(self.nh)]
     731
     732                # We are still missing several vertices ...
     733                else:
     734                    self.busy[self.stack[self.active]] = 1
     735                    self.active += 1
     736
     737                    # we begin the search of the next vertex at 0
     738                    self.stack[self.active] = -1
     739
     740            # If we found no representant for the i-th vertex, it
     741            # means that we cannot extend the current copy of H so we
     742            # update the status of stack[active] and prepare to change
     743            # the previous vertex.
     744
    601745            else:
    602                 tmp_array = g.adjacency_sequence(active, stack, i)
    603                 if is_admissible(active, tmp_array, line_h[active]):
    604                     sage_free(tmp_array)
    605                     break
    606                 else:
    607                     sage_free(tmp_array)
    608                     i += 1
    609         # If we found none, it means that we cannot extend the current copy
    610         # of H so we update the status of stack[active] and prepare to change
    611         # the previous vertex.
    612         if i >= ng:
    613             if stack[active] != -1:
    614                 busy[stack[active]] = 0
    615             stack[active] = -1
    616             active -= 1
    617         # If we have found a good representant of H's i-th vertex in G
    618         else:
    619             if stack[active] != -1:
    620                 busy[stack[active]] = 0
    621             stack[active] = i
    622             busy[stack[active]] = 1
    623             active += 1
    624             # We have found our copy!!!
    625             if active == nh:
    626                 g_vertices = G.vertices()
    627                 value = [g_vertices[stack[i]] for i in xrange(nh)]
    628                 break
    629             else:
    630                 # we begin the search of the next vertex at 0
    631                 stack[active] = -1
     746                if self.stack[self.active] != -1:
     747                    self.busy[self.stack[self.active]] = 0
     748                self.stack[self.active] = -1
     749                self.active -= 1
    632750
    633     _sig_off
     751        _sig_off
     752        raise StopIteration
    634753
    635     # Free the memory
    636     sage_free(busy)
    637     sage_free(stack)
    638     sage_free(vertices)
    639     for 0 <= i < nh:
    640         sage_free(line_h[i])
    641     sage_free(line_h)
     754    def __dealloc__(self):
     755        r"""
     756        Freeing the allocated memory.
     757        """
    642758
    643     return value
     759        # Free the memory
     760        sage_free(self.busy)
     761        sage_free(self.stack)
     762        sage_free(self.vertices)
     763        for 0 <= i < self.nh:
     764            sage_free(self.line_h_out[i])
     765        sage_free(self.line_h_out)
     766
     767        if self.directed:
     768            for 0 <= i < self.nh:
     769                sage_free(self.line_h_in[i])
     770            sage_free(self.line_h_in)
    644771
    645772cdef inline bint vectors_equal(int n, int *a, int *b):
    646773    r"""