# hw18 - Pablo Blanco
# OK to post
read(`C18.txt`);



AveLengthSP:=proc(G,a) local Paths,P,totLen,k,mu,len,SS,i,A:
 Paths:=DijP(G,a)[2]:
 totLen:=0:
 k:=0:
 A:=Array([0$G[1]]):
 for P in Paths do:
  len:=nops(P)-1:
  totLen+=len:
  k++:
  A[k]:=len:
 od:

 mu:=totLen/k:
 SS:= add((A[i]-mu)^2,i=1..G[1]): # sum of squares
 
 mu,SS/k:
end:



AveLengthStat:=proc(n,p,K,M) local GraphAvgs, G,i,j,a,VtxAvgs,BigAvg,BigVar:
 GraphAvgs:=Array([infinity$M]):
 

 for i from 1 to M do:
  G:=RWG(n,p,K):
  VtxAvgs:=Array([infinity$n]):
  for a from 1 to n do:
   VtxAvgs[a]:=AveLengthSP(G,a)[1]:
  od:
  GraphAvgs[i]:= add(VtxAvgs)/n:
 od:

 BigAvg:=add(GraphAvgs)/M:
 BigVar:=add( (GraphAvgs[j]-BigAvg)^2  ,j=1..M)/M:

 BigAvg,BigVar:
end:

# Experiment:
# Ran this once: (took about 3.7 hrs)
# AveLengthStat(100,1/4,2,10000);
# 87964001/50000000, 318700761999/2500000000000000