#Nathan Fox
#Homework 27
#I give permission for this file to be posted online

##Read old files
read(`C27.txt`):

#Help procedure
Help:=proc():
print(` Wardrop(J,R,DelayJ,y,S,Sspec) `):
end:

with(LinearAlgebra):
with(Optimization):
with(ListTools):

##Problem 1
#Wardrop(J,R,DelayJ,y,S,Sspec): inputs
#J, a list of edges (in the form of pairs [i,j]),
#R, a list of Routes (in the form [i1,i2,i3, ..., ir] where i1 is
#the start, ir the end, and i2,i3, .. i(r-1) the intermediate
#cities),
#DelayJ (of the same length as the list J), of delay expressions
#phrased in terms of the variable y.
#a variable y (see above)
#a list S of allowed pairs [start,end] (possibly with only one
#member)
#A list Sspec, of the same size as S, indicating how many cars
#have to move from start to end for any pair in S (correponding)
#
#outputs a list of length nops(R) that indicate the flows along
#each route (corresponding to the list R), and another list with
#the corresponding delay times for each route.
#
#For example, for
#the wikipedia graph before the shortcut the function call should
#be (calling START: 1, A:2, B:3, END: 4)
#Wardrop([[1,2],[2,4],[1,3],[3,4]], [[1,2,4],[1,3,4]],
#[T/100,45,45,T/100],T,[[1,4]],[4000]);
#and with the "shortcut"
#Wardrop([[1,2],[2,4],[1,3],[3,4],[2,3]], [[1,2,4],[1,3,4],
#[1,2,3,4]],[T/100,45,45,T/100],T,[[1,4]],[4000]);
#and the outputs should be (respectively)
#[2000,2000],[65,65]
#and
#[0,0,4000],[85,85,80]
Wardrop:=proc(J,R,DelayJ,y,S,Sspec)
local H,A,i,j,k,l,xv,x,z,f,part1,part2,exp1,summ:
H:=Matrix(nops(S),nops(R),0):
A:=Matrix(nops(J),nops(R),0):
xv:=Matrix(nops(R),1,0):
for i from 1 to nops(R) do
    for j from 1 to nops(S) do
        if R[i][1] = S[j][1] and R[i][nops(R[i])] = S[j][2] then
            H[j,i]:=1:
        fi:
    od:
    for j from 1 to nops(J) do
        for k from 1 to nops(R[i]) - 1 do
            if R[i][k] = J[j][1] and R[i][k+1] = J[j][2] then
                A[j,i]:=1:
                break:
            fi:
        od:
    od:
    xv[i,1]:=x[i]:
od:
xv:=convert(xv,Vector):
z:=Multiply(A,xv):
f:=Multiply(H,xv):
part1:=Minimize(add(int(DelayJ[j],y=0..z[j]),j=1..nops(J)),
    {seq(x[j]>=0,j=1..nops(R)),seq(f[j]=Sspec[j],j=1..nops(S))})[2]:
for i from 1 to nops(part1) do
    part1[i]:=op(part1[i])[2]:
od:
part2:=[seq(0, i=1..nops(part1))]:
for i from 1 to nops(R) do
    summ:=0:
    for l from 1 to nops(R[i])-1 do
        j:=Search([R[i][l],R[i][l+1]],J):
        exp1:=subs(y=z[j],DelayJ[j]):
        exp1:=subs({seq(x[k]=part1[k],k=1..nops(part1))},exp1):
        summ:=summ+exp1:
    od:
    part2[i]:=summ:
od:
part1, part2:
end:

##Problem 2
#Input before addition is
#Wardrop([[1,2],[1,3],[2,4],[3,4]], [[1,2,4],[1,3,4]],
#[10*T,T+50,T+50,10*T],T,[[1,4]],[6]);
#Output is [3,3], [83,83], as expected
#Input after addition is
#Wardrop([[1,2],[1,3],[2,4],[3,4],[2,3]], [[1,2,4],[1,3,4],
#[1,2,3,4]], [10*T,T+50,T+50,10*T,T+10],T,[[1,4]],[6]);
#Output is [2,2,2], [92,92,92], as expected

##Problem 3
#This problem is open-ended enough that I don't feel I can do it
#justice in an assignment like this.  This could easily be the
#topic of a research paper.