# 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 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

diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/base/c_graph.pyx
 a """ 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 = 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 = 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

diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/base/dense_graph.pyx
 a 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 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

diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/base/sparse_graph.pyx
 a 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 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

diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/generic_graph.py
 a is no copy (induced or otherwise) of G in self, we return None. .. NOTE:: This method also works on digraphs. ALGORITHM: Brute-force search. 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
• ## sage/graphs/generic_graph_pyx.pxd

diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/generic_graph_pyx.pxd
 a 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

diff -r 113df1e3b947 -r 20b8ec1a5d34 sage/graphs/generic_graph_pyx.pyx
 a # 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) # 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 = 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 = 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, ..., i-1 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 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 = 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, ..., i-1 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 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"""