# Sage code for Cactus network scoring

```# Function FGH( N, v, true_leaves )
# Input:
#   N is a cactus network, possibly a subnetwork of the initial network
#   v is a vertex in N
#   true_leaves is an optional parameter, used internally,
#               set of vertices that are leaves in the initial network;
#               should be omitted when N is an initial network
# Output:
#   [ F_v(x), G_v(x), H_v(x) ]

from itertools import combinations

def FGH( N, v, true_leaves = None ):

# test if N is a DAG
if not N.is_directed_acyclic():
raise ValueError, "Input graph is not a network!"

if true_leaves is None:
# N is a parent network, so we initialize true_leaves:
true_leaves = Set(N.sinks())

# define g.f. variable x
x = QQ['x'].0

if N.out_degree(v) == 0:
# v is a leaf
if v in true_leaves:
return [x,0,0]
else:
return [0,0,1]

P = Poset(N)

# look for a sink in N_v if it exists
t = v
S = [ Set(P.principal_upper_set(u)) for u in N.neighbors_out(v) ]
for u in combinations(S,2):
T = u.intersection(u)
if not T.is_empty():
# sink is found
t = P.subposet(T).minimal_elements()
break

if t == v:
# v is a regular vertex
# compute FGH at the children of v
cFGH = [ FGH(N,u,true_leaves) for u in N.neighbors_out(v) ]
# apply Theorem 3
H = prod(fgh+fgh for fgh in cFGH)
G = sum((fgh+fgh)*H/(fgh+fgh) for fgh in cFGH)
F = x*prod(fgh+fgh + (fgh+fgh)/x for fgh in cFGH) - x*H - G;

else:
# v is a source and t is a sink
t_FGH = FGH(N, t, true_leaves)

p = N.neighbors_in(t)
if len(p) != 2:
raise ValueError, "N is not a cactus network"

N.delete_edge(p,t)
# N represents L at this point
L_FGH = FGH(N, v, true_leaves)

N.delete_edge(p,t)
# N represents R at this point
R_FGH = FGH(N, v, true_leaves)

N.delete_edge(p,t)
# N represents N_{s\t} U N_t at this point
Nst_FGH = FGH(N, v, true_leaves)

# remove edges in the branching path v => t from N to form B
for q in p:
u = q
while u!=v:
w = N.neighbors_in(u)
if len(w)!=1:
raise ValueError, "Parents error w"
N.delete_edge(w,u)
u = w
# N represents B at this point
B_FGH = [ FGH(N, u, true_leaves)
for u in Set(P.closed_interval(v,p)) + Set(P.closed_interval(v,p)) ]
B_prodH = prod(fgh for fgh in B_FGH)
# apply Theorem 9
H = L_FGH + R_FGH - Nst_FGH * (t_FGH + t_FGH)
G = L_FGH + R_FGH - Nst_FGH * (t_FGH + t_FGH)
- (t_FGH + t_FGH) * B_prodH
F = L_FGH + R_FGH - Nst_FGH * (t_FGH + t_FGH)
- (t_FGH + t_FGH) * (prod((fgh+fgh)/x + fgh for fgh in B_FGH) - B_prodH)

return [F,G,H]
```