Hide keyboard shortcuts

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 

13 

14 

15def get_data(whereto): 

16 """ 

17 Récupère les données. 

18 

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) 

23 

24 

25def enumerate_events(df): 

26 """ 

27 Construit la liste des événements (vélo réposé ou retiré). 

28 

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,) 

50 

51 lastrow = row 

52 

53 

54class ParemetreCoutTrajet: 

55 

56 """ 

57 Regroupe l'ensembles des paramètres pour le calcul de la distance 

58 associé à un appariement. 

59 """ 

60 

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 

75 

76 def __str__(self): 

77 """ 

78 usuel 

79 """ 

80 return str(self.values) 

81 

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} 

88 

89 def cost(self, dh, dt, v): 

90 """ 

91 Retourne un coût, plus il est bas, 

92 plus de déplacement est probable. 

93 

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 

112 

113 

114def vitesse(c, d, params): 

115 """ 

116 Calcule la vitesse d'un déplacement. 

117 

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 

122 

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``. 

126 

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 

149 

150 

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. 

156 

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) 

176 

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) 

188 

189 

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. 

197 

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) 

206 

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] 

217 

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:] 

224 

225 default = (None, None, None, 0, 0, 0) 

226 

227 positif = [e for e in events if e[-1] > 0] 

228 negatif = [e for e in events if e[-1] < 0] 

229 

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) 

235 

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 

240 

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 

262 

263 moyenne = distance(positif, negatif, appariement_, params)[0] 

264 

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]]) 

272 

273 return mindist, moyenne, appariement_, positif, negatif 

274 

275 

276def distance_path(dfp): 

277 """ 

278 Calcule la vitesse moyenne lorsque le chemin est connu. 

279 

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 

292 

293 

294if __name__ == "__main__": # pragma: no cover 

295 

296 def main_velib(): 

297 dest = "codpart1" 

298 if not os.path.exists(dest): 

299 os.makedirs(dest) 

300 get_data(dest) 

301 

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") 

306 

307 # conversion des dates 

308 df["collect_date"] = df.apply( 

309 lambda r: str2datetime(r["collect_date"]), axis=1) 

310 

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) 

318 

319 # on calcule les événements (1 vélo apparu, 1 vélo disparu) 

320 events = list(sorted(enumerate_events(df))) 

321 

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) 

327 

328 main_velib()