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