# OK to post homework
# Aurora Hiveley, 4/4/25, Assignment 19

Help := proc(): print(`EstimateAverageMinCost(n,p,K,N), EstimateTable(n,N,a,b)`): end:
# scroll down further for a separate Help() for Eulerian Paths!


# EstimateAverageMinCost(n,p,K,N): runs MST(RWG(n,p,K)) N times (where N is a large pos. integer) 
# and outputs the fraction of those that are connected (do not return FAIL), and among them 
# the average cost of a minimal spanning trees (both in decimals)

EstimateAverageMinCost := proc(n,p,K,N) local i,T,w,totN,totW :
totN := 0:
totW := 0:

for i from 1 to N do 
    T := MST(RWG(n,p,K)):
    
    if T<>FAIL then
        totN ++: 
        w := T[2]: 
        totW := totW + w:
    fi: 
od:

[evalf(totN/N), evalf(totW/N)]: # average weight across spanning trees that actually exist
end:



# EstimateTable(n,N,a,b): inputs a pos. integer n (not too large), a large (at least 1000) pos. integer N, 
# pos. integer a and pos. integer b and outputs the list L[i][K] (i between 1 and a) (K between 1 and b)
# such that L[i][K] is the estimated average of a minimal spanning tree of a random weighted spanning tree 
# with prob. of picking an edge with prob. i/a and picking a weight of an edge from 1 to K.

EstimateTable := proc(n,N,a,b) local i,K,L:

for i from 1 to a do
    for K from 1 to b do 
        L[i][K] := EstimateAverageMinCost(n,i/a,K,N):
    od:
od:

op(L):
end:


# What did you get for EstimateTable(10,3000,5,10)?
# table([1 = table([1 = [0.2100000000, 1.890000000], 2 = [0.2143333333, 2.655000000], 3 = [0.2196666667, 3.527000000], 4 = [0.2160000000, 4.272666667], 5 = [0.2190000000, 5.116000000], 6 = [0.2003333333, 5.458000000], 7 = [0.2043333333, 6.445333333], 9 = [0.2153333333, 8.250666667], 8 = [0.2050000000, 7.153666667], 10 = [0.2140000000, 8.863666667]]), 2 = table([1 = [0.8946666667, 8.052000000], 2 = [0.8980000000, 9.526333333], 3 = [0.8913333333, 11.59133333], 4 = [0.8993333333, 13.94500000], 5 = [0.9000000000, 16.39300000], 6 = [0.9063333333, 18.87266667], 7 = [0.9006666667, 21.03900000], 9 = [0.9020000000, 25.81666667], 8 = [0.8930000000, 23.18466667], 10 = [0.9033333333, 28.13766667]]), 3 = table([1 = [0.9983333333, 8.985000000], 2 = [0.9976666667, 9.464333333], 3 = [0.9973333333, 10.76400000], 4 = [0.9976666667, 12.44166667], 5 = [0.9983333333, 14.29233333], 6 = [0.9963333333, 15.89966667], 7 = [0.9960000000, 17.64900000], 9 = [0.9976666667, 21.39233333], 8 = [0.9953333333, 19.49500000], 10 = [0.9976666667, 23.40400000]]), 4 = table([1 = [1., 9.], 2 = [1., 9.108000000], 3 = [1., 9.763333333], 4 = [1., 10.83533333], 5 = [1., 12.07733333], 6 = [1., 13.32600000], 7 = [1., 14.75366667], 9 = [1., 17.31666667], 8 = [1., 16.05633333], 10 = [1., 18.67766667]]), 5 = table([1 = [0., 0.], 2 = [0., 0.], 3 = [0., 0.], 4 = [0., 0.], 5 = [0., 0.], 6 = [0., 0.], 7 = [0., 0.], 9 = [0., 0.], 8 = [0., 0.], 10 = [0., 0.]])])

# Run it a few times and see whether the results are close (of course they shouldn't be the same, this is simulation)
# results are generally pretty close, although obviously not exact due to the randomness component 
# but 3000 trials seems sufficiently large to give an accurate overview! 





### EULERIAN PATH CODE 

HelpEP := proc(): print(`EulerianPath(G), ET1(n,G,v), Neis(G)`): end:

## example call: 
# G := [5,{{1,2},{2,3},{3,4},{4,5},{5,1}}]:
# EulerianPath(G):

# G := [5,{{1,3},{1,4},{2,4},{2,5},{3,5}}]:
# EulerianPath(G):

# G := [5,{{1,2},{1,3},{1,4},{1,5},{2,3},{2,4},{2,5},{3,4},{3,5},{4,5}}]:
# EulerianPath(G):


# EulerianPath(G): inputs a Eulerian graph G=[n,E] (see section 6 of Robin Wilson's graph theory book), 
# that (i) checks that it is indeed Eulerian (ii) outputs a Eulerian trail using either Fleuri's algorithm 
# (see there) or any other correct method.

