# HG changeset patch
# User Nathann Cohen
# Date 1299773000 28800
# Node ID d026f22e0d50bfc9ca1b06b0944bc5638a5c735e
# Parent b0a5118e9d0af16c871460c37f02d2c1dc28cc28
trac 10905  shortest path all pairs through BFS computations
diff r b0a5118e9d0a r d026f22e0d50 sage/graphs/base/c_graph.pyx
 a/sage/graphs/base/c_graph.pyx Sun Mar 06 22:14:10 2011 +0100
+++ b/sage/graphs/base/c_graph.pyx Fri Mar 11 00:03:20 2011 +0800
@@ 3010,6 +3010,123 @@
else:
return d_dist
+cpdef all_pairs_shortest_path_BFS(gg):
+ """
+ Computes the shortest path and the distances between each pair of vertices
+ through successive breadthfirstsearches
+
+ OUTPUT:
+
+ This function return a pair of dictionaries ``(distances,paths)`` where
+ ``distance[u][v]`` denotes the distance of a shortest path from `u` to
+ `v` and ``paths[u][v]`` denotes an inneighbor `w` of `v` such that
+ `dist(u,v)= 1 + dist(u,w)`.
+
+ .. WARNING::
+
+ Because this function works on matrices whose size is quadratic compared
+ to the number of vertices, it uses short variables instead of long ones
+ to divide by 2 the size in memory. This means that the current
+ implementation does not run on a graph of more than 65536 nodes (this
+ can be easily changed if necessary, but would require much more
+ memory. It may be worth writing two versions). For information, the
+ current version of the algorithm on a graph with `65536=2^{16}` nodes
+ creates in memory `2` tables on `2^{32}` short elements (2bytes each),
+ for a total of `2^{34}` bytes or `16` gigabytes.
+
+ EXAMPLE:
+
+ On a grid::
+
+ sage: g = graphs.Grid2dGraph(10,10)
+ sage: from sage.graphs.base.c_graph import all_pairs_shortest_path_BFS
+ sage: dist, path = all_pairs_shortest_path_BFS(g)
+ sage: all( dist[u][v] == g.distance(u,v) for u in g for v in g )
+ True
+ """
+ from sage.rings.infinity import Infinity
+
+ cdef CGraph cg = gg._backend._cg
+
+ cdef list vertices = cg.verts()
+ cdef int n = max(vertices)+1
+
+ if n > 1:
+ raise ValueError("The graph backend contains more than "+( 1)+" nodes")
+
+ # The vertices which have already been visited
+ cdef bitset_t seen
+ bitset_init(seen, cg.active_vertices.size)
+
+ # The list of waiting vertices, the beginning and the end of the list
+
+ cdef unsigned short * waiting_list = sage_malloc(n*sizeof(short))
+ cdef unsigned short waiting_beginning = 0
+ cdef unsigned short waiting_end = 0
+
+ cdef unsigned short source
+ cdef unsigned short v, u
+
+ # Dictionaries of dictionaries of distances/predecessors
+ cdef dict d_distances = dict()
+ cdef dict d_prec = dict()
+
+ # Temporary dictionaries of distances/predecessors
+ cdef dict tmp_distances
+ cdef dict tmp_prec
+
+ # Temporary vectors of distances/predecessors
+ cdef unsigned short * v_distances = sage_malloc(n*sizeof(short))
+ cdef unsigned short * v_prec = sage_malloc(n*sizeof(short))
+
+ for source in vertices:
+ bitset_set_first_n(seen, 0)
+ bitset_add(seen, source)
+
+ v_distances[source] = 0
+ v_prec[source] = source
+ waiting_list[0] = source
+ waiting_beginning = 0
+ waiting_end = 0
+
+ while waiting_beginning <= waiting_end:
+ v = waiting_list[waiting_beginning]
+
+ for u in cg.out_neighbors(v):
+ if not bitset_in(seen, u):
+ v_distances[u] = v_distances[v]+1
+ v_prec[u] = v
+ bitset_add(seen, u)
+ waiting_end += 1
+ waiting_list[waiting_end] = u
+
+ waiting_beginning += 1
+
+ tmp_distances = dict()
+ tmp_prec = dict()
+ for v in vertices:
+ vv = vertex_label(v, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
+
+ if bitset_in(seen, v):
+ tmp_prec[vv] = vertex_label(v_prec[v], gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
+ tmp_distances[vv] = v_distances[v]
+ else:
+ tmp_prec[vv] = None
+ tmp_distances[vv] = Infinity
+
+ vv = vertex_label(source, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
+ tmp_prec[vv] = None
+ d_prec[vv] = tmp_prec
+ d_distances[vv] = tmp_distances
+
+ bitset_free(seen)
+ sage_free(waiting_list)
+ sage_free(v_distances)
+ sage_free(v_prec)
+
+ return d_distances, d_prec
+
+
cdef class Search_iterator:
r"""
An iterator for traversing a (di)graph.
diff r b0a5118e9d0a r d026f22e0d50 sage/graphs/generic_graph.py
 a/sage/graphs/generic_graph.py Sun Mar 06 22:14:10 2011 +0100
+++ b/sage/graphs/generic_graph.py Fri Mar 11 00:03:20 2011 +0800
@@ 10092,7 +10092,7 @@
"""
return self.shortest_path_length(u, v)
 def distance_all_pairs(self, algorithm = "BFS"):
+ def distance_all_pairs(self, algorithm = "auto"):
r"""
Returns the distances between all pairs of vertices.
@@ 10101,11 +10101,15 @@
 ``"algorithm"`` (string)  two algorithms are available
* ``algorithm = "BFS"`` in which case the distances are computed
 through `n1` different breadthfirstsearch.
+ through `n` different breadthfirstsearch.
* ``algorithm = "FloydWarshall"``, in which case the
FloydWarshall algorithm is used.
+ * ``algorithm = "auto"``, in which case the FloydWarshall
+ algorithm is used for graphs on less than 20 vertices, and BFS
+ otherwise.
+
The default is ``algorithm = "BFS"``.
OUTPUT:
@@ 10127,27 +10131,22 @@
sage: all([g.distance(0,v) == distances[0][v] for v in g])
True
"""
+ if algorithm == "auto":
+ if self.order() <= 20:
+ algorithm = "FloydWarshall"
+ else:
+ algorithm = "BFS"
+
if algorithm == "BFS":
 from sage.rings.infinity import Infinity
 distances = dict([(v, self.shortest_path_lengths(v)) for v in self])

 # setting the necessary +Infinity
 cc = self.connected_components()
 for cc1 in cc:
 for cc2 in cc:
 if cc1 != cc2:
 for u in cc1:
 for v in cc2:
 distances[u][v] = Infinity

 return distances
+ from sage.graphs.base.c_graph import all_pairs_shortest_path_BFS
+ return all_pairs_shortest_path_BFS(self)[0]
elif algorithm == "FloydWarshall":
from sage.graphs.base.c_graph import floyd_warshall
return floyd_warshall(self,paths = False, distances = True)
else:
 raise ValueError("The algorithm keyword can be equal to either \"BFS\" or \"FloydWarshall\"")
+ raise ValueError("The algorithm keyword can be equal to either \"BFS\" or \"FloydWarshall\" or \"auto\"")
def eccentricity(self, v=None, dist_dict=None, with_labels=False):
"""
@@ 10935,11 +10934,10 @@
lengths[v] = len(paths[v])  1
return lengths
 def shortest_path_all_pairs(self, by_weight=False, default_weight=1):
 """
 Uses the FloydWarshall algorithm to find a shortest weighted path
 for each pair of vertices.

+ def shortest_path_all_pairs(self, by_weight=False, default_weight=1, algorithm = "auto"):
+ """
+ Computes a shortest path between each pair of vertices.
+
INPUT:
@@ 10951,6 +10949,23 @@
edges that don't have a weight (i.e., a label).
Implies ``by_weight == True``.
+
+  ``algorithm``  four options :
+
+ * ``"BFS"``  the computation is done through a BFS
+ centered on each vertex successively. Only implemented
+ when ``default_weight = 1`` and ``by_weight = False``.
+
+ * ``"FloydWarshallCython"``  through the Cython implementation of
+ the FloydWarshall algorithm.
+
+ * ``"FloydWarshallPython"``  through the Python implementation of
+ the FloydWarshall algorithm.
+
+ * ``"auto"``  use the fastest algorithm depending on the input
+ (``"BFS"`` if possible, and ``"FloydWarshallPython"`` otherwise)
+
+ This is the default value.
OUTPUT:
@@ 10962,16 +10977,19 @@
.. NOTE::
 Two version of the FloydWarshall algorithm are implemented. One is
 pure Python, the other one is Cython. The first is used whenever
 ``by_weight = True``, and can perform its computations using the
 labels on the edges for as long as they can be added together. The
 second one, used when ``by_weight = False``, is much faster but only
 deals with the topological distance (each edge is of weight 1, or
 equivalently the length of a path is its number of edges). Besides,
 the current version of this second implementation does not deal with
 graphs larger than 65536 vertices (which already represents 16GB of
 ram when running the FloydWarshall algorithm).
+ Three different implementations are actually available through this method :
+
+ * BFS (Cython)
+ * FloydWarshall (Cython)
+ * FloydWarshall (Python)
+
+ The BFS algorithm is the fastest of the three, then comes the Cython
+ implementation of FloydWarshall, and last the Python
+ implementation. The first two implementations, however, only compute
+ distances based on the topological distance (each edge is of weight
+ 1, or equivalently the length of a path is its number of
+ edges). Besides, they do not deal with graphs larger than 65536
+ vertices (which already represents 16GB of ram).
EXAMPLES::
@@ 11026,14 +11044,50 @@
2: {0: 4, 1: 2, 2: None, 3: 2, 4: 3},
3: {0: 4, 1: 2, 2: 3, 3: None, 4: 3},
4: {0: 4, 1: 0, 2: 3, 3: 4, 4: None}})
+
+ Checking the distances are equal regardless of the algorithm used::
+
+ sage: g = graphs.Grid2dGraph(5,5)
+ sage: d1, _ = g.shortest_path_all_pairs(algorithm="BFS")
+ sage: d2, _ = g.shortest_path_all_pairs(algorithm="FloydWarshallCython")
+ sage: d3, _ = g.shortest_path_all_pairs(algorithm="FloydWarshallPython")
+ sage: d1 == d2 == d3
+ True
+
+ Checking a random path is valid ::
+
+ sage: dist, path = g.shortest_path_all_pairs(algorithm="BFS")
+ sage: u,v = g.random_vertex(), g.random_vertex()
+ sage: p = [v]
+ sage: while p[0] != None:
+ ... p.insert(0,path[u][p[0]])
+ sage: len(p) == dist[u][v] + 2
+ True
"""
if default_weight != 1:
by_weight = True
 if by_weight is False:
+ if algorithm == "auto":
+ if by_weight is False:
+ algorithm = "BFS"
+ else:
+ algorithm = "FloydWarshallPython"
+
+ if algorithm == "BFS":
+ from sage.graphs.base.c_graph import all_pairs_shortest_path_BFS
+ return all_pairs_shortest_path_BFS(self)
+
+ elif algorithm == "FloydWarshallCython":
from sage.graphs.base.c_graph import floyd_warshall
return floyd_warshall(self, distances = True)
+ elif algorithm != "FloydWarshallPython":
+ raise ValueError("The algorithm keyword can only be set to "+
+ "\"auto\","+
+ " \"BFS\", "+
+ "\"FloydWarshallCython\" or "+
+ "\"FloydWarshallCython\"")
+
from sage.rings.infinity import Infinity
dist = {}
pred = {}