######################################################################
## Braess.txt Save this file as  Braess.txt  to use it,              #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read Braess.txt <Enter>                           #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: April 2021: tested for Maple 2018 `):
print(`Current Version  : May 2022: Modified by Blair Seidler and Daniela Elizondo`):
print():
print(`This is Breaess.txt, A Maple package`):
print(`accompanying Shalosh B. Ekhad and Doron Zeilberger's  article: `):
print(` Experimenting with Braess' Amazing Paradox `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/XXX .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):

print(`----------------------------------`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print():
 print(`----------------------------------`):
print():
print(`For a list of the STORY functions type: ezraS();`):
print():
print(`For help with a specific function, type "ezra(procedure_name);" `):
print():
 print(`----------------------------------`):
print():

with(combinat):



ezraS:=proc()
if args=NULL then

print(`The STORY procedures are`):
print(`FindEquV, IsBraessV, IsEquV, TrajV  `):

else
ezra(args):

fi:


end:


ezraC:=proc()
if args=NULL then

print(`The the CONTINUOUS procedures are`):
print(`FindBestC, FindEquC `):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` BestComp, BestNei,BestMove, CarMatrix, Comps, MinMaxCS, Neis, RouteScore, RouteScore1, Routes, Routes1, TotRouteScore, Traj `):
print(` ErdosRenyi, RGrid, WattsStrogatz `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Braess.txt: A Maple package for  `):

print(`The MAIN procedures are`):
print(`  FindBraess, FindEqu, IsBraess,  IsEqu, NextState, RandNet `):
print(`  BuildRandGraph, BraessGameAuto, BraessGameAutoS`):

elif nargs=1 and args[1]=BuildRandGraph then
print(`BuildRandGraph(n,p): Builds a random network for use in the Braess Game below. Every potential`):
print(`edge has either value -1 (missing edge) or X (edge to be filled in during game)`):
print(`Try:`):
print(`BuildRandGraph(5, 0.7)`):

elif nargs=1 and args[1]=BraessGameAuto then
print(`BraessGameAuto(G,y,A): Has Alice and Bob fill in edges of a network with costs assigned as`):
print(`random linear functions. Outputs winner based on whether the final network is Braess.`):
print(`Try:`):
print(`BraessGameAuto(BuildRandGraph(5, 0.7),y,10)`):

elif nargs=1 and args[1]=BraessGameAutoS then
print(`BraessGameAutoS(G,y,A): Same as BraessGameAuto(G,y,A), but without printing the steps in the game.`):
print(`Returns true if Braess, false otherwise.`):
print(`Try:`):
print(`n:=BuildRandGraph(5, 0.7); for i to 1000 do if BraessGameAutoS(n, y, 6) then print(i): fi: od:`):

elif nargs=1 and args[1]=AddMinArc then
print(`AddMinArc(G,x,A): Given a network G with linear arc weights and A cars, adds an arc of weight a*x+b`):
print(`where a <= min({linear coefficient of arc weights}) and b <= min({constant term of arc weights})`):
print(`such that the new graph has the lowest average delay time possible, with the additional constraint that the new arc does not go from the source to the sink`):

elif nargs=1 and args[1]=RandOrigNet then
print(`RandOrigNet(A,x): Creates a network on 4 vertices with arcs 1->2, 1->3, 2->4, 3->4, and random linear arc weights`):
print(`where the coefficients in each weight are between 0 and A.`):

elif nargs=1 and args[1]=SimulateAddMinArc then 
print(`SimulateAddMinArc(K,n,r,A,x,C): Generates K random networks G := RandNet(n,r,A,x) with C cars.`):
print(`Counts the times that AddMinArc(G,x,C) has a smaller average delay time than G and returns this number/K.`):

elif nargs=1 and args[1]=SimulateAddOrig then
print(`SimulateAddOrig(K,A,C): Generates K random networks G := RandOrigNet(A,x) with C cars.`):
print(`Counts how often AddMinArc(G,x,C) improves traffic and returns this number/K.`):

elif nargs=1 and args[1]=ErdosRenyi then
print(`ErdosRenyi(n,c,x,p,q) generates a graph on n vertices under the Erdos-Renyi model.`):
print(`Edges has either value -1 or x with probability q, or c with probabilit 1-q`):

elif nargs=1 and args[1]=RGrid then
print(`RGrid(n,c,x,q) generates a grid on n^2 vertices.`):
print(`Vertex i is connected to vertex i+1 and i+n.`):
print(`Edges has either value -1 or x with probability q, or c with probabilit 1-q`):

elif nargs=1 and args[1]=WattsStrogatz then
print(`WattsStrogatz(n,c,x,p,q) generates a graph on n vertices under the Watts-Strogatz model, starting with a grid.`):
print(`Edges has either value -1 or x with probability q, or c with probabilit 1-q`):

elif nargs=1 and args[1]=BestMove then
print(`BestMove(G,x,C,i): inputs a traffic network G phrased in terms of the variable x, a current composition, C, describing`):
print(`how many cars use each of the available routes according to Routes(G), and i between 1 and nops(C), finds out whether someone who currently`):
print(`uses route R[i] can reduce its travel time, followed by the new composition of it exists. Otherwise it returns FAIL. Try:`):
print(`BestMove([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0],1);`):

elif nargs=1 and args[1]=BestNei then
print(`BestNei(G,x,C): The best neigbor of state C with network G phrased in terms of x. If none of the neigbors are better`):
print(`it returns C. Try`):
print(`BestNei([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);`):

elif nargs=1 and args[1]=CarMatrix then
print(`CarMatrix(G,C): inputs  a matrix G, and a vector C indicating how many cars choose each route`):
print(`outputs the list of lists , let's call it L, such that L[i][j] is the number of cars in the`):
print(`edge between i and j. Try:`):
print(`CarMatrix([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],[2,2,2]);`):

elif nargs=1 and args[1]=Comps then
print(`Comps(A,k): the set of vectors of non-negative integers of length k that add-up to A. Try:`):
print(`Comps(6,3);`):

elif nargs=1 and args[1]=FindBraess then
print(`FindBraess(n,r,A,x,K): Finds examples of Braess paradox with a network of n nodes and outdegree r with A cars trying up to K times`):
print(`Try: `):
print(` FindBraess(5,2,20,x,100); `):


elif nargs=1 and args[1]=FindBest then
print(`FindBest(G,x,A): Find the  configuration of routes with the traffic network G and A cars,`):
print(`followed by its value, Try`):
print(`FindBest([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);`):


elif nargs=1 and args[1]=FindBestC then
print(`FindBestC(G,x,A): Find the  configuration of routes with the traffic network G and A cars, where you are allowed fractional cars`):
print(`and usuing calcluls `):
print(``):
print(`followed by its value, Try`):
print(`FindBestC([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);`):

elif nargs=1 and args[1]=FindEqu then
print(`FindEqu(G,x,A): Finds an equilibrium for the network G phrased in terms of x, with A cars. `):
print(`The output is a list of length 6`):
print(`The first is the list of routes`):
print(`The second is the composition of A, telling you how many cars take the respetive route`):
print(`The third is the list if the delay times for each route (including those where there are zero cars`):
print(`The 4th is the average delay time`):
print(`The 5th is the minimum delay time of an actual car`):
print(`The 6th is the maximum delay time of an actual car`):
print(`Try:`):
print(`FindEqu([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);`):
print(` FindEqu([[-1,1+x,51+x/10,-1],[-1,-1,10+x/10,51+x/10],[-1,-1,-1,1+x]],x,60); `):
print(` FindEqu([[-1,x/100,45,-1],[-1,-1,1/100000,45],[-1,-1,-1,x/100]],x,4000); `):


elif nargs=1 and args[1]=FindEquC then
print(`FindEquC(G,x,A): Like FindEqu(G,x,A) (q.v.) but assuming a continous quantity of cars`):
print(`WARNING: Very slow and does not give a good answer. Use with caution`):
print(`FindEquC([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);`):
print(` FindEquC([[-1,1+x,51+x/10,-1],[-1,-1,10+x/10,51+x/10],[-1,-1,-1,1+x]],x,60); `):


elif nargs=1 and args[1]=FindEquV then
print(`FindEquV(G,x,A): verbose form of FindEqu(G,x,A) (q.v.) `):
print(`Try: `):
print(`FindEquV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);`):
print(` FindEquV([[-1,1+x,51+x/10,-1],[-1,-1,10+x/10,51+x/10],[-1,-1,-1,1+x]],x,60); `):

elif nargs=1 and args[1]=IsBraess then
print(`IsBraess(G,x,A): Is G Braess? Try:`):
print(`IsBraess([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);`):
print(` IsBraess([[-1,1+x,51+x/10,-1],[-1,-1,10+x/10,51+x/10],[-1,-1,-1,1+x]],x,60); `):


elif nargs=1 and args[1]=IsBraessV then
print(`IsBraessV(G,x,A): Verbose version if IsBraess(G,x,A) (q.v. Try:`):
print(`IsBraessV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);`):
print(` IsBraessV([[-1,1+x,51+x/10,-1],[-1,-1,10+x/10,51+x/10],[-1,-1,-1,1+x]],x,60); `):

elif nargs=1 and args[1]=IsEqu then
print(`IsEqu(G,x,C): Is the state C an equilibrium? Try`):
print(`IsEqu([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[2,2,2]);`):


elif nargs=1 and args[1]=IsEquV then
print(`IsEquV(G,x,C): Verbose form of IsEqu(G,x,C). Try`):
print(`IsEquV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[2,2,2]);`):

elif nargs=1 and args[1]=MinMaxCS then
print(`MinMaxCS(C,S): Given lists C and S of the same length outputs the pair [min,max]`):
print(`of S[i] where C[i]>0. Try:`):
print(`MinMaxCS([2,-1,-1,6],[1,2,3,4]);`):

elif nargs=1 and args[1]=Neis then
print(`Neis(C): all the neighbors of the compostion C. Try`):
print(`Neis([2,-1,-1,3]);`):

elif nargs=1 and args[1]=NextState then
print(`NextState(G,x,C): Given a traffic network G phrased in terms of x and a current composition C`):
print(`finds the next state or if it is an equilibrium, it returns FAIL. Try:`):
print(`NextState([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);`):

elif nargs=1 and args[1]=RandNet then
print(`RandNet(n,r,A,x): a random network without cycles with at most outdegree r,only edges beteen i->j i<j`):
print(`and traffic density function affine linear with coefficients between 0 and A. Try:`):
print(`RandNet(4,2,20,x);`):



elif nargs=1 and args[1]=RouteScore then
print(`RouteScore(G,x,C,R): Given a network G phrased in terms of the variable , and a composition C `):
print(`outputs the vector of delay times for all the routes, following the liset Routes(G). Try:`):
print(`RouteScore([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[2,2,2]);`):

elif nargs=1 and args[1]=RouteScore1 then
print(`RouteScore1(G,x,C,R): Given a network G phrased in terms of the variable , and a composition C and a route R`):
print(`outputs its delay time. Try:`):
print(`RouteScore1([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[2,2,2],[1,2,3,4]);`):


elif nargs=1 and args[1]=Routes then
print(`Routes(G): Given a network G (assumed without cycles)  outputs all`):
print(`routes from the source,1, to the sink, n. Try:`):
print(`Routes([[-1,10*y,y+50,-1],[-1,-1,-1,y+50],[-1,-1,-1,10*y]]);`):
print(`Routes([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]]);`):

elif nargs=1 and args[1]=Routes1 then
print(`Routes1(G,i): Given a network G (assumed without cycles) and a vertex i outputs the LIST of all`):
print(`routes from i to the sink n. Try:`):
print(`Routes1([[-1,10*y,y+50,-1],[-1,-1,-1,y+50],[-1,-1,-1,10*y]],1);`):
print(`Routes1([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],1);`):


elif nargs=1 and args[1]=TotRouteScore then
print(`TotRouteScore(G,x,C,R): Given a network G phrased in terms of the variable , and a composition C `):
print(`outputs the sum of the  delay times for all the routes. Try:`):
print(`TotRouteScore([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[2,2,2]);`):

elif nargs=1 and args[1]=Traj then
print(`Traj(G,x,C): The trajectory of states starting at C. Try:`):
print(`Traj([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);`):


elif nargs=1 and args[1]=TrajV then
print(`TrajV(G,x,C): Verbose form of Traj(G,x,C) (q.v.).  Try:`):
print(`TrajV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);`):


else
 print(`There is no such thing as`, args):

fi:


end:

#Routes1(G,i): Given a network G (assumed without cycles) and a vertex i outputs all
#routes from i to the sink n. Try:
#Routes1([[-1,10*y,y+50,-1],[-1,-1,-1,y+50],[-1,-1,-1,10*y]],1);
#Routes1([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],1);
Routes1:=proc(G,i) local n,gu,j,mu,mu1:
option remember:
n:=nops(G)+1:

if i=n then
 RETURN([[n]]):
fi:

gu:=[]:

for j from 1 to n do
 if G[i][j]<>-1 then
   mu:=Routes1(G,j):
    gu:=[op(gu),seq([i,op(mu1)],mu1 in mu)]:
 fi:
od:
gu:
end:


#Routes(G): Given a network G (assumed without cycles) and a vertex i outputs all
#routes from the source 1, to the sink n. Try:
#Routes([[-1,10*y,y+50,-1],[-1,-1,-1,y+50],[-1,-1,-1,10*y]]);
#Routes([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]]);
Routes:=proc(G) option remember: Routes1(G,1):end:



#RandNet(n,r,A,x): a random network without cycles with at most outdegree r,only edges beteen i->j i<j
#and traffic density function affine linear with coefficients between 0 and A. Try:
#RandNet(4,2,20,x);
RandNet:=proc(n,r,A,x) local ra1,ra2,G,i,j,G1,S,i1:
ra1:=rand(1..A):
ra2:=rand(0..A):

G:=[]:

for i from 1 to n-1 do
S:={seq(i1,i1=i+1..n)}:
S:=randcomb(S,min(r,nops(S))):
G1:=[]:

 for j from 1 to n do
   if not member(j,S) then
    G1:=[op(G1),-1]:
   else
    G1:=[op(G1),ra1()+ra2()*x]:
  fi:
od:
G:=[op(G),G1]:   

od:

G:

end:




#IsBraess(G,x,A): Is G Braess? Try:
#
#Try:
#IsBraess([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);
IsBraess:=proc(G,x,A) local G1,i,j,g,n:

if  FindEqu(G,x,A)=FAIL then
 RETURN(FAIL):
fi:

n:=nops(G)+1:

for i from 2 to n-2 do
 for j from 1 to nops(G[i]) do
  if G[i][j]<>-1 then
     G1:=[op(1..i-1,G),[op(1..j-1,G[i]),-1,op(j+1..nops(G[i]),G[i]) ], op(i+1..nops(G),G) ]:
if  FindEqu(G1,x,A)<>FAIL then

 if FindEqu(G1,x,A)[6]<FindEqu(G,x,A)[5] then
   RETURN(true, [i,j]):
 fi:


fi:

  fi:

 od:
od:
false,FAIL:


end:

#FindBraess(n,r,A,x,K): Finds examples of Braess paradox with a network of n nodes and outdegree r with parameters A and  K (for trying)
#Try:
#FindBraess(5,2,20,x,100);
FindBraess:=proc(n,r,A,x,K) local G,i:

G:=RandNet(n,r,A,x,K):

if IsBraess(G,x,A)[1] then
 RETURN(G, IsBraess(G,x,A)[2]):
fi:

for i from 1 to K while not IsBraess(G,x,A)[1] do od:
if i=K+1 then
 RETURN(FAIL):
else
G:=RandNet(n,r,A,x,K):
RETURN([G, IsBraess(G,x,A)[2]]):
fi:

end:


#Comps(A,k): the set of vectors of non-negative integers of length k that add-up to A. Try:
#Comps(6,3);
Comps:=proc(A,k) local gu,mu,i,mu1:

if k=0 then
 if A=0 then
    RETURN({[]}):
 else
  RETURN({}):
 fi:

fi:

gu:={}:

for i from 0 to A do
 mu:=Comps(A-i,k-1):
 gu:=gu union {seq([op(mu1),i],mu1 in mu)}:
od:

gu:

end:


#CarMatrix(G,C): inputs  a matrix G, and a vector C indicating how many cars choose each route
#outputs the list of lists , let's call it L, such that L[i][j] is the number of cars in the
#edge between i and j. Try:
#CarMatrix([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],[2,2,2]);
CarMatrix:=proc(G,C) local n,i,j,T,R,i1,i2:
option remember:
 n:=nops(G)+1:
 for i  from 1 to n-1 do
  for j from 1 to n do
   T[i,j]:=0:
 od:
od:

R:=Routes(G):

if nops(C)<>nops(R) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

for i1 from 1 to nops(R) do
 for i2 from 1 to nops(R[i1])-1 do
 
    T[R[i1][i2],R[i1][i2+1]]:=T[R[i1][i2],R[i1][i2+1]]+ C[i1]:
 od:
od:
  
[seq([seq(T[i,j],j=1..n)],i=1..n-1)]:

end:
 
#RouteScore1(G,x,C,R): Given a network G phrased in terms of the variable , and a composition C and a route R
#outputs its delay time
RouteScore1:=proc(G,x,C,R)  local CM,i1:
CM:=CarMatrix(G,C):

add(subs(x=CM[R[i1]][R[i1+1]],G[R[i1]][R[i1+1]]),i1=1..nops(R)-1):

end:



#RouteScore(G,x,C): Given a network G phrased in terms of the variable , and a composition C outputs
#the vector of Route-scores corresponding to all routes using the list Routes(G)
RouteScore:=proc(G,x,C)  local R,i1:
R:=Routes(G):
[seq(RouteScore1(G,x,C,R[i1]) ,i1=1..nops(R))]:
end:


#TotRouteScore(G,x,C): Given a network G phrased in terms of the variable , and a composition C outputs
#the total of the Route-scores
TotRouteScore:=proc(G,x,C)  local R,i1:
R:=Routes(G):
add(C[i1]*RouteScore1(G,x,C,R[i1]) ,i1=1..nops(R)):
end:

#BestMove(G,x,C,i): inputs a traffic network G phrased in terms of the variable x, a current composition, C, describing
#how many cars use each of the available routes according to Routes(G), and i between 1 and nops(C), finds out whether someone who currently
#uses route R[i] can reduce its travel time, followed by the new composition of it exists. Otherwise it returns FAIL. Try:
#BestMove([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0],1);
#
BestMove:=proc(G,x,C,i) local j,C1:

if C[i]=-1 then
 RETURN(FAIL):
fi:

for j from 1 to i-1 do
  C1:=[op(1..j-1,C),C[j]+1,op(j+1..i-1,C),C[i]-1,op(i+1..nops(C),C)]:
   if RouteScore(G,x,C1)[j]<RouteScore(G,x,C)[i] then
     RETURN(C1):
   fi:
od:


for j from i+1 to nops(C) do
  C1:=[op(1..i-1,C),C[i]-1,op(i+1..j-1,C),C[j]+1,op(j+1..nops(C),C)]:
   if RouteScore(G,x,C1)[j]<RouteScore(G,x,C)[i] then
     RETURN(C1):
   fi:
od:

FAIL:

end:

#NextState(G,x,C): Given a traffic network G phrased in terms of x and a current composition C
#finds the next state or if it is an equilibrium, it returns FAIL. Try:
#NextState([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);
NextState:=proc(G,x,C) local i,C1:

for i from 1 to nops(C) do
 C1:=BestMove(G,x,C,i):
 if C1<>FAIL then
   RETURN(C1):
 fi:
od:

FAIL:

end:


#IsEqu(G,x,C): Is the state C an equilibrium? Try
#IsEqu([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[2,2,2]);
IsEqu:=proc(G,x,C) local i,j,C1:

for i from 1 to nops(C) do
if C[i]>0 then
 for j from 1 to i-1 do
  C1:=[op(1..j-1,C),C[j]+1,op(j+1..i-1,C),C[i]-1,op(i+1..nops(C),C)]:
    if RouteScore(G,x,C1)[j]<RouteScore(G,x,C)[i] then
     RETURN(false):
    fi:
 od:


 for j from i+1 to nops(C) do
  C1:=[op(1..i-1,C),C[i]-1,op(i+1..j-1,C),C[j]+1,op(j+1..nops(C),C)]:
    if RouteScore(G,x,C1)[j]<RouteScore(G,x,C)[i] then
     RETURN(false):
    fi:
 od:
fi:
od:

true:

end:



#IsEquV(G,x,C): Verbose form of IsEqu(G,x,C) (q.v.). Try:
#IsEquV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[2,2,2]);
IsEquV:=proc(G,x,C) local n,R,i,j,C1,S,CM:

R:=Routes(G):

n:=nops(G)+1:
print(``):


print(`Our traffic network has `, nops(G)+1, `cities, called`, seq(i,i=1..n), `and traffic flows from City 1 to City`, n):
print(``):
print(`There are`, convert(C,`+`), `cars, all starting at city 1 and all of them are travelling to their destination in city `, n):
print(``):
for i from 1 to n-1 do
print(`Out of City`, i, `we have  one-way roads to the following cities`):
 for j from 2 to n do
  if G[i][j]<>-1 then
   print(`From City `, i, `To City `, j, `. If there are `, x, ` cars on that road, then  the travel time is `, G[i][j] , `minutes `):
  fi:
 od:
od:


S:=RouteScore(G,x,C):
print(``):
print(`Here are the available routes, followed by the number of cars currently using them, and the time of travel for those using them `):
print(``):

for i from 1 to nops(R) do
print(`Route number`, i, ` is the the route`, R[i], `and currently the number of cars using it is`, C[i]):
od:

CM:=CarMatrix(G,C):

print(``):
print(`This implies that the car matrix is`):
print(``):
print(matrix(CM)):
print(``):

print(`This implies that the current travel times for each of the routes are  as follows`):
print(``):
for i  from 1 to nops(R) do
print(`For route`, i, `namely `, R[i], `the `, C[i], `cars using it each has travel time `, S[i], `minutes `):
od:
print(``):

print(`Let's see if any one can change  from its current route to another one and reduce his travel time `):

for i from 1 to nops(C) do

if C[i]>0 then
print(``):
print(`---------------------------------------`):
print(`---------------------------------------`):
print(``):
 print(`One of the`, C[i], `drivers that are currently using route`, R[i], `is considering switching to another route `):
print(``):

 for j from 1 to i-1 do

  C1:=[op(1..j-1,C),C[j]+1,op(j+1..i-1,C),C[i]-1,op(i+1..nops(C),C)]:

print(``):
print(`------------------------------------------------------------------------`):
print(``):

 print(`If he changes to route`, R[j], `then that route will have one more car, hence`, C1[j], `cars and the current route`, R[i]):
print(`would have one less car, namely`, C1[i]):
print(``):
print(`the new car matrix is `):
 print(``):
 print(matrix(CarMatrix(G,C1))):
 print(``):

    if RouteScore(G,x,C1)[j]<RouteScore(G,x,C)[i] then

  print(`Since the new  travel time`,  RouteScore(G,x,C1)[j], `is LESS than the previous travel time`,  RouteScore(G,x,C)[i],  `he would change to that route`):
   print(`and the configuration`, C, ` is NOT an equilibrium`):

     RETURN(false):

    else
  print(`Since the new  travel time`,  RouteScore(G,x,C1)[j], `is Lareger or equal than the previous travel time`,    RouteScore(G,x,C)[i]  ,`he has no reason to change`):
    fi:

 od:

 for j from i+1 to nops(C) do

  C1:=[op(1..i-1,C),C[i]-1,op(i+1..j-1,C),C[j]+1,op(j+1..nops(C),C)]:

 print(`If he changes to route`, R[j], `then that route will have one more car, hence`, C1[j], `cars and the current route`, R[i]):
print(``):
print(`would have one less car, namely`, C1[i]):
print(``):
print(`the new car matrix is `):
 print(``):
 print(matrix(CarMatrix(G,C1))):
 print(``):

    if RouteScore(G,x,C1)[j]<RouteScore(G,x,C)[i] then
  print(`Since the new  travel time`,  RouteScore(G,x,C1)[j], `is LESS than the previous travel time`, RouteScore(G,x,C)[i], `he would change to that route`):
   print(`and the configuration`, C, ` is NOT an equilibrium`):

     RETURN(false):
else
  print(`Since the new  travel time`,  RouteScore(G,x,C1)[j], `is Lareger or equal than the previous travel time`, RouteScore(G,x,C)[i]  , `he has no reason to change`):
    fi:

 od:
fi:

od:
true:

end:


#Traj(G,x,C): The trajectory of states starting at C. Try:
#Traj([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);
Traj:=proc(G,x,C) local gu,C1,C2:
C1:=C:
gu:=[C1]:
C2:=NextState(G,x,C1):

while C2<>FAIL do
C1:=C2:

gu:=[op(gu),C1]:

C2:=NextState(G,x,C1):
od:

gu:

end:




#FindEquSlow(G,x,A): Finds an equilibrium for the network G phrased in terms of x, with A cars. Try:
#FindEquSlow([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);
FindEquSlow:=proc(G,x,A) local S,i,R,K,av,C,i1:

R:=Routes(G):

S:=Comps(A,nops(R)):

for i from 1 to nops(S) do
 if IsEqu(G,x,S[i]) then
   C:=S[i]:
 K:=RouteScore(G,x,C):
av:=evalf(add(C[i1]*K[i1],i1=1..nops(C))/A):
   RETURN([R,S[i],K,av,op(MinMaxCS(S[i],K))]):
 fi:
od:

FAIL:
end:


#FindEqu(G,x,A): Finds an equilibrium for the network G phrased in terms of x, with A cars. Try:
#FindEqu([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);
FindEqu:=proc(G,x,A) local R,C,S,av,i:

R:=Routes(G):

C:=Traj(G,x,[0$(nops(R)-1),A])[-1]: 

if not IsEqu(G,x,C) then
 RETURN(FAIL):
fi:

S:=RouteScore(G,x,C):
av:=evalf(add(S[i]*C[i],i=1..nops(C))/A):
[R,C,S,av,op(MinMaxCS(C,S))]:

end:

#MinMaxCS(C,S): Given lists C and S of the same length outputs the pair [min,max]
#of S[i] where C[i]>0. Try:
#MinMaxCS([2,0,0,6],[1,2,3,4]);
MinMaxCS:=proc(C,S) local S1,i:
if not nops(C)=nops(S) then
 print(`bad input`):
 RETURN(FAIL):
fi:

S1:=[]:

for  i from 1 to nops(C) do
if C[i]<>0 then
 S1:=[op(S1),S[i]]:
fi:
od:
[min(S1),max(S1)]:
end:


#FindBestSlow(G,x,A): Find the  configuration of routes with the traffic network G and A cars,
#followed by its value
FindBestSlow:=proc(G,x,A) local R,S,champ,rec,hope1,C,i:
R:=Routes(G):
S:=Comps(A,nops(R)):

champ:={[S[1], RouteScore(G,x,S[1])]}:

rec:=TotRouteScore(G,x,S[1]):


for i from 2 to nops(S) do
 C:=S[i]:
 hope1:=TotRouteScore(G,x,C):
 if hope1<rec then
   champ:={[C,RouteScore(G,x,C)]}:
   rec:=hope1:
 elif hope1=rec then
  champ:=champ union {[C,RouteScore(G,x,C)]}:
 fi:
od:
[R,champ, evalf(rec/A)]:
end:


#IsBraessV(G,x,A): Verbose versiob Is G Braess? Try:
#
#Try:
#IsBraessV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);
IsBraessV:=proc(G,x,A) local R,G1,i,j,g,n:


if  FindEqu(G,x,A)=FAIL then
 RETURN(FAIL):
fi:

R:=Routes(G):
print(`The available routes are`):
print(``):
lprint(R):
print(``):
n:=nops(G)+1:

print(`The network is`):
lprint(G):

print(`The current delay times are`):

lprint(evalf(FindEqu(G,x,A)[3])):

for i from 2 to n-2 do
 for j from 1 to nops(G[i]) do
  if G[i][j]<>-1 then
  print(`trying to delete the edge from `, i, ` to `, j):

     G1:=[op(1..i-1,G),[op(1..j-1,G[i]),-1,op(j+1..nops(G[i]),G[i]) ], op(i+1..nops(G),G) ]:

print(`Now the network is`):
print(``):
lprint(G1):

if  FindEqu(G1,x,A)<>FAIL then

print(`Now the delay times are`):

lprint(evalf(FindEqu(G1,x,A)[3])):

print(``):

 if max(FindEqu(G1,x,A)[3])<min(FindEqu(G,x,A)[3]) then
   RETURN(true, [i,j]):
 fi:


fi:

  fi:

 od:
od:
false,FAIL:


end:





#BestMoveV(G,x,C,i): inputs a traffic network G phrased in terms of the variable x, a current composition, C, describing
#how many cars use each of the available routes according to Routes(G), and i between 1 and nops(C), finds out whether someone who currently
#uses route R[i] can reduce its travel time, followed by the new composition of it exists. Otherwise it returns FAIL. Try:
#BestMove([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0],1);
#
BestMoveV:=proc(G,x,C,i) local R,j,C1:

if C[i]=0 then
 RETURN(FAIL):
fi:

R:=Routes(G):

print(``):	
print(`----------------------------------------------------------`):
print(``):
print(`The current state is`, C):
print(`This means that there are`, C[i], `drivers who use the route number `,i, `namely the route`, R[i]):
print(``):
print(`Let's see whether it is a good idea for such a driver to switch to another route`):
print(``):

for j from 1 to i-1 do
print(`Let's see what happens if he switces to route number`, j, `in other words the route`, R[j]):
print(``):
  C1:=[op(1..j-1,C),C[j]+1,op(j+1..i-1,C),C[i]-1,op(i+1..nops(C),C)]:
print(``):
print(`Now the car matrix is`):
print(``):
print(matrix(CarMatrix(G,C1))):
print(``):

   if RouteScore(G,x,C1)[j]<RouteScore(G,x,C)[i] then
      print(`Indeed now the travel time of the driver who just switched to route number`, j, `is `,  RouteScore(G,x,C1)[j], `that is smaller than `):
      print(`what it was before, namely`, RouteScore(G,x,C)[i], ` so it a good idea to switch and now we are in state`, C1):
     RETURN(C1):
  else
      print(`Since the travel time of the driver who just switched to route number`, j, `is `,  RouteScore(G,x,C1)[j], `is NOT smaller `):
      print(`than what it was before, namely`, RouteScore(G,x,C)[i], ` it is NOT a good idea to switch`):

   fi:

od:



for j from i+1 to nops(C) do
  C1:=[op(1..i-1,C),C[i]-1,op(i+1..j-1,C),C[j]+1,op(j+1..nops(C),C)]:


   if RouteScore(G,x,C1)[j]<RouteScore(G,x,C)[i] then
      print(`Indeed now the travel time of the driver who just switched to route number`, j, `is `,  RouteScore(G,x,C1)[j], `that is smaller than `):
      print(`what it was before, namely`, RouteScore(G,x,C)[i], ` so it a good idea to switch and now we are in state`, C1):
     RETURN(C1):
  else
      print(`Since the travel time of the driver who just switched to route number`, j, `is `,  RouteScore(G,x,C1)[j], `is NOT smaller `):
      print(`than what it was before, namely`, RouteScore(G,x,C)[i], ` it is NOT a good idea to switch`):

   fi:
od:

print(`None of the `, C[i], `drivers who currently use route number`, i, `namely the route `, R[i], `have any incentives to swicth `):
FAIL:

end:

#NextStateV(G,x,C): Verbose form of  NextState(G,x,C) (q.v.). Try:
#NextStateV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);
NextStateV:=proc(G,x,C) local i,C1:

for i from 1 to nops(C) do
 C1:=BestMoveV(G,x,C,i):
 if C1<>FAIL then
   RETURN(C1):
 fi:
od:

FAIL:

end:


#TrajV(G,x,C): Verbose version of Traj(G,x,C) (q.v.). Try:
#TrajV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);
TrajV:=proc(G,x,C) local gu,C1,C2, CM, R, S, i, j, n:


R:=Routes(G):

n:=nops(G)+1:
print(``):


print(`Our traffic network has `, nops(G)+1, `cities, called`, seq(i,i=1..n), `and traffic flows from City 1 to City`, n):
print(``):
print(`There are`, convert(C,`+`), `cars, all starting at city 1 and all of them are travelling to their destination in city `, n):
print(``):

for i from 1 to n-1 do
print(`Out of City`, i, `we have  one-way roads to the following cities`):
 for j from 2 to n do
  if G[i][j]<>-1 then
   print(`From City `, i, `To City `, j, `. If there are `, x, ` cars on that road, then  the travel time is `, G[i][j] , `minutes `):
  fi:
 od:
od:


print(`We start with the state`, C, `and see whether at least one driver can do better each time, until we reach an equilibrium`):
print(``):
S:=RouteScore(G,x,C):
print(``):
print(`Here are the available routes, followed by the number of cars in the intial state using them, and the time of travel for those using them `):
print(``):

for i from 1 to nops(R) do
print(`Route number`, i, ` is the the route`, R[i], `and currently the number of cars using it is`, C[i]):
od:

CM:=CarMatrix(G,C):

print(``):
print(`This implies that the car matrix is`):
print(``):
print(matrix(CM)):
print(``):


C1:=C:
gu:=[C1]:
C2:=NextStateV(G,x,C1):


while C2<>FAIL do
C1:=C2:

gu:=[op(gu),C1]:

C2:=NextStateV(G,x,C1):
od:

gu:

end:






#FindEquV(G,x,A): Verbose form of FindEqu(G,x,A). Try:
#FindEquV([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,6);
FindEquV:=proc(G,x,A) local R,C,S,av,i:

R:=Routes(G):

print(`We start with state`,  [0$(nops(R)-1),A] ):


C:=TrajV(G,x,[0$(nops(R)-1),A])[-1]: 

if not IsEqu(G,x,C) then
 RETURN(FAIL):
fi:
print(` State `, C, `is an equilibrium since no driver has an incentive to switch `):
print(``):
S:=RouteScore(G,x,C):
av:=evalf(add(S[i]*C[i],i=1..nops(C))/A):

print(`To sum up `):
print(``):
print(`In the equilibrium `):
print(``):

for i from 1 to nops(R) do
print(C[i] , `cars are driving the route `, R[i], `and these drivers take`, S[i], `minutes to get make the trip`):
print(``):
print(`The average driving time is`, av ):
print(``):
od:

[R,C,S,av,op(MinMaxCS(C,S))]:

end:



#Neis(C): all the neighbors of the compostion C. Try
#Neis([2,0,0,3]);
Neis:=proc(C) local i,j,gu:
gu:={}:

for i from 1 to nops(C) do

if C[i]>0 then


for j from 1 to i-1 do
  gu:=gu union {[op(1..j-1,C),C[j]+1,op(j+1..i-1,C),C[i]-1,op(i+1..nops(C),C)]}:

 od:

for j from i+1 to nops(C) do

  gu:=gu union {[op(1..i-1,C),C[i]-1,op(i+1..j-1,C),C[j]+1,op(j+1..nops(C),C)]}:

 od:
fi:
od:
gu:
end:


#BestNei(G,x,C): The best neigbor of state C with network G phrased in terms of x. If none of the neigbors are better
#it returns C. Try
#BestNei([[-1,10*y,y+50,-1],[-1,-1,y+10,y+50],[-1,-1,-1,10*y]],y,[6,0,0]);`):
BestNei:=proc(G,x,C) local aluf, rec,gu,hope1,gu1:

aluf:=C:

rec:=TotRouteScore(G,x,C):

gu:=Neis(C):

for gu1 in gu do
 hope1:=TotRouteScore(G,x,gu1):
 if hope1<rec then
  aluf:=gu1:
  rec:=hope1:
 fi:
od:

aluf:

end:



#FindBest(G,x,A): Find the  configuration of routes with the traffic network G and A cars,
#followed by its value
FindBest:=proc(G,x,A) local R,C,C1,C2,rec:

R:=Routes(G):
C:=[A,0$(nops(Routes(G))-1)]:
C1:=BestNei(G,x,C):

while C<>C1 do
C:=C1:
C1:=BestNei(G,x,C):
od:


rec:=TotRouteScore(G,x,C1):

[R,C1, RouteScore(G,x,C1),evalf(rec/A)]:
end:


###start continuous version

#FindBestC(G,x,A): Find the  configuration of routes with the traffic network G and A cars,
#followed by its value. USING continuours cars, and calculus
FindBestC:=proc(G,x,A) local R,c,i,f,eq,C,gu:

R:=Routes(G):


C:=[seq(c[i],i=1..nops(R))]:

f:=TotRouteScore(G,x,C)/A:

eq:={add(c[i],i=1..nops(C))=A}:

gu:=Optimization[Minimize](f,eq,assume=nonnegative);

C:=subs(gu[2],C):

[R,C, RouteScore(G,x,C),gu[1]]:

end:



#FindEquC(G,x,A): Find the  equilibria of routes with the traffic network G and A cars,
#followed by its value
FindEquC:=proc(G,x,A) local R,c,i,j,eq,C,var,C1:

R:=Routes(G):


C:=[seq(c[i],i=1..nops(R))]:

eq:={add(c[i],i=1..nops(C))=A, seq(c[i]>=0,i=1..nops(R))}:

var:={seq(c[i],i=1..nops(R))}:

for i from 1 to nops(C) do
 for j from 1 to i-1 do
  C1:=[op(1..j-1,C),C[j]+1,op(j+1..i-1,C),C[i]-1,op(i+1..nops(C),C)]:
    eq:=eq union { RouteScore(G,x,C1)[j]-RouteScore(G,x,C)[i]>=0}
 od:


 for j from i+1 to nops(C) do
  C1:=[op(1..i-1,C),C[i]-1,op(i+1..j-1,C),C[j]+1,op(j+1..nops(C),C)]:
    eq:=eq union { RouteScore(G,x,C1)[j]-RouteScore(G,x,C)[i]>=0}
 od:

od:

solve(eq,var):

end:


#FindBraessRange(A,B,N) In the network [[-1,x/A,B,-1],[-1,-1,epsilon,x/A],[-1,-1,-1,x/A]] what is the
#range of number of cars from 1 to N that give the Braess paradox, i.e. removing epsilon would make it better
FindBraessRange:=proc(A,B,N) local G,G1:
G:=[[-1,x/A,B,-1],[-1,-1,1/10^9,x/A],[-1,-1,-1,x/A]]:
G1:=[[-1,x/A,B,-1],[-1,-1,-1,x/A],[-1,-1,-1,x/A]]:


end:

# AddMinArc(G,x,A): Given a network G with linear arc weights and A cars, adds an arc of weight a*x+b
# where a <= min({linear coefficient of arc weights}) and b <= min({constant term of arc weights})
# such that the new graph has the lowest average delay time possible, with the additional constraint that the new arc does not go from the source to the sink
AddMinArc := proc(G,x,A) local best, n, linear, constant, i, j, a, b, G1:
best := G:
n := nops(G) + 1:
linear := min({seq(seq(coeff(G[i,j], x), j=1..n), i=1..n-1)} minus {-1}):
constant := min({seq(seq(subs(x=0, G[i,j]), j=1..n), i=1..n-1)} minus {-1}):

for i from 1 to n-1 do:
    for j from 1 to n do:
        if i < j and G[i][j] = -1 then
            if i <> 1 or j <> n then
                for a from 0 to linear do:
                    for b from 0 to constant do:
                        G1 := [op(1..i-1,G),[op(1..j-1, G[i]), a*x+b, op(j+1..n, G[i])] ,op(i+1..n-1, G)]:
                            if FindEqu(G1,x,A)[4] < FindEqu(best, x, A)[4] then
                                best := G1:
                            fi:
                    od:
                od:
            fi:
        fi:     
    od:
od:
best:
end:

# Creates a network on 4 vertices with arcs 1->2, 1->3, 2->4, 3->4, and random linear arc weights
# where the coefficients in each weight are between 0 and A. Returns FAIL if each arc weight is 0.
RandOrigNet := proc(A,x) local ra:
ra := rand(0..A):
[[-1, ra()*x+ra(), ra()*x+ra(), -1], [-1,-1,-1, ra()*x+ra()], [-1,-1,-1, ra()*x+ra()]]:
end:

# SimulateAddMinArc(K,n,r,A,x,C): Generates a random network G := RandNet(n,r,A,x) with C cars K times and
# counts the times that AddMinArc(G,x,C) has a smaller average delay time than G. Returns this number/K.
SimulateAddMinArc := proc(K,n,r,A,x,C) local G, i, better:
better := 0:
for i from 1 to K do:
    G := RandNet(n,r,A,x):
    if AddMinArc(G,x,C) <> G then
       better += 1: fi:
od:
evalf(better/K):
end:

# SimulateAddOrig(K,A,C): Generates K random networks G := RandOrigNet(A,x) with C cars and
# counts how often AddMinArc(G,x,C) improves traffic. Returns this number/K.
SimulateAddOrig := proc(K,A,x,C) local G, i, better:
better := 0:
for i from 1 to K do:
    G := RandOrigNet(A,x):
    if AddMinArc(G,x,C) <> G then 
        better += 1:
    fi: 
od:
evalf(better/K):
end:

# BuildRandGraph(n,p): Builds a random network for use in the Braess Game below. Every potential
# edge has either value -1 (missing edge) or X (edge to be filled in during game)
BuildRandGraph:=proc(n,p) local i,j,r,G,g:
r:=rand(0..1.0):
G:=[]:
for i from 1 to n-1 do
  g:=[-1$n]:
  for j from i+1 to n do
    if r()<p then
      g[j]:=X:
    fi:
  od:
  if {op(g)}={-1} then g:=[-1$(i),X,-1$(n-i-1)]: fi:
  G:=[op(G),g]:
od:
G:
end:

# BraessGameAuto(G,y,A): Has Alice and Bob fill in edges of a network with costs assigned as
# random linear functions. Outputs winner based on whether the final network is Braess
BraessGameAuto:=proc(G,y,A) local P,p,i,j,g,r,f,n,IB:
n:=nops(G[1]):
P:=["Alice","Bob"]:
p:=1:
g:=G:
r:=rand(0..10):
for i from 1 to n-1 do
  for j from 1 to n do
    if g[i][j]<>-1 then
      f:=r()*y+r():
      printf("%s assigns %s to link [%d,%d]\n",P[p],convert(f,string),i,j):
      g[i][j]:=f:
      p:=3-p:
      printf("Current network: %s\n",convert(g,string)):
    fi:
  od:
od:
IB:=IsBraess(g,y,A):
if IB[1] then
  printf("Network is Braess - removal of %s improves flow - %s wins!\n",convert(IB[2],string),P[2]):
else
  printf("Network is not Braess - %s wins!\n",P[1]):
fi:
IB:
end:

# BraessGameAutoS(G,y,A): Same as BraessGameAuto(G,y,A), but without printing the steps in the game.
# Returns true if Braess, false otherwise.
BraessGameAutoS:=proc(G,y,A) local P,p,i,j,g,r,f,n,IB:
n:=nops(G[1]):
P:=["Alice","Bob"]:
p:=1:
g:=G:
r:=rand(0..10):
for i from 1 to n-1 do
  for j from 1 to n do
    if g[i][j]<>-1 then
      f:=r()*y+r():
#      printf("%s assigns %s to link [%d,%d]\n",P[p],convert(f,string),i,j):
      g[i][j]:=f:
      p:=3-p:
#      printf("Current network: %s\n",convert(g,string)):
    fi:
  od:
od:
IB:=IsBraess(g,y,A):
(*if IB[1] then
  printf("Network is Braess - removal of %s improves flow - %s wins!\n",convert(IB[2],string),P[2]):
else
  printf("Network is not Braess - %s wins!\n",P[1]):
fi:*)
IB[1]:
end:

# ErdosRenyi(n,c,x,p,q) generates a graph on n vertices under the Erdos-Renyi model.
# Edges are present with probability p.
# Costs on edges are either x with probability q or c with probability 1-q.
ErdosRenyi:=
    proc(n,c,x,p,q)
        local i,j:
        local ra:= rand(1..denom(p)):
        local rc:= rand(1..denom(q)):
        [seq([seq(`if`(i<j and ra()<=numer(p),`if`(rc()<=numer(q),x,c),-1),j=1..n)],i=1..n-1)]:
    end:

(*

# Verifying Erdos-Renyi being Braess.

p:=3/4:
q:=1/2:
for n from 3 to 10 do
    print(IsBraess(ErdosRenyi(n,1,x,p,q),x,ceil(n*p*q)));
od:

*)

# RGrid(n,c,x,q) generates a grid on n^2 vertices.
# Edges are directed down and to the right.
# Costs on edges are either x with probability q or c with probability 1-q.
RGrid:=
    proc(n,c,x,q)
        local i,j:
        local rc:= rand(1..denom(q)):
        local G:= [seq([seq(-1,i=1..n^2)],j=1..n^2-1)]:
        for i from 1 to n^2-1 do
            if i mod n<>0 then
                G[i][i+1]:=`if`(rc()<=numer(q),x,c):
            fi:
            if i+n<=n^2 then
                G[i][i+n]:=`if`(rc()<=numer(q),x,c):
            fi:
        od:
        G:
    end:

(*

# Checking whether grids are Braess.

p:=3/4:
q:=1/2:
for n from 3 to 10 do
    print(IsBraess(RGrid(n,1,x,p,q),x,ceil(n*p*q)));
od:

*)

# WattsStrogatz(n,c,x,p,q) generates a graph on n vertices under the Watts-Strogatz model, starting from a grid.
# Costs on edges are either x with probability q or c with probability 1-q.
WattsStrogatz:=
    proc(n,c,x,p,q)
        local i,j,k,r:
        local ra:= rand(1..denom(p)):
        local G:= RGrid(n,c,x,q):
        for i from 1 to n^2-2 do
            for j from 1 to n^2-1 do
                if G[i][j]<>-1 and ra()<=numer(p) then
                    r:=op(randcomb({seq(k,k=j+1..n^2)},1)):
                    while G[i][r]<>-1 do
                        r:=op(randcomb({seq(k,k=j+1..n^2)},1)):
                    od:
                    G[i][r]:= G[i][j]:
                    G[i][j]:=-1:
                fi:
            od:
        od:
        G:
    end:

(*

# Checking whether Watts-Strogatz networks are Braess.

p:=3/4:
q:=1/2:
for n from 3 to 10 do
    print(IsBraess(WattsStrogatz(n,1,x,p,q),x,ceil(n*p*q)));
od:

*)