Ticket #10885: trac_10885.patch

sage/graphs/base/c_graph.pyx

@@ -2851 +2851 @@
         else:
             return True
 
+def floyd_warshall(gg, paths = True, distances = False):
+    r""" 
+    Computes the shortest path/distances between all pairs of vertices.
+
+    For more information on the Floyd-Warshall algorithm, see the `Wikipedia
+    article on Floyd-Warshall
+    <http://en.wikipedia.org/wiki/Floyd%E2%80%93Warshall_algorithm>`_.
+
+    INPUT:
+
+        - ``gg`` -- the graph on which to work.
+
+        - ``paths`` (boolean) -- whether to return the dictionary of shortest
+          paths. Set to ``True`` by default.
+
+        - ``distances`` (boolean) -- whether to return the dictionary of
+          distances. Set to ``False`` by default.
+
+    OUTPUT:
+
+        Depending on the input, this function return the dictionary of paths,
+        the dictionary of distances, or 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. Let us also remember
+        that if the memory size is quadratic, the algorithm runs in cubic time.
+
+    EXAMPLES:
+
+    Shortest paths in a small grid ::
+
+        sage: g = graphs.Grid2dGraph(2,2)
+        sage: from sage.graphs.base.c_graph import floyd_warshall
+        sage: print floyd_warshall(g)
+        {(0, 1): {(0, 1): None, (1, 0): (0, 0), (0, 0): (0, 1), (1, 1): (0, 1)},
+        (1, 0): {(0, 1): (0, 0), (1, 0): None, (0, 0): (1, 0), (1, 1): (1, 0)},
+        (0, 0): {(0, 1): (0, 0), (1, 0): (0, 0), (0, 0): None, (1, 1): (0, 1)},
+        (1, 1): {(0, 1): (1, 1), (1, 0): (1, 1), (0, 0): (0, 1), (1, 1): None}}
+
+    Checking the distances are correct ::
+
+        sage: g = graphs.Grid2dGraph(5,5)
+        sage: dist,path = floyd_warshall(g, distances = True)
+        sage: all( dist[u][v] == g.distance(u,v) for u in g for v in g )
+        True
+
+    Checking a random path is valid ::
+
+        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
+    """
+    from sage.rings.infinity import Infinity    
+    cdef CGraph g = <CGraph> gg._backend._cg
+    cdef int n = max(g.verts())+1
+
+    if n > <unsigned short> -1:
+        raise ValueError("The graph backend contains more than "+(<unsigned short> -1)+" nodes")
+
+    # Dictionaries of distance, precedent element, and integers
+    cdef dict d_prec = dict()
+    cdef dict d_dist = dict()
+    cdef dict tmp
+
+    cdef int v_int
+    cdef int u_int
+    cdef int w_int
+    cdef int i
+
+    # All this just creates two tables prec[n][n] and dist[n][n]
+    cdef unsigned short * t_prec = <unsigned short *> sage_malloc(n*n*sizeof(short))
+    cdef unsigned short * t_dist = <unsigned short *> sage_malloc(n*n*sizeof(short))
+    cdef unsigned short ** prec = <unsigned short **> sage_malloc(n*sizeof(short *))
+    cdef unsigned short ** dist = <unsigned short **> sage_malloc(n*sizeof(short *))
+    prec[0] = t_prec
+    dist[0] = t_dist
+
+    for 1 <= i< n:
+        prec[i] = prec[i-1] + n
+        dist[i] = dist[i-1] + n
+
+    # Initializing prec and dist
+    memset(t_prec, 0, n*n*sizeof(short))
+    memset(t_dist, -1, n*n*sizeof(short))
+
+    # Copying the adjacency matrix (vertices at distance 1)
+    for v_int in g.verts():
+        prec[v_int][v_int] = v_int     
+        dist[v_int][v_int] =  0
+        for u_int in g.out_neighbors(v_int):
+            dist[v_int][u_int] = 1
+            prec[v_int][u_int] = v_int
+
+    # The algorithm itself.
+
+    for w_int in g.verts():
+        for v_int in g.verts():
+            for u_int in g.verts():
+
+                # If it is shorter to go from u to v through w, do it
+                if dist[v_int][u_int] > dist[v_int][w_int] + dist[w_int][u_int]:
+                    dist[v_int][u_int] = dist[v_int][w_int] + dist[w_int][u_int]
+                    prec[v_int][u_int] = prec[w_int][u_int]
+
+    # If the paths are to be returned
+    if paths:
+        for v_int in g.verts():
+            tmp = dict()
+            v = vertex_label(v_int, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
+
+            for u_int in g.verts():
+                u = vertex_label(u_int, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
+                w = (None if v == u
+                     else vertex_label(prec[v_int][u_int], gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg))
+                tmp[u] = w
+                
+            d_prec[v] = tmp
+
+    sage_free(t_prec)
+    sage_free(prec)
+
+    # If the distances are to be returned
+    if distances:
+        for v_int in g.verts():
+            tmp = dict()
+            v = vertex_label(v_int, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
+
+            for u_int in g.verts():
+                u = vertex_label(u_int, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
+
+                tmp[u] = dist[v_int][u_int] if (dist[v_int][u_int] != <unsigned short> -1) else Infinity
+                
+            d_dist[v] = tmp
+
+    sage_free(t_dist)
+    sage_free(dist)
+
+    if distances and paths:
+        return d_dist, d_prec
+    elif paths:
+        return d_prec
+    else:
+        return d_dist
+
 cdef class Search_iterator:
     r"""
     An iterator for traversing a (di)graph.

sage/graphs/generic_graph.py

@@ -10092 +10092 @@
         """
         return self.shortest_path_length(u, v)
 
-    def distance_all_pairs(self):
+    def distance_all_pairs(self, algorithm = "BFS"):
         r"""
         Returns the distances between all pairs of vertices.
 
+        INPUT:
+
+            - ``"algorithm"`` (string) -- two algorithms are available
+
+                * ``algorithm = "BFS"`` in which case the distances are computed
+                  through `n-1` different breadth-first-search.
+
+                * ``algorithm = "Floyd-Warshall"``, in which case the
+                  Floyd-Warshall algorithm is used.
+
+                The default is ``algorithm = "BFS"``.
+
         OUTPUT:
 
         A doubly indexed dictionary
@@ -10115 +10127 @@
             sage: all([g.distance(0,v) == distances[0][v] for v in g])
             True
         """
-
-        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
+        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
+
+        elif algorithm == "Floyd-Warshall":
+            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 \"Floyd-Warshall\"")
 
     def eccentricity(self, v=None, dist_dict=None, with_labels=False):
         """
@@ -10916 +10935 @@
                 lengths[v] = len(paths[v]) - 1
             return lengths
 
-    def shortest_path_all_pairs(self, by_weight=True, default_weight=1):
+    def shortest_path_all_pairs(self, by_weight=False, default_weight=1):
         """
         Uses the Floyd-Warshall algorithm to find a shortest weighted path
         for each pair of vertices.
         
-        The weights (labels) on the vertices can be anything that can be
-        compared and can be summed.
-        
-        INPUT:
-        
-        
-        -  ``by_weight`` - If False, figure distances by the
-           numbers of edges.
-        
-        -  ``default_weight`` - (defaults to 1) The default
-           weight to assign edges that don't have a weight (i.e., a label).
-        
-        
-        OUTPUT: A tuple (dist, pred). They are both dicts of dicts. The
-        first indicates the length dist[u][v] of the shortest weighted path
-        from u to v. The second is more complicated- it indicates the
-        predecessor pred[u][v] of v in the shortest path from u to v.
+        INPUT:
+        
+        
+        - ``by_weight`` - Whether to use the labels defined over the edges as
+           weights. If ``False`` (default), the distance between `u` and `v` is
+           the minimum number of edges of a path from `u` to `v`.
+        
+        - ``default_weight`` - (defaults to 1) The default weight to assign
+           edges that don't have a weight (i.e., a label).
+
+           Implies ``by_weight == True``.
+        
+        OUTPUT: 
+
+            A tuple ``(dist, pred)``. They are both dicts of dicts. The first
+            indicates the length ``dist[u][v]`` of the shortest weighted path
+            from `u` to `v`. The second is a compact representation of all the
+            paths- it indicates the predecessor ``pred[u][v]`` of `v` in the
+            shortest path from `u` to `v`.
+
+        .. NOTE::
+
+            Two version of the Floyd-Warshall 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 Floyd-Warshall algorithm).
         
         EXAMPLES::
         
             sage: G = Graph( { 0: {1: 1}, 1: {2: 1}, 2: {3: 1}, 3: {4: 2}, 4: {0: 2} }, sparse=True )
             sage: G.plot(edge_labels=True).show() # long time
-            sage: dist, pred = G.shortest_path_all_pairs()
+            sage: dist, pred = G.shortest_path_all_pairs(by_weight = True)
             sage: dist
             {0: {0: 0, 1: 1, 2: 2, 3: 3, 4: 2}, 1: {0: 1, 1: 0, 2: 1, 3: 2, 4: 3}, 2: {0: 2, 1: 1, 2: 0, 3: 1, 4: 3}, 3: {0: 3, 1: 2, 2: 1, 3: 0, 4: 2}, 4: {0: 2, 1: 3, 2: 3, 3: 2, 4: 0}}
             sage: pred
@@ -10951 +10985 @@
             sage: pred[0]
             {0: None, 1: 0, 2: 1, 3: 2, 4: 0}
         
-        So for example the shortest weighted path from 0 to 3 is obtained
-        as follows. The predecessor of 3 is pred[0][3] == 2, the
-        predecessor of 2 is pred[0][2] == 1, and the predecessor of 1 is
-        pred[0][1] == 0.
+        So for example the shortest weighted path from `0` to `3` is obtained as
+        follows. The predecessor of `3` is ``pred[0][3] == 2``, the predecessor
+        of `2` is ``pred[0][2] == 1``, and the predecessor of `1` is
+        ``pred[0][1] == 0``.
         
         ::
         
             sage: G = Graph( { 0: {1:None}, 1: {2:None}, 2: {3: 1}, 3: {4: 2}, 4: {0: 2} }, sparse=True )
-            sage: G.shortest_path_all_pairs(by_weight=False)
+            sage: G.shortest_path_all_pairs()
             ({0: {0: 0, 1: 1, 2: 2, 3: 2, 4: 1},
             1: {0: 1, 1: 0, 2: 1, 3: 2, 4: 2},
             2: {0: 2, 1: 1, 2: 0, 3: 1, 4: 2},
@@ -10970 +11004 @@
             2: {0: 1, 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}})
-            sage: G.shortest_path_all_pairs()
+            sage: G.shortest_path_all_pairs(by_weight = True)
             ({0: {0: 0, 1: 1, 2: 2, 3: 3, 4: 2},
             1: {0: 1, 1: 0, 2: 1, 3: 2, 4: 3},
             2: {0: 2, 1: 1, 2: 0, 3: 1, 4: 3},
@@ -10993 +11027 @@
             3: {0: 4, 1: 2, 2: 3, 3: None, 4: 3},
             4: {0: 4, 1: 0, 2: 3, 3: 4, 4: None}})
         """
+        if default_weight != 1:
+            by_weight = True
+
+        if by_weight is False:
+            from sage.graphs.base.c_graph import floyd_warshall
+            return floyd_warshall(self, distances = True)
+
         from sage.rings.infinity import Infinity
         dist = {}
         pred = {}