# Linear systems of graphs in Sage

Let $\Gamma$ be a graph. A divisor on $\Gamma$ is an element of the free group generated by the vertices $V$, ${\mathbb{Z}}[V]$.

We say that divisors $D$ and $D^\prime$ are linearly equivalent and write $D \sim D^\prime$ if $D^\prime-D$ is a principal divisor, i.e., if $D^\prime = D + \text{div}(f)$ for some function $f : V \rightarrow {\mathbb{Z}}$. Note that if $D$ and $D^\prime$ are linearly equivalent, they must have the same degree, since the degree of every principal divisor is 0. Divisors of degree 0 are linearly equivalent if and only if they determine the same element of the Jacobian. If $D$ is a divisor of degree 0, we denote by $[D]$ the element of the Jacobian determined by $D$. A divisor $D$ is said to be effective if $D(v) \geq 0$ for all vertices $v$. We write $D \geq 0$ to mean that $D$ is effective. The linear system associated to a divisor $D$ is the set

$|D| = \{ D^\prime \in \text{Div}(\Gamma ) : D^\prime \geq 0 \text{ and } D^\prime \sim D\},$

i.e., $|D|$ is the set of all effective divisors linearly equivalent to $D$. Note that if $D_1 \sim D_2$, then $|D_1| = |D_2|$. We note also that if $\text{deg}(D)<0$, then $|D|$ must be empty.

Sage can be used to compute the linear system of any divisor on a graph.

def linear_system(D, Gamma):
"""
Returns linear system attached to the divisor D.

EXAMPLES:
sage: Gamma2 = graphs.CubeGraph(2)
sage: Gamma1 = Gamma2.subgraph(vertices = ['00', '01'], edges = [('00', '01')])
sage: f = [['00', '01', '10', '11'], ['00', '01', '00', '01']]
sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f)
True
sage: PhiV = matrix_of_graph_morphism_vertices(Gamma1, Gamma2, f); PhiV
[1 0 1 0]
[0 1 0 1]
sage: D = vector([1,0,0,1])
sage: PhiV*D
(1, 1)
sage: linear_system(PhiV*D, Gamma1)
[(2, 0), (1, 1), (0, 2)]
sage: linear_system(D, Gamma2)
[(0, 2, 0, 0), (0, 0, 2, 0), (1, 0, 0, 1)]
sage: [PhiV*x for x in linear_system(D, Gamma2)]
[(0, 2), (2, 0), (1, 1)]

"""
Q = Gamma.laplacian_matrix()
CS = Q.column_space()
N = len(D.list())
d = sum(D.list())
#print d
lin_sys = []
if d < 0:
return lin_sys
if (d == 0) and (D in CS):
lin_sys = [CS(0)]
return lin_sys
elif (d == 0):
return lin_sys
S = IntegerModRing(d+1)^N
V = QQ^N
for v in S:
v = V(v)
#print D-v,v,D
if D-v in CS:
lin_sys.append(v)
return lin_sys


# Floyd-Warshall-Roy

The Floyd-Warshall-Roy algorithm is an algorithm for finding shortest paths in a weighted, directed graph. It allows for negative edge weights and detects a negative weight cycle if one exists. Assuming that there are no negative weight cycles, a single execution of the FWR algorithm will find the shortest paths between all pairs of vertices.

This algorithm is an example of dynamic programming, which allows one to break the computation down to simpler steps using some sort of recursive procedure. If $n=|V|$ then this is a $O(n^3)$-time algorithm. For comparison, the Bellman-Ford algorithm has complexity $O(|V|\cdot |E|)$, which is $O(n^3)$-time for “dense” graphs. However, Bellman-Ford only yields the shortest paths eminating from a single vertex. To achieve comparable output, we would need to iterate Bellman-Ford over all vertices, which would be an $O(n^4)$-time algorithm for “dense” graphs. Except for “sparse” graphs, Floyd-Warshall-Roy is better than an interated implementation of Bellman-Ford.

Here is a Sage implementation

def floyd_warshall_roy(A):
"""
Shortest paths

INPUT:

OUTPUT
dist  - a matrix of distances of shortest paths
paths - a matrix of shortest paths

EXAMPLES:
sage: A = matrix([[0,1,2,3],[0,0,2,1],[20,10,0,3],[11,12,13,0]]); A
sage: floyd_warshall_roy(A)
([[0, 1, 2, 2], [12, 0, 2, 1], [14, 10, 0, 3], [11, 12, 13, 0]],
[-1  1  2  1]
[ 3 -1  2  3]
[ 3 -1 -1  3]
[-1 -1 -1 -1])
sage: A = matrix([[0,1,2,4],[0,0,2,1],[0,0,0,5],[0,0,0,0]])
sage: floyd_warshall_roy(A)
([[0, 1, 2, 2], [+Infinity, 0, 2, 1], [+Infinity, +Infinity, 0, 5],
[+Infinity, +Infinity, +Infinity, 0]],
[-1  1  2  1]
[-1 -1  2  3]
[-1 -1 -1  3]
[-1 -1 -1 -1])
sage: A = matrix([[0,1,2,3],[0,0,2,1],[-5,0,0,3],[1,0,1,0]]); A
sage: floyd_warshall_roy(A)
Traceback (click to the left of this block for traceback)
...
ValueError: A negative edge weight cycle exists.
sage: A = matrix([[0,1,2,3],[0,0,2,1],[-1/2,0,0,3],[1,0,1,0]]); A
sage: floyd_warshall_roy(A)
([[0, 1, 2, 2], [3/2, 0, 2, 1], [-1/2, 1/2, 0, 3/2], [1/2, 3/2, 1, 0]],
[-1  1  2  1]
[ 2 -1  2  3]
[-1  0 -1  1]
[ 2  2 -1 -1])
"""
V = G.vertices()
E = [(e[0],e[1]) for e in G.edges()]
n = len(V)
dist = [[0]*n for i in range(n)]
paths = [[-1]*n for i in range(n)]
# initialization step
for i in range(n):
for j in range(n):
if (i,j) in E:
paths[i][j] = j
if i == j:
dist[i][j] = 0
elif A[i][j] != 0:
dist[i][j] = A[i][j]
else:
dist[i][j] = infinity
# iteratively finding the shortest path
for j in range(n):
for i in range(n):
if i != j:
for k in range(n):
if k != j:
if dist[i][k]>dist[i][j]+dist[j][k]:
paths[i][k] = V[j]
dist[i][k] = min(dist[i][k], dist[i][j] +dist[j][k])
for i in range(n):
if dist[i][i] < 0:
raise ValueError, "A negative edge weight cycle exists."
return dist, matrix(paths)