SUMO - Simulation of Urban MObility
NBAlgorithms.cpp
Go to the documentation of this file.
1 /****************************************************************************/
2 // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.org/sumo
3 // Copyright (C) 2012-2017 German Aerospace Center (DLR) and others.
4 /****************************************************************************/
5 //
6 // This program and the accompanying materials
7 // are made available under the terms of the Eclipse Public License v2.0
8 // which accompanies this distribution, and is available at
9 // http://www.eclipse.org/legal/epl-v20.html
10 //
11 /****************************************************************************/
18 // Algorithms for network computation
19 /****************************************************************************/
20 
21 
22 // ===========================================================================
23 // included modules
24 // ===========================================================================
25 #ifdef _MSC_VER
26 #include <windows_config.h>
27 #else
28 #include <config.h>
29 #endif
30 
31 #include <sstream>
32 #include <iostream>
33 #include <cassert>
34 #include <algorithm>
36 #include <utils/common/ToString.h>
37 #include "NBEdge.h"
38 #include "NBNodeCont.h"
39 #include "NBTypeCont.h"
40 #include "NBNode.h"
41 #include "NBAlgorithms.h"
42 
43 
44 // ===========================================================================
45 // method definitions
46 // ===========================================================================
47 // ---------------------------------------------------------------------------
48 // NBTurningDirectionsComputer
49 // ---------------------------------------------------------------------------
50 void
52  for (std::map<std::string, NBNode*>::const_iterator i = nc.begin(); i != nc.end(); ++i) {
53  computeTurnDirectionsForNode(i->second, warn);
54  }
55 }
56 
57 void
59  const std::vector<NBEdge*>& incoming = node->getIncomingEdges();
60  const std::vector<NBEdge*>& outgoing = node->getOutgoingEdges();
61  // reset turning directions since this may be called multiple times
62  for (std::vector<NBEdge*>::const_iterator k = incoming.begin(); k != incoming.end(); ++k) {
63  (*k)->setTurningDestination(0);
64  }
65  std::vector<Combination> combinations;
66  for (std::vector<NBEdge*>::const_iterator j = outgoing.begin(); j != outgoing.end(); ++j) {
67  NBEdge* outedge = *j;
68  for (std::vector<NBEdge*>::const_iterator k = incoming.begin(); k != incoming.end(); ++k) {
69  NBEdge* e = *k;
70  // @todo: check whether NBHelpers::relAngle is properly defined and whether it should really be used, here
71  const double signedAngle = NBHelpers::normRelAngle(e->getAngleAtNode(node), outedge->getAngleAtNode(node));
72  if (signedAngle > 0 && signedAngle < 177 && e->getGeometry().back().distanceTo2D(outedge->getGeometry().front()) < POSITION_EPS) {
73  // backwards curving edges can only be turnaround when there are
74  // non-default endpoints
75  continue;
76  }
77  double angle = fabs(signedAngle);
78  // std::cout << "incoming=" << e->getID() << " outgoing=" << outedge->getID() << " relAngle=" << NBHelpers::relAngle(e->getAngleAtNode(node), outedge->getAngleAtNode(node)) << "\n";
79  if (e->getFromNode() == outedge->getToNode() && angle > 120) {
80  // they connect the same nodes; should be the turnaround direction
81  // we'll assign a maximum number
82  //
83  // @todo: indeed, we have observed some pathological intersections
84  // see "294831560" in OSM/adlershof. Here, several edges are connecting
85  // same nodes. We have to do the angle check before...
86  //
87  // @todo: and well, there are some other as well, see plain import
88  // of delphi_muenchen (elmar), intersection "59534191". Not that it would
89  // be realistic in any means; we will warn, here.
90  angle += 360;
91  }
92  if (angle < 160) {
93  continue;
94  }
95  Combination c;
96  c.from = e;
97  c.to = outedge;
98  c.angle = angle;
99  combinations.push_back(c);
100  }
101  }
102  // sort combinations so that the ones with the highest angle are at the begin
103  std::sort(combinations.begin(), combinations.end(), combination_by_angle_sorter());
104  std::set<NBEdge*> seen;
105  for (std::vector<Combination>::const_iterator j = combinations.begin(); j != combinations.end(); ++j) {
106  if (seen.find((*j).from) != seen.end() || seen.find((*j).to) != seen.end()) {
107  // do not regard already set edges
108  if ((*j).angle > 360 && warn) {
109  WRITE_WARNING("Ambiguity in turnarounds computation at junction '" + node->getID() + "'.");
110  warn = false;
111  }
112  continue;
113  }
114  // mark as seen
115  seen.insert((*j).from);
116  seen.insert((*j).to);
117  // set turnaround information
118  bool onlyPossible = (*j).from->getConnections().size() != 0 && !(*j).from->isConnectedTo((*j).to);
119  //std::cout << " setTurningDestination from=" << (*j).from->getID() << " to=" << (*j).to->getID() << " onlyPossible=" << onlyPossible << "\n";
120  (*j).from->setTurningDestination((*j).to, onlyPossible);
121  }
122 }
123 
124 
125 // ---------------------------------------------------------------------------
126 // NBNodesEdgesSorter
127 // ---------------------------------------------------------------------------
128 void
130  for (std::map<std::string, NBNode*>::const_iterator i = nc.begin(); i != nc.end(); ++i) {
131  i->second->sortEdges(useNodeShape);
132  }
133 }
134 
135 
136 void
138  const std::vector<NBEdge*>::iterator& i1,
139  const std::vector<NBEdge*>::iterator& i2) {
140  NBEdge* e1 = *i1;
141  NBEdge* e2 = *i2;
142  // @todo: The difference between "isTurningDirectionAt" and "isTurnaround"
143  // is not nice. Maybe we could get rid of it if we would always mark edges
144  // as turnarounds, even if they do not have to be added, as mentioned in
145  // notes on NBTurningDirectionsComputer::computeTurnDirectionsForNode
146  if (e2->getToNode() == n && e2->isTurningDirectionAt(e1)) {
147  std::swap(*i1, *i2);
148  }
149 }
150 
151 
152 // ---------------------------------------------------------------------------
153 // NBNodeTypeComputer
154 // ---------------------------------------------------------------------------
155 void
157  for (std::map<std::string, NBNode*>::const_iterator i = nc.begin(); i != nc.end(); ++i) {
158  NBNode* n = (*i).second;
159  // the type may already be set from the data
160  if (n->myType != NODETYPE_UNKNOWN && n->myType != NODETYPE_DEAD_END) {
161  continue;
162  }
163  // check whether the node is a waterway node. Set to unregulated by default
164  bool waterway = true;
165  for (EdgeVector::const_iterator i = n->getEdges().begin(); i != n->getEdges().end(); ++i) {
166  if (!isWaterway((*i)->getPermissions())) {
167  waterway = false;
168  break;
169  }
170  }
171  if (waterway && (n->myType == NODETYPE_UNKNOWN || n->myType == NODETYPE_DEAD_END)) {
173  continue;
174  }
175 
176  // check whether the junction is not a real junction
177  if (n->myIncomingEdges.size() == 1) {
179  continue;
180  }
181  // @todo "isSimpleContinuation" should be revalidated
182  if (n->isSimpleContinuation()) {
184  continue;
185  }
186  // determine the type
188  for (EdgeVector::const_iterator i = n->myIncomingEdges.begin(); i != n->myIncomingEdges.end(); i++) {
189  for (EdgeVector::const_iterator j = i + 1; j != n->myIncomingEdges.end(); j++) {
190  // @todo "getOppositeIncoming" should probably be refactored into something the edge knows
191  if (n->getOppositeIncoming(*j) == *i && n->myIncomingEdges.size() > 2) {
192  continue;
193  }
194  // @todo check against a legal document
195  // @todo figure out when NODETYPE_PRIORITY_STOP is appropriate
196  const double s1 = (*i)->getSpeed() * (double) 3.6;
197  const double s2 = (*j)->getSpeed() * (double) 3.6;
198  const int p1 = (*i)->getPriority();
199  const int p2 = (*j)->getPriority();
200  if (fabs(s1 - s2) > (double) 9.5 || MAX2(s1, s2) >= (double) 49. || p1 != p2) {
201  type = NODETYPE_PRIORITY;
202  break;
203  }
204  }
205  }
206  // save type
207  n->myType = type;
208  }
209 }
210 
211 
212 void
214  // the type may already be set from the data
215  if (node->myType != NODETYPE_UNKNOWN && node->myType != NODETYPE_DEAD_END) {
216  }
217  // check whether the node is a waterway node. Set to unregulated by default
218  bool waterway = true;
219  for (EdgeVector::const_iterator i = node->getEdges().begin(); i != node->getEdges().end(); ++i) {
220  if (!isWaterway((*i)->getPermissions())) {
221  waterway = false;
222  break;
223  }
224  }
225  if (waterway && (node->myType == NODETYPE_UNKNOWN || node->myType == NODETYPE_DEAD_END)) {
226  node->myType = NODETYPE_NOJUNCTION;
227  }
228 
229  // check whether the junction is not a real junction
230  if (node->myIncomingEdges.size() == 1) {
231  node->myType = NODETYPE_PRIORITY;
232  }
233  // @todo "isSimpleContinuation" should be revalidated
234  if (node->isSimpleContinuation()) {
235  node->myType = NODETYPE_PRIORITY;
236  }
237  // determine the type
239  for (EdgeVector::const_iterator i = node->myIncomingEdges.begin(); i != node->myIncomingEdges.end(); i++) {
240  for (EdgeVector::const_iterator j = i + 1; j != node->myIncomingEdges.end(); j++) {
241  // @todo "getOppositeIncoming" should probably be refactored into something the edge knows
242  if (node->getOppositeIncoming(*j) == *i && node->myIncomingEdges.size() > 2) {
243  continue;
244  }
245  // @todo check against a legal document
246  // @todo figure out when NODETYPE_PRIORITY_STOP is appropriate
247  const double s1 = (*i)->getSpeed() * (double) 3.6;
248  const double s2 = (*j)->getSpeed() * (double) 3.6;
249  const int p1 = (*i)->getPriority();
250  const int p2 = (*j)->getPriority();
251  if (fabs(s1 - s2) > (double) 9.5 || MAX2(s1, s2) >= (double) 49. || p1 != p2) {
252  type = NODETYPE_PRIORITY;
253  break;
254  }
255  }
256  }
257  // save type
258  node->myType = type;
259 }
260 
261 // ---------------------------------------------------------------------------
262 // NBEdgePriorityComputer
263 // ---------------------------------------------------------------------------
264 void
266  for (std::map<std::string, NBNode*>::const_iterator i = nc.begin(); i != nc.end(); ++i) {
267  computeEdgePrioritiesSingleNode((*i).second);
268  }
269 }
270 
271 
272 void
274  // preset all junction's edge priorities to zero
275  for (EdgeVector::iterator j = node->myAllEdges.begin(); j != node->myAllEdges.end(); ++j) {
276  (*j)->setJunctionPriority(node, NBEdge::MINOR_ROAD);
277  }
278  // check if the junction is not a real junction
279  if (node->myIncomingEdges.size() == 1 && node->myOutgoingEdges.size() == 1) {
280  return;
281  }
282  // compute the priorities on junction when needed
283  if (node->getType() != NODETYPE_RIGHT_BEFORE_LEFT) {
284  setPriorityJunctionPriorities(*node);
285  }
286 }
287 
288 
289 void
291  if (n.myIncomingEdges.size() == 0 || n.myOutgoingEdges.size() == 0) {
292  return;
293  }
294  EdgeVector incoming = n.myIncomingEdges;
295  EdgeVector outgoing = n.myOutgoingEdges;
296  // what we do want to have is to extract the pair of roads that are
297  // the major roads for this junction
298  // let's get the list of incoming edges with the highest priority
299  std::sort(incoming.begin(), incoming.end(), NBContHelper::edge_by_priority_sorter());
300  EdgeVector bestIncoming;
301  NBEdge* best = incoming[0];
302  while (incoming.size() > 0 && samePriority(best, incoming[0])) {
303  bestIncoming.push_back(*incoming.begin());
304  incoming.erase(incoming.begin());
305  }
306  // now, let's get the list of best outgoing
307  assert(outgoing.size() != 0);
308  sort(outgoing.begin(), outgoing.end(), NBContHelper::edge_by_priority_sorter());
309  EdgeVector bestOutgoing;
310  best = outgoing[0];
311  while (outgoing.size() > 0 && samePriority(best, outgoing[0])) { //->getPriority()==best->getPriority()) {
312  bestOutgoing.push_back(*outgoing.begin());
313  outgoing.erase(outgoing.begin());
314  }
315  // special case: user input makes mainDirection unambiguous
316  const bool mainDirectionExplicit = (
317  bestIncoming.size() == 1 && n.myIncomingEdges.size() <= 2
318  && (incoming.size() == 0 || bestIncoming[0]->getPriority() > incoming[0]->getPriority())
319  && bestOutgoing.size() == 1 && n.myOutgoingEdges.size() <= 2
320  && (outgoing.size() == 0 || bestOutgoing[0]->getPriority() > outgoing[0]->getPriority())
321  && !bestIncoming[0]->isTurningDirectionAt(bestOutgoing[0]));
322  // now, let's compute for each of the best incoming edges
323  // the incoming which is most opposite
324  // the outgoing which is most opposite
325  EdgeVector::iterator i;
326  std::map<NBEdge*, NBEdge*> counterIncomingEdges;
327  std::map<NBEdge*, NBEdge*> counterOutgoingEdges;
328  incoming = n.myIncomingEdges;
329  outgoing = n.myOutgoingEdges;
330  for (i = bestIncoming.begin(); i != bestIncoming.end(); ++i) {
331  std::sort(incoming.begin(), incoming.end(), NBContHelper::edge_opposite_direction_sorter(*i, &n, true));
332  counterIncomingEdges[*i] = *incoming.begin();
333  std::sort(outgoing.begin(), outgoing.end(), NBContHelper::edge_opposite_direction_sorter(*i, &n, true));
334  counterOutgoingEdges[*i] = *outgoing.begin();
335  }
336  //std::cout << "n=" << n.getID() << " best=" << best->getID() << " bestIncoming=" << toString(bestIncoming) << "\n incoming=" << toString(incoming) << "\n outgoing=" << toString(outgoing) << "\n mainExplicit=" << mainDirectionExplicit << " counterBest=" << counterIncomingEdges.find(bestIncoming[0])->second->getID() << "\n";
337  // ok, let's try
338  // 1) there is one best incoming road
339  if (bestIncoming.size() == 1) {
340  // let's mark this road as the best
341  NBEdge* best1 = extractAndMarkFirst(n, bestIncoming);
342  if (!mainDirectionExplicit && counterIncomingEdges.find(best1) != counterIncomingEdges.end()) {
343  // ok, look, what we want is the opposit of the straight continuation edge
344  // but, what if such an edge does not exist? By now, we'll determine it
345  // geometrically
346  NBEdge* s = counterIncomingEdges.find(best1)->second;
347  const double minAngleDiff = GeomHelper::getMinAngleDiff(best1->getAngleAtNode(&n), s->getAngleAtNode(&n));
348  if (minAngleDiff > 180 - 45
349  || (minAngleDiff > 75 && s->getPriority() == best1->getPriority() && hasDifferentPriorities(incoming, best1))) {
351  }
352  }
353  assert(bestOutgoing.size() != 0);
354  // mark the best outgoing as the continuation
355  sort(bestOutgoing.begin(), bestOutgoing.end(), NBContHelper::edge_similar_direction_sorter(best1));
356  // assign extra priority if the priorities are unambiguous (regardless of geometry)
357  best1 = extractAndMarkFirst(n, bestOutgoing);
358  if (!mainDirectionExplicit && counterOutgoingEdges.find(best1) != counterOutgoingEdges.end()) {
359  NBEdge* s = counterOutgoingEdges.find(best1)->second;
360  if (GeomHelper::getMinAngleDiff(best1->getAngleAtNode(&n), s->getAngleAtNode(&n)) > 180 - 45) {
361  s->setJunctionPriority(&n, 1);
362  }
363  }
364  return;
365  }
366 
367  // ok, what we want to do in this case is to determine which incoming
368  // has the best continuation...
369  // This means, when several incoming roads have the same priority,
370  // we want a (any) straight connection to be more priorised than a turning
371  double bestAngle = 0;
372  NBEdge* bestFirst = 0;
373  NBEdge* bestSecond = 0;
374  bool hadBest = false;
375  for (i = bestIncoming.begin(); i != bestIncoming.end(); ++i) {
376  EdgeVector::iterator j;
377  NBEdge* t1 = *i;
378  double angle1 = t1->getAngleAtNode(&n) + 180;
379  if (angle1 >= 360) {
380  angle1 -= 360;
381  }
382  for (j = i + 1; j != bestIncoming.end(); ++j) {
383  NBEdge* t2 = *j;
384  double angle2 = t2->getAngleAtNode(&n) + 180;
385  if (angle2 >= 360) {
386  angle2 -= 360;
387  }
388  double angle = GeomHelper::getMinAngleDiff(angle1, angle2);
389  if (!hadBest || angle > bestAngle) {
390  bestAngle = angle;
391  bestFirst = *i;
392  bestSecond = *j;
393  hadBest = true;
394  }
395  }
396  }
397  bestFirst->setJunctionPriority(&n, 1);
398  sort(bestOutgoing.begin(), bestOutgoing.end(), NBContHelper::edge_similar_direction_sorter(bestFirst));
399  if (bestOutgoing.size() != 0) {
400  extractAndMarkFirst(n, bestOutgoing);
401  }
402  bestSecond->setJunctionPriority(&n, 1);
403  sort(bestOutgoing.begin(), bestOutgoing.end(), NBContHelper::edge_similar_direction_sorter(bestSecond));
404  if (bestOutgoing.size() != 0) {
405  extractAndMarkFirst(n, bestOutgoing);
406  }
407 }
408 
409 
410 NBEdge*
411 NBEdgePriorityComputer::extractAndMarkFirst(NBNode& n, std::vector<NBEdge*>& s, int prio) {
412  if (s.size() == 0) {
413  return 0;
414  }
415  NBEdge* ret = s.front();
416  s.erase(s.begin());
417  ret->setJunctionPriority(&n, prio);
418  return ret;
419 }
420 
421 
422 bool
423 NBEdgePriorityComputer::samePriority(const NBEdge* const e1, const NBEdge* const e2) {
424  if (e1 == e2) {
425  return true;
426  }
427  if (e1->getPriority() != e2->getPriority()) {
428  return false;
429  }
430  if ((int) e1->getSpeed() != (int) e2->getSpeed()) {
431  return false;
432  }
433  return (int) e1->getNumLanes() == (int) e2->getNumLanes();
434 }
435 
436 
437 bool
439  if (edges.size() < 2) {
440  return false;
441  }
442  int prio = edges[0] == excluded ? edges[1]->getPriority() : edges[0]->getPriority();
443  for (auto e : edges) {
444  if (e != excluded && e->getPriority() != prio) {
445  return true;
446  }
447  }
448  return false;
449 }
450 
451 
453  // reorder based on getAngleAtNodeToCenter
454  myOrdering = ordering;
455  sort(myOrdering.begin(), myOrdering.end(), NBContHelper::edge_by_angle_to_nodeShapeCentroid_sorter(node));
456  // let the first edge remain the first
457  rotate(myOrdering.begin(), std::find(myOrdering.begin(), myOrdering.end(), ordering.front()), myOrdering.end());
458 }
459 
460 
461 /****************************************************************************/
462 
static double getMinAngleDiff(double angle1, double angle2)
Returns the minimum distance (clockwise/counter-clockwise) between both angles.
Definition: GeomHelper.cpp:167
std::map< std::string, NBNode * >::const_iterator begin() const
Returns the pointer to the begin of the stored nodes.
Definition: NBNodeCont.h:117
std::map< std::string, NBNode * >::const_iterator end() const
Returns the pointer to the end of the stored nodes.
Definition: NBNodeCont.h:122
static double normRelAngle(double angle1, double angle2)
ensure that reverse relAngles (>=179.999) always count as turnarounds (-180)
Definition: NBHelpers.cpp:66
SumoXMLNodeType myType
The type of the junction.
Definition: NBNode.h:746
int getPriority() const
Returns the priority of the edge.
Definition: NBEdge.h:419
The representation of a single edge during network building.
Definition: NBEdge.h:70
Class to sort edges by their angle in relation to the given edge.
Definition: NBContHelper.h:177
T MAX2(T a, T b)
Definition: StdDefs.h:73
static void computeSingleNodeType(NBNode *node)
Computes a single node type.
const std::string & getID() const
Returns the id.
Definition: Named.h:65
static void sortNodesEdges(NBNodeCont &nc, bool useNodeShape=false)
Sorts a node&#39;s edges clockwise regarding driving direction.
Stores the information about the angle between an incoming ("from") and an outgoing ("to") edge...
Definition: NBAlgorithms.h:74
#define WRITE_WARNING(msg)
Definition: MsgHandler.h:199
const EdgeVector & getOutgoingEdges() const
Returns this node&#39;s outgoing edges (The edges which start at this node)
Definition: NBNode.h:254
static bool hasDifferentPriorities(const EdgeVector &edges, const NBEdge *excluded)
return whether the priorite attribute can be used to distinguish the edges
EdgeVector myAllEdges
Vector of incoming and outgoing edges.
Definition: NBNode.h:734
static void swapWhenReversed(const NBNode *const n, const std::vector< NBEdge *>::iterator &i1, const std::vector< NBEdge *>::iterator &i2)
Assures correct order for same-angle opposite-direction edges.
int getNumLanes() const
Returns the number of lanes.
Definition: NBEdge.h:412
static void computeEdgePriorities(NBNodeCont &nc)
Computes edge priorities within a node.
bool isWaterway(SVCPermissions permissions)
Returns whether an edge with the given permission is a waterway edge.
const EdgeVector & getEdges() const
Returns all edges which participate in this node (Edges that start or end at this node) ...
Definition: NBNode.h:259
#define POSITION_EPS
Definition: config.h:175
double getAngleAtNode(const NBNode *const node) const
Returns the angle of the edge&#39;s geometry at the given node.
Definition: NBEdge.cpp:1611
crossing_by_junction_angle_sorter(const NBNode *node, const EdgeVector &ordering)
EdgeVector myIncomingEdges
Vector of incoming edges.
Definition: NBNode.h:728
static void computeTurnDirectionsForNode(NBNode *node, bool warn)
Computes turnaround destinations for all incoming edges of the given nodes (if any) ...
static bool samePriority(const NBEdge *const e1, const NBEdge *const e2)
Returns whether both edges have the same priority.
double getSpeed() const
Returns the speed allowed on this edge.
Definition: NBEdge.h:506
const PositionVector & getGeometry() const
Returns the geometry of the edge.
Definition: NBEdge.h:602
EdgeVector myOutgoingEdges
Vector of outgoing edges.
Definition: NBNode.h:731
static void computeEdgePrioritiesSingleNode(NBNode *node)
Computes edge priorities within a single node.
const EdgeVector & getIncomingEdges() const
Returns this node&#39;s incoming edges (The edges which yield in this node)
Definition: NBNode.h:249
static void computeNodeTypes(NBNodeCont &nc)
Computes node types.
SumoXMLNodeType
Numbers representing special SUMO-XML-attribute values for representing node- (junction-) types used ...
std::vector< NBEdge * > EdgeVector
container for (sorted) edges
Definition: NBCont.h:40
void setJunctionPriority(const NBNode *const node, int prio)
Sets the junction priority of the edge.
Definition: NBEdge.cpp:1601
SumoXMLNodeType getType() const
Returns the type of this node.
Definition: NBNode.h:266
bool isTurningDirectionAt(const NBEdge *const edge) const
Returns whether the given edge is the opposite direction to this edge.
Definition: NBEdge.cpp:2450
Represents a single node (junction) during network building.
Definition: NBNode.h:74
static void computeTurnDirections(NBNodeCont &nc, bool warn=true)
Computes turnaround destinations for all edges (if exist)
static NBEdge * extractAndMarkFirst(NBNode &n, std::vector< NBEdge *> &s, int prio=1)
Sets the priorites in case of a priority junction.
NBEdge * getOppositeIncoming(NBEdge *e) const
returns the opposite incoming edge of certain edge
Definition: NBNode.cpp:1220
Sorts "Combination"s by decreasing angle.
Definition: NBAlgorithms.h:84
static void setPriorityJunctionPriorities(NBNode &n)
Sets the priorites in case of a priority junction.
NBNode * getFromNode() const
Returns the origin node of the edge.
Definition: NBEdge.h:426
Container for nodes during the netbuilding process.
Definition: NBNodeCont.h:66
NBNode * getToNode() const
Returns the destination node of the edge.
Definition: NBEdge.h:433
bool isSimpleContinuation(bool checkLaneNumbers=true) const
check if node is a simple continuation
Definition: NBNode.cpp:431