# Sage Floyd algorithm in Cython

Hi

Classical Floyd algorithm computes all shortest path in a weighted graph and this algorithm has a $O(n^3)$ complexity (-> it can be slow). So i'm a bit surprised to see that the Cython version provided in Sage works only if by_weight==False meaning that all weights are 1 by default (its computes the transitive closure of the graph).

As a consequence, does it mean the only way to compute the "true" all-shortest-path in Sage is by using the pure Python implementation ?

edit retag close merge delete

Sort by » oldest newest most voted

I have searched more ;)

From Sagemath source code, we can spot the following that leaves no room for hope:

    Traceback (most recent call last):
...
ValueError: Algorithm 'Floyd-Warshall-Cython' does not work with weights.


But all is not lost. For each of Numpy (via the NetworkX library) and Scipy, there is an efficient implementation of Floyd algorithm (for Scipy, the implementation is as efficient as pure C code). The following code shows how to use them within SageMath:

import networkx as nx
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import floyd_warshall
import numpy as np

def wedges2matrix(wedges, n, loop=True):

M=[[float('inf')]*n for _ in range(n)]
for a, b, w in wedges:
M[a][b]=w
if loop:
for a in range(n):
M[a][a]=0
return M

def floyd_scipy(wedges, n):
M=wedges2matrix(wedges, n)
return floyd_warshall(np.array(M))

def floyd_np(wedges, n):
G = nx.DiGraph()
M=wedges2matrix(wedges, n)
dist=nx.floyd_warshall_numpy(G, nodelist=range(n))
return dist

n=10 #nb of vertices
wedges=[(9, 3, 38.0), (6, 0, 6.0), (0, 2, 53.0), (4, 7, 66.0),
(1, 9, 84.0), (8, 5, 24.0), (2, 0, 73.0), (4, 2, 79.0),
(6, 7, 78.0), (2, 4, 90.0), (7, 2, 81.0), (7, 0, 57.0),
(1, 3, 89.0), (0, 3, 17.0)]

A=floyd_np(wedges, n)
B=floyd_scipy(wedges, n)

print A
print "-------------"
print B


outputting identical matrices.

This is an opportunity to benchmark Sage implementation:

from time import clock

def dist2matrix(dist,n):
M=matrix(n, n, [[float('inf')]*n for _ in range(n)])
for a in dist:
for b in dist[a]:
M[a,b]=dist[a][b]
for i in range(n):
M[i,i]=0
return M

def random_edges(n, density, max_weight=100):
M=n*(n-1)//2
m=int(density*M)
edges=set()
wedges=[]
while len(edges)<m:
L=(randrange(n),randrange(n))
if L!=L and L not in edges:
w=float(randrange(max_weight))
wedges.append(L+(w,))
return wedges

def floyd(wedges, n):
M=wedges2matrix(wedges, n)
for k in range(n):
A=M[k]
for i in range(n):
B=M[i]
for j in range(n):
if B[k] + A[j]<B[j]:
B[j]=B[k] + A[j]
return M

def floyd_sage(wedges, n):
G = DiGraph(wedges)
dist, _ = G.shortest_path_all_pairs(by_weight = True, algorithm='Floyd-Warshall-Python')
return dist

n=300
density=0.33

wedges=random_edges(n, density)

start=clock()
M=floyd_scipy(wedges, n)
end=clock()
print("Scipy Floyd \t: %.2f" %(end-start))

start=clock()
M=floyd_np(wedges, n)
end=clock()
print("Numpy Floyd  \t: %.2f" %(end-start))

start=clock()
A=floyd(wedges, n)
end=clock()
print("Usual Floyd  \t: %.2f" %(end-start))

start=clock()
dist=floyd_sage(wedges, n)
end=clock()
print("Sage Floyd \t: %.2f" %(end-start))

matrix(RR, A)==dist2matrix(dist,n)


outputting on my machine:

Scipy Floyd     : 0.03
Numpy Floyd     : 0.11
Usual Floyd     : 2.52
Sage Floyd  : 10.43
True


It seems that Sage pure Python implementation is not very efficient...

Edit: final check was missing.

more

1

Are you welling to add algorithm='Floyd-Warshall-scipy', algorithm='Floyd-Warshall-numpy' to the method shortest_path_all_pairs in Sage and even improve the Sage pure python implementation? That would be awesome.

Thanks for the suggestion, I will reflect on it, but not much time at the moment. I have never done this before so could be time consuming to dive in this process. More over, the Floyd'algorithm Python implementation above needs some wrapping to match Sage interface, perpaps causing lower performance, this needs some testing before.

scipy is a dependency of Sage anyways, so there doesn't seem to be any reason for Sage to maintain its own inferior implementation of this algorithm when at most it could just wrap the one from scipy.