######################################################################
## 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(`Version  : April 2021:  `):
print():
print(`This is Breaess.txt, A Maple package`):
print(`accompanying Shalsoh 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 `):

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 `):



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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,10*y]],y,6);`):
print(` FindEqu([[0,1+x,51+x/10,0],[0,0,10+x/10,51+x/10],[0,0,0,1+x]],x,60); `):
print(` FindEqu([[0,x/100,45,0],[0,0,1/100000,45],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,10*y]],y,6);`):
print(` FindEquC([[0,1+x,51+x/10,0],[0,0,10+x/10,51+x/10],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,10*y]],y,6);`):
print(` FindEquV([[0,1+x,51+x/10,0],[0,0,10+x/10,51+x/10],[0,0,0,1+x]],x,60); `):

elif nargs=1 and args[1]=IsBraess then
print(`IsBraess(G,x,A): Is G Braess? Try:`):
print(`IsBraess([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,10*y]],y,6);`):
print(` IsBraess([[0,1+x,51+x/10,0],[0,0,10+x/10,51+x/10],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,10*y]],y,6);`):
print(` IsBraessV([[0,1+x,51+x/10,0],[0,0,10+x/10,51+x/10],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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,0,0,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,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,0,y+50],[0,0,0,10*y]]);`):
print(`Routes([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,0,y+50],[0,0,0,10*y]],1);`):
print(`Routes1([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,0,y+50],[0,0,0,10*y]],1);
#Routes1([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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]<>0 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([[0,10*y,y+50,0],[0,0,0,y+50],[0,0,0,10*y]]);
#Routes([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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),0]:
   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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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]<>0 then
     G1:=[op(1..i-1,G),[op(1..j-1,G[i]),0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,10*y]],y,[6,0,0],1);
#
BestMove:=proc(G,x,C,i) local j,C1:

if C[i]=0 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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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]<>0 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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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]<>0 then
  print(`trying to delete the edge from `, i, ` to `, j):

     G1:=[op(1..i-1,G),[op(1..j-1,G[i]),0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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]<>0 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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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([[0,10*y,y+50,0],[0,0,y+10,y+50],[0,0,0,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 [[0,x/A,B,0],[0,0,epsilon,x/A],[0,0,0,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:=[[0,x/A,B,0],[0,0,1/10^9,x/A],[0,0,0,x/A]]:
G1:=[[0,x/A,B,0],[0,0,0,x/A],[0,0,0,x/A]]:


end: