1 /****************************************************************************/
2 // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo
3 // Copyright (C) 2001-2019 German Aerospace Center (DLR) and others.
4 // This program and the accompanying materials
5 // are made available under the terms of the Eclipse Public License v2.0
6 // which accompanies this distribution, and is available at
7 // http://www.eclipse.org/legal/epl-v20.html
8 // SPDX-License-Identifier: EPL-2.0
9 /****************************************************************************/
10 /// @file    ROMAAssignments.cpp
11 /// @author  Yun-Pang Floetteroed
12 /// @author  Laura Bieker
13 /// @author  Michael Behrisch
14 /// @date    Feb 2013
15 /// @version $Id$
16 ///
17 // Assignment methods
18 /****************************************************************************/
19 
20 
21 // ===========================================================================
22 // included modules
23 // ===========================================================================
24 #include <config.h>
25 
26 #include <vector>
27 #include <algorithm>
28 #include <utils/common/SUMOTime.h>
29 #include <utils/distribution/Distribution_Points.h>
30 #include <utils/router/RouteCostCalculator.h>
31 #include <utils/router/SUMOAbstractRouter.h>
32 #include <router/ROEdge.h>
33 #include <router/RONet.h>
34 #include <router/RORoute.h>
35 #include <od/ODMatrix.h>
36 #include "ROMAEdge.h"
37 #include "ROMAAssignments.h"
38 
39 
40 // ===========================================================================
41 // static member variables
42 // ===========================================================================
43 std::map<const ROEdge* const, double> ROMAAssignments::myPenalties;
44 
45 
46 // ===========================================================================
47 // method definitions
48 // ===========================================================================
49 
ROMAAssignments(const SUMOTime begin,const SUMOTime end,const bool additiveTraffic,const double adaptionFactor,const int maxAlternatives,RONet & net,ODMatrix & matrix,SUMOAbstractRouter<ROEdge,ROVehicle> & router)50 ROMAAssignments::ROMAAssignments(const SUMOTime begin, const SUMOTime end, const bool additiveTraffic,
51                                  const double adaptionFactor, const int maxAlternatives, RONet& net, ODMatrix& matrix,
52                                  SUMOAbstractRouter<ROEdge, ROVehicle>& router)
53     : myBegin(begin), myEnd(end), myAdditiveTraffic(additiveTraffic), myAdaptionFactor(adaptionFactor),
54       myMaxAlternatives(maxAlternatives), myNet(net), myMatrix(matrix), myRouter(router) {
55     myDefaultVehicle = new ROVehicle(SUMOVehicleParameter(), nullptr, net.getVehicleTypeSecure(DEFAULT_VTYPE_ID), &net);
56 }
57 
58 
~ROMAAssignments()59 ROMAAssignments::~ROMAAssignments() {
60     delete myDefaultVehicle;
61 }
62 
63 // based on the definitions in PTV-Validate and in the VISUM-Cologne network
64 double
getCapacity(const ROEdge * edge)65 ROMAAssignments::getCapacity(const ROEdge* edge) {
66     if (edge->isTazConnector()) {
67         return 0;
68     }
69     const int roadClass = -edge->getPriority();
70     // TODO: differ road class 1 from the unknown road class 1!!!
71     if (edge->getNumLanes() == 0) {
72         // TAZ have no cost
73         return 0;
74     } else if (roadClass == 0 || roadClass == 1)  {
75         return edge->getNumLanes() * 2000.; //CR13 in table.py
76     } else if (roadClass == 2 && edge->getSpeedLimit() <= 11.) {
77         return edge->getNumLanes() * 1333.33; //CR5 in table.py
78     } else if (roadClass == 2 && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 16.) {
79         return edge->getNumLanes() * 1500.; //CR3 in table.py
80     } else if (roadClass == 2 && edge->getSpeedLimit() > 16.) {
81         return edge->getNumLanes() * 2000.; //CR13 in table.py
82     } else if (roadClass == 3 && edge->getSpeedLimit() <= 11.) {
83         return edge->getNumLanes() * 800.; //CR5 in table.py
84     } else if (roadClass == 3 && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 13.) {
85         return edge->getNumLanes() * 875.; //CR5 in table.py
86     } else if (roadClass == 3 && edge->getSpeedLimit() > 13. && edge->getSpeedLimit() <= 16.) {
87         return edge->getNumLanes() * 1500.; //CR4 in table.py
88     } else if (roadClass == 3 && edge->getSpeedLimit() > 16.) {
89         return edge->getNumLanes() * 1800.; //CR13 in table.py
90     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() <= 5.) {
91         return edge->getNumLanes() * 200.; //CR7 in table.py
92     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 5. && edge->getSpeedLimit() <= 7.) {
93         return edge->getNumLanes() * 412.5; //CR7 in table.py
94     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 7. && edge->getSpeedLimit() <= 9.) {
95         return edge->getNumLanes() * 600.; //CR6 in table.py
96     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 9. && edge->getSpeedLimit() <= 11.) {
97         return edge->getNumLanes() * 800.; //CR5 in table.py
98     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 13.) {
99         return edge->getNumLanes() * 1125.; //CR5 in table.py
100     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 13. && edge->getSpeedLimit() <= 16.) {
101         return edge->getNumLanes() * 1583.; //CR4 in table.py
102     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 16. && edge->getSpeedLimit() <= 18.) {
103         return edge->getNumLanes() * 1100.; //CR3 in table.py
104     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 18. && edge->getSpeedLimit() <= 22.) {
105         return edge->getNumLanes() * 1200.; //CR3 in table.py
106     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 22. && edge->getSpeedLimit() <= 26.) {
107         return edge->getNumLanes() * 1300.; //CR3 in table.py
108     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 26.) {
109         return edge->getNumLanes() * 1400.; //CR3 in table.py
110     }
111     return edge->getNumLanes() * 800.; //CR5 in table.py
112 }
113 
114 
115 // based on the definitions in PTV-Validate and in the VISUM-Cologne network
116 double
capacityConstraintFunction(const ROEdge * edge,const double flow) const117 ROMAAssignments::capacityConstraintFunction(const ROEdge* edge, const double flow) const {
118     if (edge->isTazConnector()) {
119         return 0;
120     }
121     const int roadClass = -edge->getPriority();
122     const double capacity = getCapacity(edge);
123     // TODO: differ road class 1 from the unknown road class 1!!!
124     if (edge->getNumLanes() == 0) {
125         // TAZ have no cost
126         return 0;
127     } else if (roadClass == 0 || roadClass == 1)  {
128         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.3)) * 2.); //CR13 in table.py
129     } else if (roadClass == 2 && edge->getSpeedLimit() <= 11.) {
130         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
131     } else if (roadClass == 2 && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 16.) {
132         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
133     } else if (roadClass == 2 && edge->getSpeedLimit() > 16.) {
134         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.3)) * 2.); //CR13 in table.py
135     } else if (roadClass == 3 && edge->getSpeedLimit() <= 11.) {
136         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
137     } else if (roadClass == 3 && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 13.) {
138         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
139     } else if (roadClass == 3 && edge->getSpeedLimit() > 13. && edge->getSpeedLimit() <= 16.) {
140         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.7 * (flow / (capacity * 1.)) * 2.); //CR4 in table.py
141     } else if (roadClass == 3 && edge->getSpeedLimit() > 16.) {
142         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.3)) * 2.); //CR13 in table.py
143     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() <= 5.) {
144         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.5)) * 3.); //CR7 in table.py
145     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 5. && edge->getSpeedLimit() <= 7.) {
146         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.5)) * 3.); //CR7 in table.py
147     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 7. && edge->getSpeedLimit() <= 9.) {
148         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.8)) * 3.); //CR6 in table.py
149     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 9. && edge->getSpeedLimit() <= 11.) {
150         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
151     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 13.) {
152         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
153     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 13. && edge->getSpeedLimit() <= 16.) {
154         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.7 * (flow / (capacity * 1.)) * 2.); //CR4 in table.py
155     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 16. && edge->getSpeedLimit() <= 18.) {
156         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
157     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 18. && edge->getSpeedLimit() <= 22.) {
158         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
159     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 22. && edge->getSpeedLimit() <= 26.) {
160         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
161     } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 26.) {
162         return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
163     }
164     return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
165 }
166 
167 
168 bool
addRoute(const ConstROEdgeVector & edges,std::vector<RORoute * > & paths,std::string routeId,double prob)169 ROMAAssignments::addRoute(const ConstROEdgeVector& edges, std::vector<RORoute*>& paths, std::string routeId, double prob) {
170     std::vector<RORoute*>::iterator p;
171     for (p = paths.begin(); p != paths.end(); p++) {
172         if (edges == (*p)->getEdgeVector()) {
173             break;
174         }
175     }
176     if (p == paths.end()) {
177         paths.push_back(new RORoute(routeId, 0., prob, edges, nullptr, std::vector<SUMOVehicleParameter::Stop>()));
178         return true;
179     }
180     (*p)->addProbability(prob);
181     std::iter_swap(paths.end() - 1, p);
182     return false;
183 }
184 
185 
186 const ConstROEdgeVector
computePath(ODCell * cell,const SUMOTime time,const double probability,SUMOAbstractRouter<ROEdge,ROVehicle> * router)187 ROMAAssignments::computePath(ODCell* cell, const SUMOTime time, const double probability, SUMOAbstractRouter<ROEdge, ROVehicle>* router) {
188     const ROEdge* const from = myNet.getEdge(cell->origin + (cell->originIsEdge ? "" : "-source"));
189     if (from == nullptr) {
190         throw ProcessError("Unknown origin '" + cell->origin + "'.");
191     }
192     const ROEdge* const to = myNet.getEdge(cell->destination + (cell->destinationIsEdge ? "" : "-sink"));
193     if (to == nullptr) {
194         throw ProcessError("Unknown destination '" + cell->destination + "'.");
195     }
196     ConstROEdgeVector edges;
197     if (router == nullptr) {
198         router = &myRouter;
199     }
200     if (myMaxAlternatives > 0 && (int)cell->pathsVector.size() < myMaxAlternatives) {
201         router->compute(from, to, myDefaultVehicle, time, edges);
202         if (addRoute(edges, cell->pathsVector, cell->origin + cell->destination + toString(cell->pathsVector.size()), probability)) {
203             return edges;
204         }
205     } else {
206         double minCost = std::numeric_limits<double>::max();
207         RORoute* minRoute = nullptr;
208         for (RORoute* const p : cell->pathsVector) {
209             const double cost = router->recomputeCosts(edges, myDefaultVehicle, time);
210             if (cost < minCost) {
211                 minCost = cost;
212                 minRoute = p;
213             }
214         }
215         minRoute->addProbability(probability);
216     }
217     return ConstROEdgeVector();
218 }
219 
220 
221 void
getKPaths(const int kPaths,const double penalty)222 ROMAAssignments::getKPaths(const int kPaths, const double penalty) {
223     for (ODCell* const c : myMatrix.getCells()) {
224         myPenalties.clear();
225         for (int k = 0; k < kPaths; k++) {
226             for (const ROEdge* const e : computePath(c)) {
227                 myPenalties[e] += penalty;
228             }
229         }
230     }
231     myPenalties.clear();
232 }
233 
234 
235 void
resetFlows()236 ROMAAssignments::resetFlows() {
237     const double begin = STEPS2TIME(MIN2(myBegin, myMatrix.getCells().front()->begin));
238     for (std::map<std::string, ROEdge*>::const_iterator i = myNet.getEdgeMap().begin(); i != myNet.getEdgeMap().end(); ++i) {
239         ROMAEdge* edge = static_cast<ROMAEdge*>(i->second);
240         edge->setFlow(begin, STEPS2TIME(myEnd), 0.);
241         edge->setHelpFlow(begin, STEPS2TIME(myEnd), 0.);
242     }
243 }
244 
245 
246 void
incremental(const int numIter,const bool verbose)247 ROMAAssignments::incremental(const int numIter, const bool verbose) {
248     SUMOTime lastBegin = -1;
249     std::vector<int> intervals;
250     int count = 0;
251     for (const ODCell* const c : myMatrix.getCells()) {
252         if (c->begin != lastBegin) {
253             intervals.push_back(count);
254             lastBegin = c->begin;
255         }
256         count++;
257     }
258     lastBegin = -1;
259     for (std::vector<int>::const_iterator offset = intervals.begin(); offset != intervals.end(); offset++) {
260         std::vector<ODCell*>::const_iterator cellsEnd = myMatrix.getCells().end();
261         if (offset != intervals.end() - 1) {
262             cellsEnd = myMatrix.getCells().begin() + (*(offset + 1));
263         }
264         const SUMOTime intervalStart = (*(myMatrix.getCells().begin() + (*offset)))->begin;
265         if (verbose) {
266             WRITE_MESSAGE(" starting interval " + time2string(intervalStart));
267         }
268         std::map<const ROMAEdge*, double> loadedTravelTimes;
269         for (std::map<std::string, ROEdge*>::const_iterator i = myNet.getEdgeMap().begin(); i != myNet.getEdgeMap().end(); ++i) {
270             ROMAEdge* edge = static_cast<ROMAEdge*>(i->second);
271             if (edge->hasLoadedTravelTime(STEPS2TIME(intervalStart))) {
272                 loadedTravelTimes[edge] = edge->getTravelTime(myDefaultVehicle, STEPS2TIME(intervalStart));
273             }
274         }
275         for (int t = 0; t < numIter; t++) {
276             if (verbose) {
277                 WRITE_MESSAGE("  starting iteration " + toString(t));
278             }
279             std::string lastOrigin = "";
280             int workerIndex = 0;
281             for (std::vector<ODCell*>::const_iterator i = myMatrix.getCells().begin() + (*offset); i != cellsEnd; i++) {
282                 ODCell* const c = *i;
283                 const double linkFlow = c->vehicleNumber / numIter;
284                 const SUMOTime begin = myAdditiveTraffic ? myBegin : c->begin;
285 #ifdef HAVE_FOX
286                 if (myNet.getThreadPool().size() > 0) {
287                     if (lastOrigin != c->origin) {
288                         workerIndex++;
289                         if (workerIndex == myNet.getThreadPool().size()) {
290                             workerIndex = 0;
291                         }
292                         myNet.getThreadPool().add(new RONet::BulkmodeTask(false), workerIndex);
293                         lastOrigin = c->origin;
294                         myNet.getThreadPool().add(new RoutingTask(*this, c, begin, linkFlow), workerIndex);
295                         myNet.getThreadPool().add(new RONet::BulkmodeTask(true), workerIndex);
296                     } else {
297                         myNet.getThreadPool().add(new RoutingTask(*this, c, begin, linkFlow), workerIndex);
298                     }
299                     continue;
300                 }
301 #endif
302                 if (lastOrigin != c->origin) {
303                     myRouter.setBulkMode(false);
304                     lastOrigin = c->origin;
305                 }
306                 computePath(c, begin, linkFlow);
307                 myRouter.setBulkMode(true);
308             }
309 #ifdef HAVE_FOX
310             if (myNet.getThreadPool().size() > 0) {
311                 myNet.getThreadPool().waitAll();
312             }
313 #endif
314             for (std::vector<ODCell*>::const_iterator i = myMatrix.getCells().begin() + (*offset); i != cellsEnd; i++) {
315                 ODCell* const c = *i;
316                 const double linkFlow = c->vehicleNumber / numIter;
317                 const SUMOTime begin = myAdditiveTraffic ? myBegin : c->begin;
318                 const SUMOTime end = myAdditiveTraffic ? myEnd : c->end;
319                 const double intervalLengthInHours = STEPS2TIME(end - begin) / 3600.;
320                 const ConstROEdgeVector& edges = c->pathsVector.back()->getEdgeVector();
321                 for (ConstROEdgeVector::const_iterator e = edges.begin(); e != edges.end(); e++) {
322                     ROMAEdge* edge = static_cast<ROMAEdge*>(myNet.getEdge((*e)->getID()));
323                     const double newFlow = edge->getFlow(STEPS2TIME(begin)) + linkFlow;
324                     edge->setFlow(STEPS2TIME(begin), STEPS2TIME(end), newFlow);
325                     double travelTime = capacityConstraintFunction(edge, newFlow / intervalLengthInHours);
326                     if (lastBegin >= 0 && myAdaptionFactor > 0.) {
327                         if (loadedTravelTimes.count(edge) != 0) {
328                             travelTime = loadedTravelTimes[edge] * myAdaptionFactor + (1. - myAdaptionFactor) * travelTime;
329                         } else {
330                             travelTime = edge->getTravelTime(myDefaultVehicle, STEPS2TIME(lastBegin)) * myAdaptionFactor + (1. - myAdaptionFactor) * travelTime;
331                         }
332                     }
333                     edge->addTravelTime(travelTime, STEPS2TIME(begin), STEPS2TIME(end));
334                 }
335             }
336         }
337         lastBegin = intervalStart;
338     }
339 }
340 
341 
342 void
sue(const int maxOuterIteration,const int maxInnerIteration,const int kPaths,const double penalty,const double tolerance,const std::string)343 ROMAAssignments::sue(const int maxOuterIteration, const int maxInnerIteration, const int kPaths, const double penalty, const double tolerance, const std::string /* routeChoiceMethod */) {
344     getKPaths(kPaths, penalty);
345     std::map<const double, double> intervals;
346     if (myAdditiveTraffic) {
347         intervals[STEPS2TIME(myBegin)] = STEPS2TIME(myEnd);
348     } else {
349         for (const ODCell* const c : myMatrix.getCells()) {
350             intervals[STEPS2TIME(c->begin)] = STEPS2TIME(c->end);
351         }
352     }
353     for (int outer = 0; outer < maxOuterIteration; outer++) {
354         for (int inner = 0; inner < maxInnerIteration; inner++) {
355             for (const ODCell* const c : myMatrix.getCells()) {
356                 const SUMOTime begin = myAdditiveTraffic ? myBegin : c->begin;
357                 const SUMOTime end = myAdditiveTraffic ? myEnd : c->end;
358                 // update path cost
359                 for (std::vector<RORoute*>::const_iterator j = c->pathsVector.begin(); j != c->pathsVector.end(); ++j) {
360                     RORoute* r = *j;
361                     r->setCosts(myRouter.recomputeCosts(r->getEdgeVector(), myDefaultVehicle, 0));
362 //                    std::cout << std::setprecision(20) << r->getID() << ":" << r->getCosts() << std::endl;
363                 }
364                 // calculate route utilities and probabilities
365                 RouteCostCalculator<RORoute, ROEdge, ROVehicle>::getCalculator().calculateProbabilities(c->pathsVector, myDefaultVehicle, 0);
366                 // calculate route flows
367                 for (std::vector<RORoute*>::const_iterator j = c->pathsVector.begin(); j != c->pathsVector.end(); ++j) {
368                     RORoute* r = *j;
369                     const double pathFlow = r->getProbability() * c->vehicleNumber;
370                     // assign edge flow deltas
371                     for (ConstROEdgeVector::const_iterator e = r->getEdgeVector().begin(); e != r->getEdgeVector().end(); e++) {
372                         ROMAEdge* edge = static_cast<ROMAEdge*>(myNet.getEdge((*e)->getID()));
373                         edge->setHelpFlow(STEPS2TIME(begin), STEPS2TIME(end), edge->getHelpFlow(STEPS2TIME(begin)) + pathFlow);
374                     }
375                 }
376             }
377             // calculate new edge flows and check for stability
378             int unstableEdges = 0;
379             for (std::map<const double, double>::const_iterator i = intervals.begin(); i != intervals.end(); ++i) {
380                 const double intervalLengthInHours = STEPS2TIME(i->second - i->first) / 3600.;
381                 for (std::map<std::string, ROEdge*>::const_iterator e = myNet.getEdgeMap().begin(); e != myNet.getEdgeMap().end(); ++e) {
382                     ROMAEdge* edge = static_cast<ROMAEdge*>(e->second);
383                     const double oldFlow = edge->getFlow(i->first);
384                     double newFlow = oldFlow;
385                     if (inner == 0 && outer == 0) {
386                         newFlow += edge->getHelpFlow(i->first);
387                     } else {
388                         newFlow += (edge->getHelpFlow(i->first) - oldFlow) / (inner + 1);
389                     }
390                     //                if not lohse:
391                     if (newFlow > 0.) {
392                         if (fabs(newFlow - oldFlow) / newFlow > tolerance) {
393                             unstableEdges++;
394                         }
395                     } else if (newFlow == 0.) {
396                         if (oldFlow != 0. && (fabs(newFlow - oldFlow) / oldFlow > tolerance)) {
397                             unstableEdges++;
398                         }
399                     } else { // newFlow < 0.
400                         unstableEdges++;
401                         newFlow = 0.;
402                     }
403                     edge->setFlow(i->first, i->second, newFlow);
404                     const double travelTime = capacityConstraintFunction(edge, newFlow / intervalLengthInHours);
405                     edge->addTravelTime(travelTime, i->first, i->second);
406                     edge->setHelpFlow(i->first, i->second, 0.);
407                 }
408             }
409             // if stable break
410             if (unstableEdges == 0) {
411                 break;
412             }
413             // additional stability check from python script: if notstable < math.ceil(net.geteffEdgeCounts()*0.005) or notstable < 3: stable = True
414         }
415         // check for a new route, if none available, break
416         // several modifications about when a route is new and when to break are in the original script
417         bool newRoute = false;
418         for (ODCell* const c : myMatrix.getCells()) {
419             newRoute |= !computePath(c).empty();
420         }
421         if (!newRoute) {
422             break;
423         }
424     }
425     // final round of assignment
426     for (const ODCell* const c : myMatrix.getCells()) {
427         // update path cost
428         for (std::vector<RORoute*>::const_iterator j = c->pathsVector.begin(); j != c->pathsVector.end(); ++j) {
429             RORoute* r = *j;
430             r->setCosts(myRouter.recomputeCosts(r->getEdgeVector(), myDefaultVehicle, 0));
431         }
432         // calculate route utilities and probabilities
433         RouteCostCalculator<RORoute, ROEdge, ROVehicle>::getCalculator().calculateProbabilities(c->pathsVector, myDefaultVehicle, 0);
434         // calculate route flows
435         for (std::vector<RORoute*>::const_iterator j = c->pathsVector.begin(); j != c->pathsVector.end(); ++j) {
436             RORoute* r = *j;
437             r->setProbability(r->getProbability() * c->vehicleNumber);
438         }
439     }
440 }
441 
442 
443 double
getPenalizedEffort(const ROEdge * const e,const ROVehicle * const v,double t)444 ROMAAssignments::getPenalizedEffort(const ROEdge* const e, const ROVehicle* const v, double t) {
445     const std::map<const ROEdge* const, double>::const_iterator i = myPenalties.find(e);
446     return i == myPenalties.end() ? e->getEffort(v, t) : e->getEffort(v, t) + i->second;
447 }
448 
449 
450 double
getPenalizedTT(const ROEdge * const e,const ROVehicle * const v,double t)451 ROMAAssignments::getPenalizedTT(const ROEdge* const e, const ROVehicle* const v, double t) {
452     const std::map<const ROEdge* const, double>::const_iterator i = myPenalties.find(e);
453     return i == myPenalties.end() ? e->getTravelTime(v, t) : e->getTravelTime(v, t) + i->second;
454 }
455 
456 
457 double
getTravelTime(const ROEdge * const e,const ROVehicle * const v,double t)458 ROMAAssignments::getTravelTime(const ROEdge* const e, const ROVehicle* const v, double t) {
459     return e->getTravelTime(v, t);
460 }
461 
462 
463 #ifdef HAVE_FOX
464 // ---------------------------------------------------------------------------
465 // ROMAAssignments::RoutingTask-methods
466 // ---------------------------------------------------------------------------
467 void
run(FXWorkerThread * context)468 ROMAAssignments::RoutingTask::run(FXWorkerThread* context) {
469     myAssign.computePath(myCell, myBegin, myLinkFlow, &static_cast<RONet::WorkerThread*>(context)->getVehicleRouter());
470 }
471 #endif
472