# hw7 Pablo Blanco
# OK to post


# Note: some of this code is re-used from other projects and so uses some extra package(s).
#with(GraphTheory): #replaced Neighbors() with Neis from C3.txt.
with(ArrayTools):
randomize():


# Estimates the average number of leaves in a random tree on n vertices by running K trials.
#
# Experiment results. K=100: average was 37.27. 
# K=1000: 37.521
# Did not run very fast. Takes about 15 sec to count the leaves of 10 trees.
EstimateAverageNumberOfLeaves:=proc(n,K) local i,T,total:
 total:=0: # total number of leaves so far
 for i from 1 to K do:
  total:=total+NumberOfLeaves(RandomTree(n)):
 od:

 total/K:
end:


# counts the number of leaves in a labeled tree T. T=[n,E].
NumberOfLeaves:=proc(T) local n,E,Nhbrs,e,lfct,A:
 n:=T[1]:
 E:=T[2]:
 A:=Array([0$(n)]): #degree sequence
 Nhbrs:=Neis([n,E]):        ### use neighbors
 lfct:=n: # leaf counter. Assume every vertex is a leaf and then decrease accordingly.

 for e in E do:
  if A[e[1]] = 1 then:
   lfct--:
  fi:
  if A[e[2]] = 1 then:
   lfct--:
  fi:

  A[e[1]]++:
  A[e[2]]++:
 od:
 lfct:
end:

# constructs a uniformly random spanning tree, returns its set of EDGES
RandomTree:=proc(n) local E,Next,i,j,N,snt,TrEd:
 E:=RandomGraph(n,1):
 Next:=RandomTreeWithRoot(n,E,1):
 
 if type(Next, Array) then 
  N:=NumElems(Next):    
 else:
   N:=nops(Next):
 fi: 
 
 TrEd:= Array([0$(N-1)]):  #
 snt:=false:
 j:=0:
 for i from 1 to N do: 
  if snt or Next[i]<>0 then:
   j++:
   TrEd[j]:= {i,Next[i]}:
  else:
   snt:=true:
  fi:

 od:

[n,{seq(TrEd[i],i=1..N-1)}]:
end:

# generates random N-vertex graph. returns N and edge-set. Takes in an optional second argument, p (probability) which is assumed to be rational.
RandomGraph:=proc(N,p:=1/2) local E,i,j:
E:={}:

for i from 1 to N-1 do:
 for j from i+1 to N do:
   if rand(1..denom(p))()<=numer(p) then
    E:= E union {{i,j}}:
   fi:
 od:
od:


E
end:

#Given a connected graph n,E, constructs a random spanning tree
RandomTreeWithRoot:=proc(n,E,r) local i,InTree, Next,Nhbrs,u,tCount:
 InTree:=Array([false$n]): #
 InTree[r]:=true:
 Next:=Array([0$n]): 
 Nhbrs:=Neis([n,E]):  # use neighbors
 
 for i from 1 to n do:
  u:=i:
  while not InTree[u] do:
    Next[u]:=RandomSuccessor(Nhbrs[u]):
    u:=Next[u]:
  od:
  u:=i:
  while not InTree[u] do:
   InTree[u]:=true:
   u:=Next[u]:
  od:
 od:
 
 Next:
end:

# picks a random successor, given a vertex u and its neighborhood N, picks a random neighbor
RandomSuccessor:=proc(N):
 
 if type(N, Array) then                       #
  RETURN(N[rand(1..NumElems(N))()]):         #  
 fi:                                         #
 N[rand(1..nops(N))()]:
end:


### from C3.txt
#Neis(G): inputs a graph G=[n,E] and outputs a list of length n, N, such that
#N[i] is the set of neighbors of vertex i
Neis:=proc(G) local n,E,N,i,e:
 n:=G[1]: 
 E:=G[2]:
 for i from 1 to n do
  N[i]:={}:
 od:
 for e in E do
  N[e[1]]:=N[e[1]] union {e[2]}:
  N[e[2]]:=N[e[2]] union {e[1]}:
 od:
 [seq(N[i],i=1..n)]:
end: