Ticket #10885: trac_10885.patch

File trac_10885.patch, 14.1 KB (added by ncohen, 10 years ago)
  • sage/graphs/base/c_graph.pyx

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1299446050 -3600
    # Node ID 954991d687ee27a6abc8cbbbb02cd96c7177563e
    # Parent  361a4ad7d52c69b64ae2e658ffd0820af0d87e93
    Ticket #10885 - Cython implementation of Floyd-Warshall algorithm
    
    diff -r 361a4ad7d52c -r 954991d687ee sage/graphs/base/c_graph.pyx
    a b  
    28512851        else:
    28522852            return True
    28532853
     2854def floyd_warshall(gg, paths = True, distances = False):
     2855    r"""
     2856    Computes the shortest path/distances between all pairs of vertices.
     2857
     2858    For more information on the Floyd-Warshall algorithm, see the `Wikipedia
     2859    article on Floyd-Warshall
     2860    <http://en.wikipedia.org/wiki/Floyd%E2%80%93Warshall_algorithm>`_.
     2861
     2862    INPUT:
     2863
     2864        - ``gg`` -- the graph on which to work.
     2865
     2866        - ``paths`` (boolean) -- whether to return the dictionary of shortest
     2867          paths. Set to ``True`` by default.
     2868
     2869        - ``distances`` (boolean) -- whether to return the dictionary of
     2870          distances. Set to ``False`` by default.
     2871
     2872    OUTPUT:
     2873
     2874        Depending on the input, this function return the dictionary of paths,
     2875        the dictionary of distances, or a pair of dictionaries
     2876        ``(distances,paths)`` where ``distance[u][v]`` denotes the distance of a
     2877        shortest path from `u` to `v` and ``paths[u][v]`` denotes an inneighbor
     2878        `w` of `v` such that `dist(u,v)= 1 + dist(u,w)`.
     2879
     2880    .. WARNING::
     2881
     2882        Because this function works on matrices whose size is quadratic compared
     2883        to the number of vertices, it uses short variables instead of long ones
     2884        to divide by 2 the size in memory. This means that the current
     2885        implementation does not run on a graph of more than 65536 nodes (this
     2886        can be easily changed if necessary, but would require much more
     2887        memory. It may be worth writing two versions). For information, the
     2888        current version of the algorithm on a graph with `65536=2^{16}` nodes
     2889        creates in memory `2` tables on `2^{32}` short elements (2bytes each),
     2890        for a total of `2^{34}` bytes or `16` gigabytes. Let us also remember
     2891        that if the memory size is quadratic, the algorithm runs in cubic time.
     2892
     2893    EXAMPLES:
     2894
     2895    Shortest paths in a small grid ::
     2896
     2897        sage: g = graphs.Grid2dGraph(2,2)
     2898        sage: from sage.graphs.base.c_graph import floyd_warshall
     2899        sage: print floyd_warshall(g)
     2900        {(0, 1): {(0, 1): None, (1, 0): (0, 0), (0, 0): (0, 1), (1, 1): (0, 1)},
     2901        (1, 0): {(0, 1): (0, 0), (1, 0): None, (0, 0): (1, 0), (1, 1): (1, 0)},
     2902        (0, 0): {(0, 1): (0, 0), (1, 0): (0, 0), (0, 0): None, (1, 1): (0, 1)},
     2903        (1, 1): {(0, 1): (1, 1), (1, 0): (1, 1), (0, 0): (0, 1), (1, 1): None}}
     2904
     2905    Checking the distances are correct ::
     2906
     2907        sage: g = graphs.Grid2dGraph(5,5)
     2908        sage: dist,path = floyd_warshall(g, distances = True)
     2909        sage: all( dist[u][v] == g.distance(u,v) for u in g for v in g )
     2910        True
     2911
     2912    Checking a random path is valid ::
     2913
     2914        sage: u,v = g.random_vertex(), g.random_vertex()
     2915        sage: p = [v]
     2916        sage: while p[0] != None:
     2917        ...     p.insert(0,path[u][p[0]])
     2918        sage: len(p) == dist[u][v] + 2
     2919        True
     2920    """
     2921    from sage.rings.infinity import Infinity   
     2922    cdef CGraph g = <CGraph> gg._backend._cg
     2923    cdef int n = max(g.verts())+1
     2924
     2925    if n > <unsigned short> -1:
     2926        raise ValueError("The graph backend contains more than "+(<unsigned short> -1)+" nodes")
     2927
     2928    # Dictionaries of distance, precedent element, and integers
     2929    cdef dict d_prec = dict()
     2930    cdef dict d_dist = dict()
     2931    cdef dict tmp
     2932
     2933    cdef int v_int
     2934    cdef int u_int
     2935    cdef int w_int
     2936    cdef int i
     2937
     2938    # All this just creates two tables prec[n][n] and dist[n][n]
     2939    cdef unsigned short * t_prec = <unsigned short *> sage_malloc(n*n*sizeof(short))
     2940    cdef unsigned short * t_dist = <unsigned short *> sage_malloc(n*n*sizeof(short))
     2941    cdef unsigned short ** prec = <unsigned short **> sage_malloc(n*sizeof(short *))
     2942    cdef unsigned short ** dist = <unsigned short **> sage_malloc(n*sizeof(short *))
     2943    prec[0] = t_prec
     2944    dist[0] = t_dist
     2945
     2946    for 1 <= i< n:
     2947        prec[i] = prec[i-1] + n
     2948        dist[i] = dist[i-1] + n
     2949
     2950    # Initializing prec and dist
     2951    memset(t_prec, 0, n*n*sizeof(short))
     2952    memset(t_dist, -1, n*n*sizeof(short))
     2953
     2954    # Copying the adjacency matrix (vertices at distance 1)
     2955    for v_int in g.verts():
     2956        prec[v_int][v_int] = v_int     
     2957        dist[v_int][v_int] =  0
     2958        for u_int in g.out_neighbors(v_int):
     2959            dist[v_int][u_int] = 1
     2960            prec[v_int][u_int] = v_int
     2961
     2962    # The algorithm itself.
     2963
     2964    for w_int in g.verts():
     2965        for v_int in g.verts():
     2966            for u_int in g.verts():
     2967
     2968                # If it is shorter to go from u to v through w, do it
     2969                if dist[v_int][u_int] > dist[v_int][w_int] + dist[w_int][u_int]:
     2970                    dist[v_int][u_int] = dist[v_int][w_int] + dist[w_int][u_int]
     2971                    prec[v_int][u_int] = prec[w_int][u_int]
     2972
     2973    # If the paths are to be returned
     2974    if paths:
     2975        for v_int in g.verts():
     2976            tmp = dict()
     2977            v = vertex_label(v_int, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
     2978
     2979            for u_int in g.verts():
     2980                u = vertex_label(u_int, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
     2981                w = (None if v == u
     2982                     else vertex_label(prec[v_int][u_int], gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg))
     2983                tmp[u] = w
     2984               
     2985            d_prec[v] = tmp
     2986
     2987    sage_free(t_prec)
     2988    sage_free(prec)
     2989
     2990    # If the distances are to be returned
     2991    if distances:
     2992        for v_int in g.verts():
     2993            tmp = dict()
     2994            v = vertex_label(v_int, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
     2995
     2996            for u_int in g.verts():
     2997                u = vertex_label(u_int, gg._backend.vertex_ints, gg._backend.vertex_labels, gg._backend._cg)
     2998
     2999                tmp[u] = dist[v_int][u_int] if (dist[v_int][u_int] != <unsigned short> -1) else Infinity
     3000               
     3001            d_dist[v] = tmp
     3002
     3003    sage_free(t_dist)
     3004    sage_free(dist)
     3005
     3006    if distances and paths:
     3007        return d_dist, d_prec
     3008    elif paths:
     3009        return d_prec
     3010    else:
     3011        return d_dist
     3012
    28543013cdef class Search_iterator:
    28553014    r"""
    28563015    An iterator for traversing a (di)graph.
  • sage/graphs/generic_graph.py

    diff -r 361a4ad7d52c -r 954991d687ee sage/graphs/generic_graph.py
    a b  
    1009210092        """
    1009310093        return self.shortest_path_length(u, v)
    1009410094
    10095     def distance_all_pairs(self):
     10095    def distance_all_pairs(self, algorithm = "BFS"):
    1009610096        r"""
    1009710097        Returns the distances between all pairs of vertices.
    1009810098
     10099        INPUT:
     10100
     10101            - ``"algorithm"`` (string) -- two algorithms are available
     10102
     10103                * ``algorithm = "BFS"`` in which case the distances are computed
     10104                  through `n-1` different breadth-first-search.
     10105
     10106                * ``algorithm = "Floyd-Warshall"``, in which case the
     10107                  Floyd-Warshall algorithm is used.
     10108
     10109                The default is ``algorithm = "BFS"``.
     10110
    1009910111        OUTPUT:
    1010010112
    1010110113        A doubly indexed dictionary
     
    1011510127            sage: all([g.distance(0,v) == distances[0][v] for v in g])
    1011610128            True
    1011710129        """
    10118 
    10119         from sage.rings.infinity import Infinity
    10120         distances = dict([(v, self.shortest_path_lengths(v)) for v in self])
    10121        
    10122         # setting the necessary +Infinity
    10123         cc = self.connected_components()
    10124         for cc1 in cc:
    10125             for cc2 in cc:
    10126                 if cc1 != cc2:
    10127                     for u in cc1:
    10128                         for v in cc2:
    10129                             distances[u][v] = Infinity
    10130 
    10131         return distances
     10130        if algorithm == "BFS":
     10131            from sage.rings.infinity import Infinity
     10132            distances = dict([(v, self.shortest_path_lengths(v)) for v in self])
     10133       
     10134            # setting the necessary +Infinity
     10135            cc = self.connected_components()
     10136            for cc1 in cc:
     10137                for cc2 in cc:
     10138                    if cc1 != cc2:
     10139                        for u in cc1:
     10140                            for v in cc2:
     10141                                distances[u][v] = Infinity
     10142
     10143            return distances
     10144
     10145        elif algorithm == "Floyd-Warshall":
     10146            from sage.graphs.base.c_graph import floyd_warshall
     10147            return floyd_warshall(self,paths = False, distances = True)
     10148
     10149        else:
     10150            raise ValueError("The algorithm keyword can be equal to either \"BFS\" or \"Floyd-Warshall\"")
    1013210151
    1013310152    def eccentricity(self, v=None, dist_dict=None, with_labels=False):
    1013410153        """
     
    1091610935                lengths[v] = len(paths[v]) - 1
    1091710936            return lengths
    1091810937
    10919     def shortest_path_all_pairs(self, by_weight=True, default_weight=1):
     10938    def shortest_path_all_pairs(self, by_weight=False, default_weight=1):
    1092010939        """
    1092110940        Uses the Floyd-Warshall algorithm to find a shortest weighted path
    1092210941        for each pair of vertices.
    1092310942       
    10924         The weights (labels) on the vertices can be anything that can be
    10925         compared and can be summed.
    10926        
    10927         INPUT:
    10928        
    10929        
    10930         -  ``by_weight`` - If False, figure distances by the
    10931            numbers of edges.
    10932        
    10933         -  ``default_weight`` - (defaults to 1) The default
    10934            weight to assign edges that don't have a weight (i.e., a label).
    10935        
    10936        
    10937         OUTPUT: A tuple (dist, pred). They are both dicts of dicts. The
    10938         first indicates the length dist[u][v] of the shortest weighted path
    10939         from u to v. The second is more complicated- it indicates the
    10940         predecessor pred[u][v] of v in the shortest path from u to v.
     10943        INPUT:
     10944       
     10945       
     10946        - ``by_weight`` - Whether to use the labels defined over the edges as
     10947           weights. If ``False`` (default), the distance between `u` and `v` is
     10948           the minimum number of edges of a path from `u` to `v`.
     10949       
     10950        - ``default_weight`` - (defaults to 1) The default weight to assign
     10951           edges that don't have a weight (i.e., a label).
     10952
     10953           Implies ``by_weight == True``.
     10954       
     10955        OUTPUT:
     10956
     10957            A tuple ``(dist, pred)``. They are both dicts of dicts. The first
     10958            indicates the length ``dist[u][v]`` of the shortest weighted path
     10959            from `u` to `v`. The second is a compact representation of all the
     10960            paths- it indicates the predecessor ``pred[u][v]`` of `v` in the
     10961            shortest path from `u` to `v`.
     10962
     10963        .. NOTE::
     10964
     10965            Two version of the Floyd-Warshall algorithm are implemented. One is
     10966            pure Python, the other one is Cython. The first is used whenever
     10967            ``by_weight = True``, and can perform its computations using the
     10968            labels on the edges for as long as they can be added together. The
     10969            second one, used when ``by_weight = False``, is much faster but only
     10970            deals with the topological distance (each edge is of weight 1, or
     10971            equivalently the length of a path is its number of edges). Besides,
     10972            the current version of this second implementation does not deal with
     10973            graphs larger than 65536 vertices (which already represents 16GB of
     10974            ram when running the Floyd-Warshall algorithm).
    1094110975       
    1094210976        EXAMPLES::
    1094310977       
    1094410978            sage: G = Graph( { 0: {1: 1}, 1: {2: 1}, 2: {3: 1}, 3: {4: 2}, 4: {0: 2} }, sparse=True )
    1094510979            sage: G.plot(edge_labels=True).show() # long time
    10946             sage: dist, pred = G.shortest_path_all_pairs()
     10980            sage: dist, pred = G.shortest_path_all_pairs(by_weight = True)
    1094710981            sage: dist
    1094810982            {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}}
    1094910983            sage: pred
     
    1095110985            sage: pred[0]
    1095210986            {0: None, 1: 0, 2: 1, 3: 2, 4: 0}
    1095310987       
    10954         So for example the shortest weighted path from 0 to 3 is obtained
    10955         as follows. The predecessor of 3 is pred[0][3] == 2, the
    10956         predecessor of 2 is pred[0][2] == 1, and the predecessor of 1 is
    10957         pred[0][1] == 0.
     10988        So for example the shortest weighted path from `0` to `3` is obtained as
     10989        follows. The predecessor of `3` is ``pred[0][3] == 2``, the predecessor
     10990        of `2` is ``pred[0][2] == 1``, and the predecessor of `1` is
     10991        ``pred[0][1] == 0``.
    1095810992       
    1095910993        ::
    1096010994       
    1096110995            sage: G = Graph( { 0: {1:None}, 1: {2:None}, 2: {3: 1}, 3: {4: 2}, 4: {0: 2} }, sparse=True )
    10962             sage: G.shortest_path_all_pairs(by_weight=False)
     10996            sage: G.shortest_path_all_pairs()
    1096310997            ({0: {0: 0, 1: 1, 2: 2, 3: 2, 4: 1},
    1096410998            1: {0: 1, 1: 0, 2: 1, 3: 2, 4: 2},
    1096510999            2: {0: 2, 1: 1, 2: 0, 3: 1, 4: 2},
     
    1097011004            2: {0: 1, 1: 2, 2: None, 3: 2, 4: 3},
    1097111005            3: {0: 4, 1: 2, 2: 3, 3: None, 4: 3},
    1097211006            4: {0: 4, 1: 0, 2: 3, 3: 4, 4: None}})
    10973             sage: G.shortest_path_all_pairs()
     11007            sage: G.shortest_path_all_pairs(by_weight = True)
    1097411008            ({0: {0: 0, 1: 1, 2: 2, 3: 3, 4: 2},
    1097511009            1: {0: 1, 1: 0, 2: 1, 3: 2, 4: 3},
    1097611010            2: {0: 2, 1: 1, 2: 0, 3: 1, 4: 3},
     
    1099311027            3: {0: 4, 1: 2, 2: 3, 3: None, 4: 3},
    1099411028            4: {0: 4, 1: 0, 2: 3, 3: 4, 4: None}})
    1099511029        """
     11030        if default_weight != 1:
     11031            by_weight = True
     11032
     11033        if by_weight is False:
     11034            from sage.graphs.base.c_graph import floyd_warshall
     11035            return floyd_warshall(self, distances = True)
     11036
    1099611037        from sage.rings.infinity import Infinity
    1099711038        dist = {}
    1099811039        pred = {}