From 2b3a2e10ec226f1610b9c39abd20f0899a34a652 Mon Sep 17 00:00:00 2001 From: Francois Fleuret Date: Wed, 19 Dec 2012 14:23:13 +0100 Subject: [PATCH] Seems to work. Lot of debugging garbage still in there. --- mtp_graph.cc | 265 +++++++++++++++++++++++++++++++++++++++++---------- mtp_graph.h | 20 +++- 2 files changed, 233 insertions(+), 52 deletions(-) diff --git a/mtp_graph.cc b/mtp_graph.cc index 5371732..12017c2 100644 --- a/mtp_graph.cc +++ b/mtp_graph.cc @@ -50,6 +50,8 @@ public: int last_change; // Used to mark which edges have already been // processed in some methods + Vertex **heap_position; + Vertex(); inline void add_leaving_edge(Edge *e); @@ -95,6 +97,10 @@ void Vertex::del_leaving_edge(Edge *e) { ////////////////////////////////////////////////////////////////////// +static int compare_vertex(const void *v1, const void *v2) { + return (*((Vertex **) v1))->last_change - (*((Vertex **) v2))->last_change; +} + MTPGraph::MTPGraph(int nb_vertices, int nb_edges, int *vertex_from, int *vertex_to, int source, int sink) { @@ -103,8 +109,8 @@ MTPGraph::MTPGraph(int nb_vertices, int nb_edges, _edges = new Edge[_nb_edges]; _vertices = new Vertex[_nb_vertices]; - _front = new Vertex *[_nb_vertices]; - _new_front = new Vertex *[_nb_vertices]; + _heap = new Vertex *[_nb_vertices]; + _dp_order = new Vertex *[_nb_vertices]; _source = &_vertices[source]; _sink = &_vertices[sink]; @@ -116,15 +122,31 @@ MTPGraph::MTPGraph(int nb_vertices, int nb_edges, _edges[e].terminal_vertex = _vertices + vertex_to[e]; } + for(int v = 0; v < _nb_vertices; v++) { + _heap[v] = &_vertices[v]; + _vertices[v].heap_position = &_heap[v]; + } + paths = 0; nb_paths = 0; + + if(is_dag()) { + // Here the last_change field of every vertex tells us how many + // iterations of DP we need to reach it. Hence we only have to + // process the vertex in that order. + for(int v = 0; v < _nb_vertices; v++) { _dp_order[v] = &_vertices[v]; } + qsort(_dp_order, _nb_vertices, sizeof(Vertex *), compare_vertex); + } else { + cerr << __FILE__ << ": This graph is not a DAG." << endl; + abort(); + } } MTPGraph::~MTPGraph() { delete[] _vertices; + delete[] _dp_order; + delete[] _heap; delete[] _edges; - delete[] _front; - delete[] _new_front; for(int p = 0; p < nb_paths; p++) delete paths[p]; delete[] paths; } @@ -213,40 +235,170 @@ int MTPGraph::is_dag() { Vertex *v; Edge *e; - // We put everybody in the front + Vertex **active = new Vertex *[_nb_vertices]; + + // We put everybody in the active for(int k = 0; k < _nb_vertices; k++) { - _vertices[k].last_change = -1; - _front[k] = &_vertices[k]; + _vertices[k].last_change = 0; + active[k] = &_vertices[k]; } - int iteration = 0; - int front_size = _nb_vertices, pred_front_size; + int iteration = 1; + int nb_active = _nb_vertices, pred_nb_active; do { // We set the last_change field of all the vertices with incoming // edges to the current iteration value - for(int f = 0; f < front_size; f++) { - v = _front[f]; + for(int f = 0; f < nb_active; f++) { + v = active[f]; for(e = v->leaving_edges; e; e = e->next_leaving_edge) { e->terminal_vertex->last_change = iteration; } } - pred_front_size = front_size; - front_size = 0; + pred_nb_active = nb_active; + nb_active = 0; // We keep all the vertices with incoming nodes - for(int f = 0; f < pred_front_size; f++) { - v = _front[f]; + for(int f = 0; f < pred_nb_active; f++) { + v = active[f]; if(v->last_change == iteration) { - _front[front_size++] = v; + active[nb_active++] = v; } } iteration++; - } while(front_size < pred_front_size); + } while(nb_active < pred_nb_active); + + delete[] active; + + return nb_active == 0; +} + +#ifdef DEBUG +void print_heap(int first, int heap_size, Vertex **heap, int depth, Vertex *vertices) { + // if(depth == 0) cout << "** START_HEAP ********************************************************" << endl; + // if(first < heap_size) { + // print_heap((first + 1) * 2 - 1, heap_size, heap, depth + 1, vertices); + // for(int d = 0; d < depth; d++) cout << " "; + // cout << "[" << heap[first] - vertices << "] "; + // if(heap[first]->distance_from_source == FLT_MAX) cout << "inf"; + // else cout << heap[first]->distance_from_source; + // cout << endl; + // print_heap((first + 1) * 2, heap_size, heap, depth + 1, vertices); + // } + // if(depth == 0) cout << "** END_HEAP **********************************************************" << endl; +} + +void check_heap(Vertex **heap, int heap_size, Vertex *vertices, int nb_vertices) { + Vertex **p, **h; + int *used = new int[nb_vertices]; + for(int v = 0; v < nb_vertices; v++) used[v] = 0; + for(int k = 0; k < heap_size; k++) { + used[heap[k] - vertices]++; + h = heap + k; + if(h != (*h)->heap_position) abort(); + if(k > 0) { + p = heap + (h - heap + 1) / 2 - 1; + if((*h)->distance_from_source < (*p)->distance_from_source) abort(); + } + } + for(int v = 0; v < nb_vertices; v++) { + cout << used[v]; + if(used[v] > 1) abort(); + } + cout << endl; + delete[] used; +} +#endif - return front_size == 0; +void MTPGraph::decrease_distance_in_heap(Vertex *v) { +#ifdef DEBUG + // cout << "START decrease_distance_in_heap" << endl; + // print_heap(0, _heap_size, _heap, 0, _vertices); +#endif + Vertex **p, **h; + // There is some beauty in that + h = v->heap_position; + while(h > _heap && + (p = _heap + (h - _heap + 1) / 2 - 1, + (*p)->distance_from_source > (*h)->distance_from_source)) { +// #warning REMOVE + // cout << "SWAP [" << (*p) - _vertices << " | " << (*h) - _vertices << "]" << endl; + // if((*p) - _vertices == 6 && (*h) - _vertices == 96) abort(); + swap(*p, *h); + swap((*p)->heap_position, (*h)->heap_position); + h = p; + } +#ifdef DEBUG + // check_heap(_heap, _heap_size, _vertices, _nb_vertices); + // print_heap(0, _heap_size, _heap, 0, _vertices); +#endif +} + +void MTPGraph::increase_distance_in_heap(Vertex *v) { +#ifdef DEBUG + // cout << "START increase_distance_in_heap" << endl; + // print_heap(0, _heap_size, _heap, 0, _vertices); +#endif + Vertex **c1, **c2, **h; + // There is some beauty in that + h = v->heap_position; + while(c1 = _heap + 2 * (h - _heap + 1) - 1, c2 = c1 + 1, + (c1 < _heap + _heap_size && + (*c1)->distance_from_source < (*h)->distance_from_source) + || + (c2 < _heap + _heap_size && + (*c2)->distance_from_source < (*h)->distance_from_source) + ) { + if(c1 < _heap + _heap_size && + !(c2 < _heap + _heap_size && + (*c2)->distance_from_source < (*c1)->distance_from_source)){ + swap(*c1, *h); + swap((*c1)->heap_position, (*h)->heap_position); + h = c1; + } else { + swap(*c2, *h); + swap((*c2)->heap_position, (*h)->heap_position); + h = c2; + } + } +#ifdef DEBUG + // check_heap(_heap, _heap_size, _vertices, _nb_vertices); + // print_heap(0, _heap_size, _heap, 0, _vertices); +#endif +} + +void MTPGraph::dp_distance_propagation() { + Vertex *v, *tv; + Edge *e; + scalar_t d; + + for(int k = 0; k < _nb_vertices; k++) { + _vertices[k].distance_from_source = FLT_MAX; + _vertices[k].pred_edge_toward_source = 0; + } + + _source->distance_from_source = 0; + + for(int k = 0; k < _nb_vertices; k++) { + v = _dp_order[k]; + for(e = v->leaving_edges; e; e = e->next_leaving_edge) { + d = v->distance_from_source + e->positivized_length; + tv = e->terminal_vertex; + if(d < tv->distance_from_source) { + tv->distance_from_source = d; + tv->pred_edge_toward_source = e; + decrease_distance_in_heap(tv); + } + } + } + +#ifdef DEBUG + for(int k = 0; k < _nb_vertices; k++) { + if(_vertices[k].distance_from_source == FLT_MAX) abort(); + } +#endif } // This method does not change the edge occupation. It only set @@ -254,7 +406,6 @@ int MTPGraph::is_dag() { // pred_edge_toward_source. void MTPGraph::find_shortest_path() { - Vertex **tmp_front; Vertex *v, *tv; Edge *e; scalar_t d; @@ -264,7 +415,7 @@ void MTPGraph::find_shortest_path() { cout << "find_shortest_path: DAG -> ok" << endl; } else { for(int e = 0; e < _nb_edges; e++) { - if(_edges[e].positivized_length < 0) abort(); + if(_edges[e].positivized_length < 0) abort(); } cout << "find_shortest_path: All positivized_length are positive -> ok" << endl; } @@ -273,40 +424,54 @@ void MTPGraph::find_shortest_path() { for(int k = 0; k < _nb_vertices; k++) { _vertices[k].distance_from_source = FLT_MAX; _vertices[k].pred_edge_toward_source = 0; - _vertices[k].last_change = -1; } - int iteration = 0; - - int front_size = 0, new_front_size; - _front[front_size++] = _source; + _heap_size = _nb_vertices; _source->distance_from_source = 0; + decrease_distance_in_heap(_source); do { - new_front_size = 0; +#ifdef DEBUG + for(int k = 0; k < _heap_size; k++) { + if(_heap[0]->distance_from_source > _heap[k]->distance_from_source) abort(); + } + // cout << "OK!" << endl; +#endif - for(int f = 0; f < front_size; f++) { - v = _front[f]; - for(e = v->leaving_edges; e; e = e->next_leaving_edge) { - d = v->distance_from_source + e->positivized_length; - tv = e->terminal_vertex; - if(d < tv->distance_from_source) { - tv->distance_from_source = d; - tv->pred_edge_toward_source = e; - if(tv->last_change < iteration) { - _new_front[new_front_size++] = tv; - tv->last_change = iteration; - } - } + v = _heap[0]; + _heap_size--; + Vertex **a = _heap, **b = _heap + _heap_size; + swap(*a, *b); swap((*a)->heap_position, (*b)->heap_position); + increase_distance_in_heap(_heap[0]); + + for(e = v->leaving_edges; e; e = e->next_leaving_edge) { + d = v->distance_from_source + e->positivized_length; + tv = e->terminal_vertex; + if(d < tv->distance_from_source) { + ASSERT(tv->heap_position - _heap < _heap_size); + tv->distance_from_source = d; + tv->pred_edge_toward_source = e; + decrease_distance_in_heap(tv); + } +#ifdef DEBUG + for(int k = 0; k < _heap_size; k++) { + if(_heap[0]->distance_from_source > _heap[k]->distance_from_source) abort(); } +#endif } + } while(_heap_size > 0); - tmp_front = _new_front; _new_front = _front; _front = tmp_front; - - front_size = new_front_size; - - iteration++; - } while(front_size > 0); +#ifdef DEBUG + for(int k = 0; k < _nb_vertices; k++) { + v = &_vertices[k]; + for(e = v->leaving_edges; e; e = e->next_leaving_edge) { + d = v->distance_from_source + e->positivized_length; + tv = e->terminal_vertex; + if(d < tv->distance_from_source) abort(); + } + } +#endif + cout << "DONE!" << endl; } void MTPGraph::find_best_paths(scalar_t *lengths) { @@ -320,10 +485,12 @@ void MTPGraph::find_best_paths(scalar_t *lengths) { _edges[e].positivized_length = _edges[e].length; } - // We call find_shortest_path here to set properly the distances to - // the source, so that we can make all the edge lengths positive at - // the first iteration. - find_shortest_path(); + // // We call find_shortest_path here to set properly the distances to + // // the source, so that we can make all the edge lengths positive at + // // the first iteration. + // find_shortest_path(); + + dp_distance_propagation(); do { update_positivized_lengths(); diff --git a/mtp_graph.h b/mtp_graph.h index f28e460..5842671 100644 --- a/mtp_graph.h +++ b/mtp_graph.h @@ -49,6 +49,13 @@ class MTPGraph { // the total correction when compiled in VERBOSE mode. void force_positivized_lengths(); + void decrease_distance_in_heap(Vertex *v); + void increase_distance_in_heap(Vertex *v); + + // Visit the vertices according to _dp_order and simply update their + // distance to the source + void dp_distance_propagation(); + // Set in every vertex pred_edge_toward_source correspondingly to // the path of shortest length. The current implementation is not // Dijkstra's! @@ -59,17 +66,24 @@ class MTPGraph { // nodes met along the path, and computes path->length properly. int retrieve_one_path(Edge *e, Path *path); - // Returns if the graph is a DAG + // Returns if the graph is a DAG, and set the last_change field of + // each vertex to the maximum number of iterations required to reach + // it int is_dag(); - Vertex **_front, **_new_front; - int _nb_vertices, _nb_edges; Vertex *_source, *_sink; Edge *_edges; Vertex *_vertices; + // For the shortest path search + + Vertex **_heap; + int _heap_size; + + Vertex **_dp_order; + public: // These variables are filled when retrieve_disjoint_paths is called -- 2.20.1