# 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[0].intersection(u[1])
if not T.is_empty():
# sink is found
t = P.subposet(T).minimal_elements()[0]
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[0]+fgh[2] for fgh in cFGH)
G = sum((fgh[0]+fgh[1])*H/(fgh[0]+fgh[2]) for fgh in cFGH)
F = x*prod(fgh[0]+fgh[2] + (fgh[0]+fgh[1])/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[0],t)
# N represents L at this point
L_FGH = FGH(N, v, true_leaves)
N.add_edge(p[0],t)
N.delete_edge(p[1],t)
# N represents R at this point
R_FGH = FGH(N, v, true_leaves)
N.delete_edge(p[0],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[0],u)
u = w[0]
# N represents B at this point
B_FGH = [ FGH(N, u, true_leaves)
for u in Set(P.closed_interval(v,p[0])) + Set(P.closed_interval(v,p[1])) ]
B_prodH = prod(fgh[2] for fgh in B_FGH)
# apply Theorem 9
H = L_FGH[2] + R_FGH[2] - Nst_FGH[2] * (t_FGH[0] + t_FGH[2])
G = L_FGH[1] + R_FGH[1] - Nst_FGH[1] * (t_FGH[0] + t_FGH[2])
- (t_FGH[0] + t_FGH[1]) * B_prodH
F = L_FGH[0] + R_FGH[0] - Nst_FGH[0] * (t_FGH[0] + t_FGH[2])
- (t_FGH[0] + t_FGH[1]) * (prod((fgh[0]+fgh[1])/x + fgh[2] for fgh in B_FGH) - B_prodH)
return [F,G,H]

### Like this:

Like Loading...