Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2"""
3@file
4@brief Une solution au problème proposée :
5:ref:`Reconstruction de trajectoire velib <l-codingparty1>`
6"""
7import os
8import random
9import pandas
10from pyensae.datasource import download_data
11from manydataapi.velib import DataCollectJCDecaux as DataVelibCollect
12from pyquickhelper.loghelper import str2datetime
15def get_data(whereto):
16 """
17 Récupère les données.
19 @param whereto destination
20 """
21 download_data('velib_synthetique.zip', website='xdtd', whereTo=whereto)
22 download_data('besancon.df.txt.zip', website='xdtd', whereTo=whereto)
25def enumerate_events(df):
26 """
27 Construit la liste des événements (vélo réposé ou retiré).
29 @param df DataFrame
30 @return énumère événéments (ierator) ("file", "collect_date", "name", "lat", "lng", +1 ou -1)
31 """
32 df = df[["file", "collect_date", "name", "lat", "lng",
33 "available_bike_stands", "available_bikes"]]
34 df = df.sort_values(["name", "file", "collect_date"])
35 lastrow = None
36 for row in df.values:
37 if lastrow is not None and lastrow[2] == row[2]: # pylint: disable=E1136
38 # nombre de place en plus
39 d1 = row[-2] - lastrow[-2] # pylint: disable=E1136
40 # nombre de vélos en plus
41 d2 = row[-1] - lastrow[-1] # pylint: disable=E1136
42 if d1 != 0:
43 step = d1 // abs(d1)
44 for _ in range(1, abs(d1) + 1):
45 yield tuple(row[:-2]) + (step,)
46 elif d2 != 0:
47 step = d2 // abs(d2)
48 for _ in range(1, abs(d2) + 1):
49 yield tuple(row[:-2]) + (step,)
51 lastrow = row
54class ParemetreCoutTrajet:
56 """
57 Regroupe l'ensembles des paramètres pour le calcul de la distance
58 associé à un appariement.
59 """
61 def __init__(self, max_speed=50, high_speed=25, low_speed=5,
62 high_time=0.75, low_time=0.1):
63 """
64 @param max_speed au-delà, c'est une vitesse impossible
65 @param high_speed au-delà, c'est la vitesse d'un cycliste
66 @param low_speed en-deça, il vaut mieux marcher
67 @param high_time au-delà, c'est une randonnée
68 @param low_time en-deça, pourquoi un vélo
69 """
70 self.max_speed = max_speed
71 self.high_speed = high_speed
72 self.low_speed = low_speed
73 self.high_time = high_time
74 self.low_time = low_time
76 def __str__(self):
77 """
78 usuel
79 """
80 return str(self.values)
82 @property
83 def values(self):
84 """
85 Retourne les valeurs dans un dictionnaire.
86 """
87 return {k: getattr(self, k) for k in self.__dict__ if "time" in k or "speed" in k}
89 def cost(self, dh, dt, v):
90 """
91 Retourne un coût, plus il est bas,
92 plus de déplacement est probable.
94 @param dh distance
95 @param dt durée
96 @param v vitesse
97 @return coût
98 """
99 c = 0
100 if v > self.max_speed:
101 c += 1e3
102 if v > self.high_speed:
103 c += (v - self.high_speed) ** 2 * 10
104 if v < self.low_speed:
105 c += (v - self.low_speed) ** 2 * 100
106 dt = dt.total_seconds() / 3600
107 if dt > self.high_time:
108 c += (dt - self.high_time) ** 2 * 10
109 if dt < self.low_time:
110 c += (dt - self.low_time) ** 2 * 10
111 return c
114def vitesse(c, d, params):
115 """
116 Calcule la vitesse d'un déplacement.
118 @param c tuple ("file", "collect_date", "name", "lat", "lng", +1 ou -1)
119 @param d tuple ("file", "collect_date", "name", "lat", "lng", +1 ou -1)
120 @param params ParemetreCoutTrajet
121 @return vitesse, cout
123 La fonction retourne une valeur aberrante si le temps entre les deux événements est négatifs.
124 C'est une configuration impossible : on ne peut reposer un vélo avant de l'avoir retiré.
125 La valeur aberrante est ``1e8``.
127 Il reste un cas pour lequel, je ne sais pas encore quelle valeur donner :
128 il s'agit des demi-appariements : un vélo rétiré mais jamais reposé et réciproquement.
129 """
130 if c[0] is None or d[0] is None:
131 # cas des vélos perdus
132 if c[0] is None:
133 if d[0] is None:
134 return None
135 else:
136 return 0.0, 0.0 # je ne sais pas trop quoi mettre
137 else:
138 return 0.0, 0.0 # je ne sais pas trop quoi mettre
139 else:
140 lat1, lng1 = c[3], c[4]
141 lat2, lng2 = d[3], d[4]
142 dh = DataVelibCollect.distance_haversine(lat1, lng1, lat2, lng2)
143 dt = d[1] - c[1]
144 if dt.total_seconds() <= 0:
145 return 1e8, 1e8 # infini
146 v = dh / (dt.total_seconds() / 3600)
147 cost = params.cost(dh, dt, v)
148 return v, cost
151def distance(positif, negatif, app, params):
152 """
153 Calcule une distance pour un appariement conçu ici comme
154 la variance de la vitesse de chaque déplacement + la somme
155 des coûts de chaque appariement retourné par la fonction @see fn vitesse.
157 @param positif vélos pris (ou l'inverse)
158 @param negatif vélos remis (ou l'inverse)
159 @param app appariement (list de tuple (i,j))
160 @param params @see cl ParemetreCoutTrajet
161 @return tuple: vitesse moyenne (sans appariements négatifs),
162 distance, vitesse moyenne avec, nombre d'appariements négatifs
163 """
164 val = []
165 nb_max = []
166 cost = 0
167 for i, j in app:
168 p = positif[i]
169 n = negatif[j]
170 v, d = vitesse(p, n, params)
171 cost += d
172 if v > params.max_speed:
173 nb_max.append(v)
174 if v is not None:
175 val.append(v)
177 mean = sum(val) / len(val)
178 # on enlève les appariements négatifs
179 cor = [v_ for v_ in val if v_ < 1e8]
180 if len(cor) == 0:
181 raise RuntimeError("Pas d'appariements positifs: {0}".format(val))
182 mean_cor = sum(cor) / len(cor)
183 dev = sum((x - mean) ** 2 for x in val) / len(val)
184 return mean_cor, \
185 cost + dev ** 0.5 * (1 + len(nb_max)), \
186 mean, \
187 len(val) - len(cor)
190def appariement(events, iter=1000, params=ParemetreCoutTrajet(), fLOG=print): # pylint: disable=W0622
191 """
192 On veut apparier les événemens -1 aux événemens +1.
193 On s'attend aux colonnes suivantes:
194 *"file"*, *"collect_date"*, *"name"*, *"lat"*, *"lng"*, *+1* ou *-1*.
195 La fonction part des deux séries d'évéments rangés par ordre croissant,
196 puis il permute deux appariements tant que cela diminue l'erreur d'appariement.
198 @param events list d'événements produits par la fonction
199 @see fn enumerate_events
200 @param iter nombre d'itérations
201 @param params ParemetreCoutTrajet
202 @param fLOG logging function
203 @return tuple (mindist, moyenne, appariement, positif, negatif)
204 """
205 events = sorted(events)
207 # on élimine tous les -1 du début (forcément des vélos pris avant la
208 # période d'étude)
209 lasti = 0
210 for i, ev in enumerate(events):
211 lasti = i
212 if ev[-1] == 1:
213 break
214 i = lasti
215 if i > 0:
216 del events[:i]
218 # pareil de l'autre côté
219 for i, ev in enumerate(reversed(events)):
220 if ev[-1] == -1:
221 break
222 if i > 0:
223 del events[len(events) - i:]
225 default = (None, None, None, 0, 0, 0)
227 positif = [e for e in events if e[-1] > 0]
228 negatif = [e for e in events if e[-1] < 0]
230 # on veut autant d'événements de chaque côté
231 while len(positif) > len(negatif):
232 negatif.append(default)
233 while len(positif) < len(negatif):
234 positif.append(default)
236 appariement_ = [(i, i) for i in range(0, len(positif))]
237 vit, mindist, vitav, nbneg = distance( # pylint: disable=W0612
238 positif, negatif, appariement_, params)
239 nbchange = 0
241 for it in range(0, iter):
242 if it % 10 == 0:
243 fLOG("iteration ", it, ": app-", nbneg, "/", len(appariement_),
244 "min", mindist, "vitesse ", vit, " nbchange", nbchange)
245 nbchange = 0
246 for _ in range(0, len(appariement_)):
247 i = random.randint(0, len(appariement_) - 1)
248 j = random.randint(0, len(appariement_) - 1)
249 if i == j:
250 continue
251 ki, kj = appariement_[i], appariement_[j]
252 appariement_[i] = (ki[0], kj[1])
253 appariement_[j] = (kj[0], ki[1])
254 v, dist, vt, nbneg = distance( # pylint: disable=W0612
255 positif, negatif, appariement_, params)
256 if dist < mindist:
257 mindist = dist
258 vit = v
259 nbchange += 1
260 else:
261 appariement_[i], appariement_[j] = ki, kj
263 moyenne = distance(positif, negatif, appariement_, params)[0]
265 # def dd(a, b):
266 # try:
267 # return b - a
268 # except Exception:
269 # return None
270 # for a in appariement_:
271 # fLOG(a,dd(positif [a[0]][1],negatif[a[1]][1]), vitesse(positif [a[0]], negatif[a[1]], params), positif [a[0]],"-->",negatif[a[1]])
273 return mindist, moyenne, appariement_, positif, negatif
276def distance_path(dfp):
277 """
278 Calcule la vitesse moyenne lorsque le chemin est connu.
280 @param dfp liste des chemins
281 @return moyenne, stddev
282 """
283 dfp = dfp.copy()
284 dfp["speed"] = dfp["dist"] / dfp["hours"]
285 dfp["dist"] = dfp.apply(
286 lambda r: DataVelibCollect.distance_haversine(r["lat0"], r["lng0"],
287 r["lat1"], r["lng1"]),
288 axis=1)
289 mean_ = sum(dfp["speed"]) / len(dfp)
290 std_ = sum((x - mean_) ** 2 for x in dfp["speed"]) / len(dfp)
291 return mean_, std_ ** 0.5
294if __name__ == "__main__": # pragma: no cover
296 def main_velib():
297 dest = "codpart1"
298 if not os.path.exists(dest):
299 os.makedirs(dest)
300 get_data(dest)
302 # récupère les données
303 jeu = os.path.join(dest, "besancon.df.txt")
304 jeu = os.path.join(dest, "out_simul_bike_nb1_sp10_data.txt")
305 df = pandas.read_csv(jeu, sep="\t", encoding="utf8")
307 # conversion des dates
308 df["collect_date"] = df.apply(
309 lambda r: str2datetime(r["collect_date"]), axis=1)
311 # on regarde s'il existe le fichier des trajectoires
312 path = jeu.replace("_data.", "_path.")
313 if path != jeu and os.path.exists(path):
314 dfp = pandas.read_csv(path, sep="\t")
315 dfp = dfp[dfp["beginend"] == "end"]
316 mean, std = distance_path(dfp)
317 print("expected: vitesse moyenne ", mean, " stddev ", std)
319 # on calcule les événements (1 vélo apparu, 1 vélo disparu)
320 events = list(sorted(enumerate_events(df)))
322 params = ParemetreCoutTrajet()
323 print(params)
324 mindist, moyenne, appariement_, positif, negatif = appariement( # pylint: disable=W0612
325 events, iter=200, params=params)
326 print("vitesse moyenne", moyenne)
328 main_velib()