-
Notifications
You must be signed in to change notification settings - Fork 0
/
Social_Optimum.py
executable file
·158 lines (133 loc) · 5.79 KB
/
Social_Optimum.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
__author__ = "Jerome Thai"
__email__ = "[email protected]"
'''
This module is frank-wolfe algorithm using an all-or-nothing assignment
based on igraph package
'''
import numpy as np
from process_data import construct_igraph, construct_od
from AoN_igraph import all_or_nothing
#Profiling the code
import timeit
#This finds the potential for social optimum case
def potential_Social_Optimum(graph ,f):
#Finds the potential in the case of social optimum
links = int(np.max(graph[:,0])+1)
#The potential in social optimum is a0*f+a1*f^2+a2*f^3+a3*f^4+a4*a^5
x = np.power(f.reshape((links,1)), np.array([1,2,3,4,5]))
return np.sum(np.einsum('ij,ij->i', x, graph[:,3:]))
def line_search(f, res=20):
# on a grid of 2^res points bw 0 and 1, find global minimum
# of continuous convex function
d = 1./(2**res-1)
l, r = 0, 2**res-1
while r-l > 1:
if f(l*d) <= f(l*d+d): return l*d
if f(r*d-d) >= f(r*d): return r*d
# otherwise f(l) > f(l+d) and f(r-d) < f(r)
m1, m2 = (l+r)/2, 1+(l+r)/2
if f(m1*d) < f(m2*d): r = m1
if f(m1*d) > f(m2*d): l = m2
if f(m1*d) == f(m2*d): return m1*d
return l*d
def total_free_flow_cost(g, od):
return np.array(g.es["weight"]).dot(all_or_nothing(g, od))
#Calculates the total travel cost/time
def total_cost(graph, f, grad):
#Since the cost function equals to t(f) = a0+a1*f+a2*f^2+a3*f^3+a4*f^4 (where f is the flow)
#the travel cost, f*t(f) = a0*f+ a1*f^2 + a2*f^3+ a3*f^4 + a4*f^5
x = np.power(f.reshape((f.shape[0],1)), np.array([1,2,3,4,5])) # x is a matrix containing f,f^2, f^3, f^4, f^5
tCost = np.sum(np.einsum('ij,ij->i', x, graph[:,3:])) # Multiply matrix x with coefficients a0, a1, a2, a3 and a4
return tCost
#This function allow each user to account for their impact on the global travel cost
#thus allowing to find the social optimum
def search_direction_social_optimum(f, graph, g, od):
# computes the Frank-Wolfe step
# g is just a canvas containing the link information and to be updated with
# the most recent edge costs
x = np.power(f.reshape((f.shape[0],1)), np.array([0,1,2,3,4]))
#import pdb; pdb.set_trace()
#When we add the price of anarchy, the cost function becomes a0+ 2*a2*x+ 3*a3*x^2+ 4*a4*x^3+ 5*a5*x^4
#Matrix coefficients saves the 1, 2, 3, 4, and 5 coefficients in front of the ai*x terms
coefficients = np.array([1,2,3,4,5])
onesArray = np.ones((f.shape[0],5))
coefficients = onesArray * coefficients.transpose()
#y stands for the flow multiplied by the coefficients
#import pdb; pdb.set_trace()
y = np.einsum('ij, ij->ij', x, coefficients)
grad = np.einsum('ij,ij->i', y, graph[:,3:])
g.es["weight"] = grad.tolist()
#Calculating All_or_nothing
L = all_or_nothing(g, od)
return L, grad
def solver_social_optimum(graph, demand, g=None, od=None, past=10, max_iter=100, eps=1e-16, \
q=50, display=0, stop=1e-8):
'''
this is an adaptation of Fukushima's algorithm
graph: numpy array of the format [[link_id from to a0 a1 a2 a3 a4]]
demand: mumpy arrau of the format [[o d flow]]
g: igraph object constructed from graph
od: od in the format {from: ([to], [rate])}
past: search direction is the mean over the last 'past' directions
max_iter: maximum number of iterations
esp: used as a stopping criterium if some quantities are too close to 0
q: first 'q' iterations uses open loop step sizes 2/(i+2)
display: controls the display of information in the terminal
stop: stops the algorithm if the error is less than 'stop'
'''
assert past <= q, "'q' must be bigger or equal to 'past'"
if g is None:
g = construct_igraph(graph)
if od is None:
od = construct_od(demand)
f = np.zeros(graph.shape[0],dtype="float64") # initial flow assignment is null
fs = np.zeros((graph.shape[0],past),dtype="float64") #not sure what fs does
K = total_free_flow_cost(g, od)
# why this?
if K < eps:
K = np.sum(demand[:,2])
elif display >= 1:
print 'average free-flow travel time', K / np.sum(demand[:,2])
for i in range(max_iter):
if display >= 1:
if i <= 1:
print 'iteration: {}'.format(i+1)
else:
print 'iteration: {}, error: {}'.format(i+1, error)
# search direction using social optimum
L, grad = search_direction_social_optimum(f, graph, g, od)
fs[:,i%past] = L
w = L - f
if i >= 1:
error = -grad.dot(w) / K
# if error < stop and error > 0.0:
if error < stop:
if display >= 1: print 'stop with error: {}'.format(error)
return f
if i > q:
# step 3 of Fukushima
v = np.sum(fs,axis=1) / min(past,i+1) - f
norm_v = np.linalg.norm(v,1)
if norm_v < eps:
if display >= 1: print 'stop with norm_v: {}'.format(norm_v)
return f
norm_w = np.linalg.norm(w,1)
if norm_w < eps:
if display >= 1: print 'stop with norm_w: {}'.format(norm_w)
return f
# step 4 of Fukushima
gamma_1 = grad.dot(v) / norm_v
gamma_2 = grad.dot(w) / norm_w
if gamma_2 > -eps:
if display >= 1: print 'stop with gamma_2: {}'.format(gamma_2)
return f
d = v if gamma_1 < gamma_2 else w
# step 5 of Fukushima
s = line_search(lambda a: potential_Social_Optimum(graph, f+a*d))
if s < eps:
if display >= 1: print 'stop with step_size: {}'.format(s)
return f
f = f + s*d
else:
f = f + 2. * w/(i+2.)
return f