X-Git-Url: https://www.fleuret.org/cgi-bin/gitweb/gitweb.cgi?p=mtp.git;a=blobdiff_plain;f=mtp_graph.cc;h=2cbb33848a3589007528a6e3b1098357dbf3d3df;hp=c17917ee675a3ed61eed111ba4965b9cfa0c0647;hb=HEAD;hpb=9d423de49725cbb204657653340f8585a0557f55

diff --git a/mtp_graph.cc b/mtp_graph.cc
index c17917e..2cbb338 100644
--- a/mtp_graph.cc
+++ b/mtp_graph.cc
@@ -97,36 +97,40 @@ void Vertex::del_leaving_edge(Edge *e) {
 
 void Vertex::decrease_distance_in_heap(Vertex **heap) {
   Vertex **p, **h;
-  // There is some beauty in that
   h = heap_slot;
-  while(h > heap &&
-        (p = heap + (h - heap + 1) / 2 - 1,
-         (*p)->distance_from_source > (*h)->distance_from_source)) {
+  while(1) {
+    if(h <= heap) break;
+    p = heap + ((h - heap + 1) >> 1) - 1;
+    if((*p)->distance_from_source <= distance_from_source) break;
+    swap((*p)->heap_slot, heap_slot);
     swap(*p, *h);
-    swap((*p)->heap_slot, (*h)->heap_slot);
     h = p;
   }
 }
 
 void Vertex::increase_distance_in_heap(Vertex **heap, Vertex **heap_bottom) {
   Vertex **c1, **c2, **h;
-  // omg, that's beautiful
   h = heap_slot;
-  while(c1 = heap + 2 * (h - heap) + 1,
-        c1 < heap_bottom &&
-        (c2 = c1 + 1,
-         (*c1)->distance_from_source < (*h)->distance_from_source
-         ||
-         (c2 < heap_bottom && (*c2)->distance_from_source < (*h)->distance_from_source)
-         )) {
-    if(c2 < heap_bottom && (*c2)->distance_from_source <= (*c1)->distance_from_source) {
-      swap(*c2, *h);
-      swap((*c2)->heap_slot, (*h)->heap_slot);
-      h = c2;
+  while(1) {
+    c1 = heap + 2 * (h - heap) + 1;
+    if(c1 >= heap_bottom) break;
+    c2 = c1 + 1;
+    if((*c1)->distance_from_source < distance_from_source) {
+      if(c2 < heap_bottom && (*c2)->distance_from_source < (*c1)->distance_from_source) {
+        swap((*c2)->heap_slot, heap_slot);
+        swap(*c2, *h);
+        h = c2;
+      } else {
+        swap((*c1)->heap_slot, heap_slot);
+        swap(*c1, *h);
+        h = c1;
+      }
     } else {
-      swap(*c1, *h);
-      swap((*c1)->heap_slot, (*h)->heap_slot);
-      h = c1;
+      if(c2 < heap_bottom && (*c2)->distance_from_source < distance_from_source) {
+        swap((*c2)->heap_slot, heap_slot);
+        swap(*c2, *h);
+        h = c2;
+      } else break;
     }
   }
 }
@@ -237,7 +241,7 @@ void MTPGraph::force_positivized_lengths() {
     }
   }
 #ifdef VERBOSE
-  cerr << __FILE__ << ": residual_error " << residual_error << " max_error " << residual_error << endl;
+  cerr << __FILE__ << ": residual_error " << residual_error << " max_error " << max_error << endl;
 #endif
 }
 
@@ -271,6 +275,7 @@ void MTPGraph::dp_compute_distances() {
 // pred_edge_toward_source.
 
 void MTPGraph::find_shortest_path() {
+  int heap_size;
   Vertex *v, *tv, **last_slot;
   Edge *e;
   scalar_t d;
@@ -280,20 +285,20 @@ void MTPGraph::find_shortest_path() {
     _vertices[k].pred_edge_toward_source = 0;
   }
 
-  _heap_size = _nb_vertices;
+  heap_size = _nb_vertices;
   _source->distance_from_source = 0;
   _source->decrease_distance_in_heap(_heap);
 
-  do {
+  while(heap_size > 1) {
     // Get the closest to the source
     v = _heap[0];
 
     // Remove it from the heap (swap it with the last_slot in the heap, and
     // update the distance of that one)
-    _heap_size--;
-    last_slot = _heap + _heap_size;
+    heap_size--;
+    last_slot = _heap + heap_size;
     swap(*_heap, *last_slot); swap((*_heap)->heap_slot, (*last_slot)->heap_slot);
-    _heap[0]->increase_distance_in_heap(_heap, _heap + _heap_size);
+    (*_heap)->increase_distance_in_heap(_heap, last_slot);
 
     // Now update the neighbors of the node currently closest to the
     // source
@@ -301,13 +306,13 @@ void MTPGraph::find_shortest_path() {
       d = v->distance_from_source + e->positivized_length;
       tv = e->terminal_vertex;
       if(d < tv->distance_from_source) {
-        ASSERT(tv->heap_slot - _heap < _heap_size);
+        ASSERT(tv->heap_slot < last_slot);
         tv->distance_from_source = d;
         tv->pred_edge_toward_source = e;
         tv->decrease_distance_in_heap(_heap);
       }
     }
-  } while(_heap_size > 0);
+  }
 }
 
 void MTPGraph::find_best_paths(scalar_t *lengths) {
@@ -374,7 +379,7 @@ void MTPGraph::find_best_paths(scalar_t *lengths) {
   }
 }
 
-int MTPGraph::retrieve_one_path(Edge *e, Path *path) {
+int MTPGraph::retrieve_one_path(Edge *e, Path *path, int *used_edges) {
   Edge *f, *next = 0;
   int l = 0, nb_occupied_next;
 
@@ -391,7 +396,9 @@ int MTPGraph::retrieve_one_path(Edge *e, Path *path) {
 
     nb_occupied_next = 0;
     for(f = e->terminal_vertex->leaving_edge_list_root; f; f = f->next_leaving_edge) {
-      if(f->occupied) { nb_occupied_next++; next = f; }
+      if(f->occupied && !used_edges[f - _edges]) {
+        nb_occupied_next++; next = f;
+      }
     }
 
 #ifdef DEBUG
@@ -399,13 +406,10 @@ int MTPGraph::retrieve_one_path(Edge *e, Path *path) {
       cerr << __FILE__ << ": retrieve_one_path: Non-sink end point." << endl;
       abort();
     }
-
-    else if(nb_occupied_next > 1) {
-      cerr << __FILE__ << ": retrieve_one_path: Non node-disjoint paths." << endl;
-      abort();
-    }
 #endif
 
+    if(path) { used_edges[next - _edges] = 1; }
+
     e = next;
   }
 
@@ -419,23 +423,15 @@ int MTPGraph::retrieve_one_path(Edge *e, Path *path) {
 
 //////////////////////////////////////////////////////////////////////
 
-static int compare_vertices_on_distance(const void *v1, const void *v2) {
-  scalar_t delta =
-    (*((Vertex **) v1))->distance_from_source -
-    (*((Vertex **) v2))->distance_from_source;
-  if(delta < 0) return -1;
-  else if(delta > 0) return 1;
-  else return 0;
-}
-
 void MTPGraph::compute_dp_ordering() {
   Vertex *v;
   Edge *e;
   int ntv;
 
-  // This method computes for each node the length of the longest link
-  // from the source, and orders the node in _dp_order according to
-  // it. It aborts if the graph is not a DAG.
+  // This method orders the nodes by putting first the ones with no
+  // predecessors, then going on adding nodes whose predecessors have
+  // all been already added. Computing the distances from the source
+  // by visiting nodes in that order is equivalent to DP.
 
   int *nb_predecessors = new int[_nb_vertices];
 
@@ -459,12 +455,15 @@ void MTPGraph::compute_dp_ordering() {
     }
   }
 
-  scalar_t rank = 1;
   while(already_processed < front) {
+    // Here, nodes before already_processed can be ignored, nodes
+    // before front were set to 0 predecessors during the previous
+    // iteration. During this new iteration, we have to visit the
+    // successors of these ones only, since they are the only ones
+    // which may end up with no predecessors.
     new_front = front;
     while(already_processed < front) {
       v = *(already_processed++);
-      v->distance_from_source = rank;
       for(e = v->leaving_edge_list_root; e; e = e->next_leaving_edge) {
         ntv = int(e->terminal_vertex - _vertices);
         nb_predecessors[ntv]--;
@@ -475,7 +474,6 @@ void MTPGraph::compute_dp_ordering() {
       }
     }
     front = new_front;
-    rank++;
   }
 
   if(already_processed < _dp_order + _nb_vertices) {
@@ -484,9 +482,6 @@ void MTPGraph::compute_dp_ordering() {
   }
 
   delete[] nb_predecessors;
-
-  for(int v = 0; v < _nb_vertices; v++) { _dp_order[v] = &_vertices[v]; }
-  qsort(_dp_order, _nb_vertices, sizeof(Vertex *), compare_vertices_on_distance);
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -494,6 +489,7 @@ void MTPGraph::compute_dp_ordering() {
 void MTPGraph::retrieve_disjoint_paths() {
   Edge *e;
   int p, l;
+  int *used_edges;
 
   for(int p = 0; p < nb_paths; p++) delete paths[p];
   delete[] paths;
@@ -504,14 +500,21 @@ void MTPGraph::retrieve_disjoint_paths() {
   }
 
   paths = new Path *[nb_paths];
+  used_edges = new int[_nb_edges];
+  for(int e = 0; e < _nb_edges; e++) {
+    used_edges[e] = 0;
+  }
 
   p = 0;
   for(e = _source->leaving_edge_list_root; e; e = e->next_leaving_edge) {
-    if(e->occupied) {
-      l = retrieve_one_path(e, 0);
+    if(e->occupied && !used_edges[e - _edges]) {
+      l = retrieve_one_path(e, 0, used_edges);
       paths[p] = new Path(l);
-      retrieve_one_path(e, paths[p]);
+      retrieve_one_path(e, paths[p], used_edges);
+      used_edges[e - _edges] = 1;
       p++;
     }
   }
+
+  delete[] used_edges;
 }