#OK to post hw
#Michael Yen, 11/21/20, Assignment 21

#Lecture 21: Due Nov. 22, 9:00pm. Email ShaloshBEKhad@gmail.com an attachment

#hw21FirstLast.txt

#Indicate whether it is OK to post

#1.How many labeled connected graphs with 50 vertices and 70 edges are?
f := taylor(log(add((1 + a)^(n*(n - 1)/2)*x^n/n!, n = 0 .. 51)), x = 0, 51);
coeff(expand(coeff(f, x, 50)*50!), a, 70);
#557570125262936299552095051452832858951958306952392379552719516748979534894703561052813706583449695070191616000000

#2. How many labeled graphs with i connected components with 50 vertices and 70 edges are for i=2,3,4,5? 
#Explain how you do it, by using the taylor command in Maple.
seq(coeff(expand(coeff(taylor(log(add((1 + a)^(n*(n - 1)/2)*x^n/n!, n = 0 .. 51))^i/i!, x = 0, 51), x, 50)*50!), a, 70),i=2..5);

#The EGF of the molecule is the EGF of the atoms^k/k!.
#So, I just used the same answer as question 1 but made sure to do f^k/k!.

#i=2:
#2267762282072205271778791270680596849921516338574571823883456519281730035101965414722739023408444177334403072000000,
#i=3:
#4084614615513849564059280183743190588234377301848961646892154560205839868330317696952259164363869748482983526400000,
#i=4:
#4336520950646656327224993268880264842114984745013581637572959168612209444462829083123209142308423023328834795520000, 
#i=5:
#3046643850794698885807347271650949921058341236747610399052555568087884858386903803880100124793197196674954096998400

#3. In the homework assigment for Lecture 18 you were asked to write a procedure
#EstimateCutOff(n,M)

#by simulating it M times, and using as cutoff 1.05. Procedure FindCutOff(n,co) in

#M21.txt

#does it EXACTLY, for an arbitrary parameter co (foremerly I asked you to take 1.05).

#Comparee the output of the simulated procedure with the exact one, i.e. compare your EstimateCutOff(30,1000) 
#with my exact FindCutOff(30,1.05)

#Here is my EstimateCutOff from Hw18
EstimateCutOff:=proc(n,M) local k, ave, beforeAve, start, fin:
start:=0:
fin:=binomial(30,2):
k:=ceil((start+fin)/2):
ave:=AveNuCC(n,k,M):
beforeAve:=AveNuCC(n,k-1,M):
while ave>=1.05 or beforeAve<1.05 and start < fin do:
 if ave>=1.05 then
  start:=k:
  k:=ceil((start+fin)/2):
 else 
  fin:=k:
  k:=ceil((start+fin)/2):
 fi:
 ave:=AveNuCC(n,k,M):
 beforeAve:=AveNuCC(n,k-1,M):
od:
k:
end:

#EstimateCutOff(30,1000) = 82, taken from hw18's answer

#FindCutOff(30,1.05) = 84

#I would say they are pretty close, my estimation does a good job of guessing the actual cutoff.

#4. Try to do the optional challenge from the homework of Lecture 18 exactly, using evalf(FindCutOff(n,1.05)/
#n), for n= 20,30,40,50,60. Do you see a trend?

#Here is the optional challenge from hw18.
#[Optional challenge] Running it for various n from 100 to 200, and M=1000 (or possibly larger) Can you 
#conjecture an expression for that cut-off in terms of n?
evalf(FindCutOff(n,1.05)/n), for n= 20,30,40,50,60.

#20: 2.45
#30: 2.8
#40: 3.
#50: 3.16
#60: 3.3

#The trend is that the cutoff increases as n increases, but the increase gets continuously smaller.
#I don't see an exact expression to predict the cut-off in terms of n.

#5. Generalize procedure
#EstProbConn(n,m,N):

#and write a procedure

#EstProbKcomps(n,m,k,N):

#that estimates the probability that a random graph with n vertices and m edges has exactly k components, by 
#sampling N such random graphs.

#Note that EstProbKcomps(n,m,1,N) is exactly the same as EstProbConn(n,m,N) (except, since we are doing 
#simulations, you get different answers each time)

EstProbKcomps:=proc(n,m,k,N) local i,co, G:

co:=0:

for i from 1 to N do

G:=RandGr(n,m):

if nops(CCs(G))=k then
 co:=co+1:
fi:

od:

evalf(co/N):
end:

#6. Generalize procedure
#ProbConn(n,m):

#To write a procedure

#ProbKcomponents(n,m,k)

#that uses the method of exponential generating functions to find the EXACT probability that a random graph #with n vertices and m edges has k components. Note that ProbKcomponents(n,m,1) is EXACTLY the same as #ProbConn(m,m,1)

#[Hints: (i) you may need to write another procedure generalizing one of those that I wrote and that was used 
#in ProbConn(n,m). (ii) the egf of labeled graphs with n vertices and k components with keeping track of the 
#number of edges (i.e. the weight-enumerator according to the weight anumber of edges) is

#(log(Sum((1+a)^(n*(n-1)/2)*x^n/n!,n=0..infinity))^k/k!

#How does ProbKcomponents(40,50,2) differ from EstProbKcomps(40,50,2,300)?

WtEdConGclever2:=proc(N,a,k) local f,i,x:
option remember:
f:=log(add((1+a)^binomial(i,2)*x^i/i!,i=0..N+1))^k/k!: #just modified this line right here
f:=taylor(f,x=0,N+1):
[seq(expand(coeff(f,x,i)*i!),i=1..N)]:
end:

ProbKcomponents:=proc(n,m,k)  local a:
evalf(coeff(WtEdConGclever2(n,a,k)[n],a,m)/binomial(n*(n-1)/2,m)):
end:

ProbKcomponents(40,50,2)
#0.09801230118
EstProbKcomps(40,50,2,300)
#0.1066666667
%-%%
#0.00865436552
#0.00865436552/0.09801230118*100 = 8.829876879%, the percent error is not that bad but it could be better.