# HG changeset patch
# User Nathann Cohen
# Date 1278157947 7200
# Node ID 20b8ec1a5d34e2633a8d6ce057f4f19b0165fd86
# Parent 113df1e3b94772d3d5cf45f213f48a0c480a25ec
trac 9420  SubgraphSearch class
diff r 113df1e3b947 r 20b8ec1a5d34 sage/graphs/base/c_graph.pxd
 a/sage/graphs/base/c_graph.pxd Thu Jul 01 13:20:00 2010 +0200
+++ b/sage/graphs/base/c_graph.pxd Sat Jul 03 13:52:27 2010 +0200
@@ 28,7 +28,8 @@
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)
diff r 113df1e3b947 r 20b8ec1a5d34 sage/graphs/base/c_graph.pyx
 a/sage/graphs/base/c_graph.pyx Thu Jul 01 13:20:00 2010 +0200
+++ b/sage/graphs/base/c_graph.pyx Sat Jul 03 13:52:27 2010 +0200
@@ 840,38 +840,100 @@
"""
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 ith element
 is set to 1 iff ``v`` is adjacent to ``vertices[i]``.
+ Returns a list of ``n`` integers, whose ith 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 = 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 ith 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 = 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``.
diff r 113df1e3b947 r 20b8ec1a5d34 sage/graphs/base/dense_graph.pyx
 a/sage/graphs/base/dense_graph.pyx Thu Jul 01 13:20:00 2010 +0200
+++ b/sage/graphs/base/dense_graph.pyx Sat Jul 03 13:52:27 2010 +0200
@@ 625,15 +625,15 @@
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,7 +658,7 @@
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])
diff r 113df1e3b947 r 20b8ec1a5d34 sage/graphs/base/sparse_graph.pyx
 a/sage/graphs/base/sparse_graph.pyx Thu Jul 01 13:20:00 2010 +0200
+++ b/sage/graphs/base/sparse_graph.pyx Sat Jul 03 13:52:27 2010 +0200
@@ 1306,15 +1306,15 @@
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,7 +1341,7 @@
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])
diff r 113df1e3b947 r 20b8ec1a5d34 sage/graphs/generic_graph.py
 a/sage/graphs/generic_graph.py Thu Jul 01 13:20:00 2010 +0200
+++ b/sage/graphs/generic_graph.py Sat Jul 03 13:52:27 2010 +0200
@@ 8008,6 +8008,10 @@
is no copy (induced or otherwise) of ``G`` in ``self``, we return
``None``.
+ .. NOTE::
+
+ This method also works on digraphs.
+
ALGORITHM:
Bruteforce search.
@@ 8058,20 +8062,145 @@
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:
+
+ Bruteforce 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, ..., k1\}` such that `ij \in
+ T_k` iif `i 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 bruteforce 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 = 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 ith vertex of H.
 #
 # stack[i] = 1 means that i is not represented
 # ... yet!
 cdef int *stack = 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 = 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, ..., i1
 cdef int **line_h = sage_malloc(nh * sizeof(int *))
 for 0 <= i < nh:
 line_h[i] = 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 nonvoid 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 ith 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
+ reinitialization 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 ith 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 = sage_malloc(self.ng * sizeof(int))
+ self.stack = sage_malloc(self.nh * sizeof(int))
+
+ # vertices is equal to range(nh), as an int *variable
+ self.vertices = 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, ..., i1
+ self.line_h_out = sage_malloc(self.nh * sizeof(int *))
+ for 0 <= i < self.nh:
+ self.line_h_out[i] = 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 = sage_malloc(self.nh * sizeof(int *))
+ for 0 <= i < self.nh:
+ self.line_h_in[i] = 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 nonvoid 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 ith 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 ith 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.nh1:
+ 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 ith 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 ith 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"""