# OK to post homework
# Aurora Hiveley, 3/31/25, Assignment 18

Help := proc(): print(`AveLengthSP(G,a), AveLengthStat(n,p,K,M)`): end:


# AveLengthSP(G,a): inputs a weighted graph and a starting vertex a and uses DijP(G,a)[2] to compute the average length of a shortest path from a.

AveLengthSP := proc(G,a) local p,P:
P := DijP(G,a)[2]:
add(nops(p), p in P)/nops(P):
end:



# AveLengthStat(n,p,K,M): picks a random weighted graph according to RWG(n,p,K) M times, finds for each AveLengthSP(G,a) 
# and finds the average and variance for each run.

AveLengthStat := proc(n,p,K,M) local L,L1,G,i,avg,var,a,l:
L := []:

for i from 1 to M do 
    G := RWG(n,p,K):
    L1 := []:

    for a from 1 to n do 
        L1 := [op(L1), AveLengthSP(G,a)]:
    od:

    L := [op(L), add(l,l in L1)/n]:
od:

avg := add(l,l in L)/nops(L):
var := sqrt(add((l-avg)^2, l in L)/nops(L)):

evalf([avg,var]):
end:


# Run AveLengthStat(100,1/4,2,10000) five times. What did you get? Are they similar?
# time(AveLengthStat(100,1/4,2,10)) was 11.000
# time(AveLengthStat(100,1/4,2,100)) was 120.828
# by my math, time(AveLengthStat(100,1/4,2,10000)) ~ 10000 seconds, or ~ 2.78 hours, 
# so i'm only going to run it with M = 1000, which took 40 mins for all five trials together 

# [2.759949200, 0.01146021637]
# [2.760042500, 0.01145728169]
# [2.759840700, 0.01145451149]
# [2.759119400, 0.01151063698]
# [2.758823900, 0.01111790173]






















### copied from C18.txt 
#C18.txt: March 31, 2025 Graph algithms
Help18:=proc(): print(` RWG(n,p,K), Dij(G,a) , DijP(G,a) `): end:

#RWG(n,p,K): inputs a pos. integer n and outputs a triple [n,E,DT]
#where the vertices are {1,...,n} the edges are of the form {i,j} and
#DT is table such DT[{i,j}] is the weight of the edge {i,j}
RWG:=proc(n,p,K) local G,DT,e,ra,E:
ra:=rand(1..K):

G:=RG(n,p):
E:=G[2]:
for e in E do
DT[e]:=ra():
od:
[n,E,op(DT)]:
end:


#Dij(G,a): inputs a weighted graph G=[n,E,DT] and outputs and a starting 
#vertex a outputs a list L of length n such that L[i] is the MINIMAL DISTANCE
#min. distance
#from a to i. In particular L[a]=0
Dij:=proc(G,a) local n,E,DT,Nei,e,i,Uvis,Tt,T,Vis,cu,Nei1,v,cha,rec,i1:
n:=G[1]:
E:=G[2]:
DT:=G[3]:
#N is a table of length n such that N[i] is the set of VERTICES j such that {i,j} belongs to E
for i from 1 to n do
 Nei[i]:={}:
od:
for e in E do
 Nei[e[1]]:=Nei[e[1]] union {e[2]}:
 Nei[e[2]]:=Nei[e[2]] union {e[1]}:
od:
 
#Uvis is he currently set of still unvisoted vertices
Uvis:={seq(i,i=1..n)}:
#T[i] is the permanent label (the min. distance from a to i
#Tt[i]: is the current best upper bound
for i from 1 to n do
 Tt[i]:=infinity:
od:


T[a]:=0:
Uvis:=Uvis minus {a}:
Vis:=[a]:

while Uvis<>{} do
cu:=Vis[nops(Vis)]:
Nei1:=Nei[cu]:

for v in Nei1 do
if T[cu]+DT[{cu,v}]<Tt[v] then
 Tt[v]:=T[cu]+DT[{cu,v}]:
fi:
od:

cha:=Uvis[1]:
rec:=Tt[cha]:
for i1 from 2 to nops(Uvis) do
  if Tt[Uvis[i1]]<rec then
     cha:=Uvis[i1]:
     rec:=Tt[cha]:
  fi:
od:
T[cha]:=Tt[cha]:
Vis:=[op(Vis), cha]:

Uvis:=Uvis minus {cha}:
od:
[seq(T[i],i=1..n)]:
end:

#DijP(G,a): inputs a weighted graph G=[n,E,DT] and outputs and a starting 
#vertex a outputs a list L of length n such that L[i] is the MINIMAL DISTANCE
#it also returns a list P such that P[i] is the path from a to i that realizes the
#min. distance
#from a to i. In particular L[a]=0
DijP:=proc(G,a) local n,E,DT,Nei,e,i,Uvis,Tt,T,Vis,cu,Nei1,v,cha,rec,i1,Paths, TempPaths:
n:=G[1]:
E:=G[2]:
DT:=G[3]:
#N is a table of length n such that N[i] is the set of VERTICES j such that {i,j} belongs to E
for i from 1 to n do
 Nei[i]:={}:
od:
for e in E do
 Nei[e[1]]:=Nei[e[1]] union {e[2]}:
 Nei[e[2]]:=Nei[e[2]] union {e[1]}:
od:
 
#Uvis is he currently set of still unvisoted vertices
Uvis:={seq(i,i=1..n)}:
#T[i] is the permanent label (the min. distance from a to i
#Tt[i]: is the current best upper bound
for i from 1 to n do
 Tt[i]:=infinity:
od:

Paths[a]:=[a]:

T[a]:=0:
Uvis:=Uvis minus {a}:
Vis:=[a]:

while Uvis<>{} do
cu:=Vis[nops(Vis)]:
Nei1:=Nei[cu]:

for v in Nei1 do
if T[cu]+DT[{cu,v}]<Tt[v] then
 Tt[v]:=T[cu]+DT[{cu,v}]:
 TempPaths[v]:=[op(Paths[cu]),v]:
fi:
od:

cha:=Uvis[1]:
rec:=Tt[cha]:
for i1 from 2 to nops(Uvis) do
  if Tt[Uvis[i1]]<rec then
     cha:=Uvis[i1]:
     rec:=Tt[cha]:
  fi:
od:
T[cha]:=Tt[cha]:
Vis:=[op(Vis), cha]:
Paths[cha]:=TempPaths[cha]:

Uvis:=Uvis minus {cha}:
od:
[seq(T[i],i=1..n)], [seq(Paths[i],i=1..n)]:
end:

#C2.txt: Jan. 27, 2025
Help2:=proc(): print(`LC(p), RG(n,p), `):end:
with(combinat):
#LC(p): inputs a rational number between 0 and 1 and outputs true with prob. p
LC:=proc(p) local a,b,ra:
if not (type(p,fraction) and p>=0 and p<=1) then
 RETURN(FAIL):
fi:
a:=numer(p):
b:=denom(p):
ra:=rand(1..b)():
if ra<=a then
 true:
else
 false:
fi:
end:

RG:=proc(n,p) local E,i,j:
E:={}:
for i from 1 to n do
 for j from i+1 to n do
   if LC(p) then
    E:=E union {{i,j}}:
   fi:
 od:
od:
[n,E]:
end:
####end from C2.txt