# 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]

Advertisements