EulerianPath := proc(G) local n,N,i,u,v,w,T,EA,EV:
n := G[1]:

# check that G is Eulerian 
# G Eulerian iff every vertex has even degree 

N := Neis(G):

for i from 1 to n do 
    if nops(N[i]) mod 2<>0 then 
        RETURN(FAIL):
    fi:
od:


# construct Eulerian trail 
EA := G[2]:   # available edges for trail 
EV := {seq(1..n)}:      # vertices in G 

# can choose any starting vertex for trail, without loss i'll start with vertex 1 
# could also do v := rand(1..n)(): for some variety
v := 1: # starting vertex 
T := [v]:  # trail

while EA<>{} do 
    u := ET1(n,[EV, EA], v):    # next vertex in trail 

    if u=0 then         # no next vertex was selected in the trail, so something went wrong 
        RETURN(FAIL):
    fi:

    T := [op(T),u]:
    EA := EA minus {{u,v}}:
    
    if CC([n,EA],u)={} then     # u is isolated after {u,v} removed 
        EV := EV minus {u}: 
    fi:

    v := u:
od:

T:
end:



# selects next vertex in eulerian trail for G = [V,E] with current vertex v

ET1 := proc(n,G,v) local E,E1,N,u,w,i,C1,C2:
N := Neis(G)[v]: # current neighbors of v 
E := G[2]:

u := 0:     # placeholder vertex, will become next vertex in trail 
i := 1:     # starting index for neighbors of u

while u=0 and i <= nops(N) do 
    w := N[i]:  # candidate next vertex in trail

    E1 := E minus {{w,v}}:  # edge set if {u,w} is selected as next edge and hence {u,w} is removed

    ## CC(G,a) = set of vertices connected to a in G = [n,E]
    C1 := nops( CC([n,E],v) ):     # number of vert connected to v before {v,w} removed
    C2 := nops( CC([n,E1],v) ):    # number of vert connected to v after {v,w} removed
        
    if C1 = C2 then     # removal does not disconnect graph, so {u,w} not a bridge 
        u := w:
    fi:

    i++:
od:

# exits while loop if all edges {u,v} are bridges. then can pick any to be the next edge in trail 
if u=0 then 
    u := N[1]:
fi:

u:
end:



# Neis(G): the table of neighbors of the vertices in G = [V,E]

Neis := proc(G) local v,V,e,E,N:
V := G[1]:
E := G[2]:

for v in V do
    N[v]:={}:
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:

op(N):
end:









### copied from C19.txt 

#C19.txt April 3, 2025
Help19:=proc(): print(`CC(G,a), MST1(G, partT, AveE) , MST(G) `):end:

#MST1(G,partT,AvE): inputs a weighted graph G and performs ONE step in the Kruskal
#algoritm by looking at the cheapest member of AvE and trying to join it to
#partT if it does NOT create a cycle, either way we remove it from AvE
MST1:=proc(G,partT,AvE) local n, E,DT,e,AvE1,AvE1d,i1:
n:=G[1]:
E:=G[2]:
DT:=G[3]:

if partT minus E<>{} or AvE minus E<>{} then
 RETURN(FAIL):
fi:
AvE1:=convert(AvE,list):
AvE1d:=[seq(DT[AvE1[i1]],i1=1..nops(AvE1))]:
#e is the cheapest edge
e:=AvE1[min[index](AvE1d)]:
if member(e[2],CC([n,partT],e[1])) then
   RETURN(G,partT,AvE minus {e}):
else
RETURN(G,partT union {e} ,AvE minus {e}):
fi:
FAIL:
end:

#MST(G): the minimal spanning tree of the weighted graph G=[n,E,DT]
#and the its cost (the sum of the costs of the edges of the cheapest tree)
MST:=proc(G) local n,E,A,i,DT,e:

n:=G[1]:
E:=G[2]:
DT:=G[3]:
if not CC([n,E],1)={seq(i,i=1..n)} then
 # print(`not connected`):
 RETURN(FAIL):
fi:

A:=G,{},E:
while nops(A[2])<n-1 and E<>{}  do
A:=MST1(A):
od:
[n,A[2]], add(DT[e], e in A[2]):
end:

#CC(G,a): inputs a graph G=[n,E] and a vertex a and outputs the set of vertices connected to
#a.For example CC([5,{{1,2},{2,3},{4,5}}],1) is {1,2,3}.
CC:=proc(G,a) local n,E,i,N, e,S,s,S1,S2:
n:=G[1]:
E:=G[2]:
#N[i]: the set of neighbors of vertex i
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:
S:={a}:
#S is a set and we want the set of all the neighbors of S, a typical set of neighbors is N[s]
#and we want the union of all of them
S1:=S union {seq( op(N[s])   , s in S)}:
while S<>S1 do
S2:=S1 union {seq( op(N[s])   , s in S1)}:
S:=S1:
S1:=S2:
od:
S1:
end:

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