Ticket #10885: trac_10885.patch
File trac_10885.patch, 14.1 KB (added by , 12 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 2851 2851 else: 2852 2852 return True 2853 2853 2854 def 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 2854 3013 cdef class Search_iterator: 2855 3014 r""" 2856 3015 An iterator for traversing a (di)graph. -
sage/graphs/generic_graph.py
diff -r 361a4ad7d52c -r 954991d687ee sage/graphs/generic_graph.py
a b 10092 10092 """ 10093 10093 return self.shortest_path_length(u, v) 10094 10094 10095 def distance_all_pairs(self ):10095 def distance_all_pairs(self, algorithm = "BFS"): 10096 10096 r""" 10097 10097 Returns the distances between all pairs of vertices. 10098 10098 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 10099 10111 OUTPUT: 10100 10112 10101 10113 A doubly indexed dictionary … … 10115 10127 sage: all([g.distance(0,v) == distances[0][v] for v in g]) 10116 10128 True 10117 10129 """ 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\"") 10132 10151 10133 10152 def eccentricity(self, v=None, dist_dict=None, with_labels=False): 10134 10153 """ … … 10916 10935 lengths[v] = len(paths[v]) - 1 10917 10936 return lengths 10918 10937 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): 10920 10939 """ 10921 10940 Uses the Floyd-Warshall algorithm to find a shortest weighted path 10922 10941 for each pair of vertices. 10923 10942 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). 10941 10975 10942 10976 EXAMPLES:: 10943 10977 10944 10978 sage: G = Graph( { 0: {1: 1}, 1: {2: 1}, 2: {3: 1}, 3: {4: 2}, 4: {0: 2} }, sparse=True ) 10945 10979 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) 10947 10981 sage: dist 10948 10982 {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}} 10949 10983 sage: pred … … 10951 10985 sage: pred[0] 10952 10986 {0: None, 1: 0, 2: 1, 3: 2, 4: 0} 10953 10987 10954 So for example the shortest weighted path from 0 to 3 is obtained10955 as follows. The predecessor of 3 is pred[0][3] == 2, the10956 predecessor of 2 is pred[0][2] == 1, and the predecessor of 1is10957 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``. 10958 10992 10959 10993 :: 10960 10994 10961 10995 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() 10963 10997 ({0: {0: 0, 1: 1, 2: 2, 3: 2, 4: 1}, 10964 10998 1: {0: 1, 1: 0, 2: 1, 3: 2, 4: 2}, 10965 10999 2: {0: 2, 1: 1, 2: 0, 3: 1, 4: 2}, … … 10970 11004 2: {0: 1, 1: 2, 2: None, 3: 2, 4: 3}, 10971 11005 3: {0: 4, 1: 2, 2: 3, 3: None, 4: 3}, 10972 11006 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) 10974 11008 ({0: {0: 0, 1: 1, 2: 2, 3: 3, 4: 2}, 10975 11009 1: {0: 1, 1: 0, 2: 1, 3: 2, 4: 3}, 10976 11010 2: {0: 2, 1: 1, 2: 0, 3: 1, 4: 3}, … … 10993 11027 3: {0: 4, 1: 2, 2: 3, 3: None, 4: 3}, 10994 11028 4: {0: 4, 1: 0, 2: 3, 3: 4, 4: None}}) 10995 11029 """ 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 10996 11037 from sage.rings.infinity import Infinity 10997 11038 dist = {} 10998 11039 pred = {}