diff options
Diffstat (limited to 'blossom5-v2.05.src/GEOM/GPMmain.cpp')
| -rw-r--r-- | blossom5-v2.05.src/GEOM/GPMmain.cpp | 148 |
1 files changed, 148 insertions, 0 deletions
diff --git a/blossom5-v2.05.src/GEOM/GPMmain.cpp b/blossom5-v2.05.src/GEOM/GPMmain.cpp new file mode 100644 index 0000000..d2bcb3c --- /dev/null +++ b/blossom5-v2.05.src/GEOM/GPMmain.cpp @@ -0,0 +1,148 @@ +#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GeomPerfectMatching.h"
+#include "GPMkdtree.h"
+#include "../timer.h"
+
+
+
+GeomPerfectMatching::REAL GeomPerfectMatching::SolveComplete()
+{
+ if (node_num != node_num_max) { printf("ComputeCost() cannot be called before all points have been added!\n"); exit(1); }
+
+ PointId p, q;
+ int e = 0, E = node_num*(node_num-1)/2;
+ PerfectMatching* pm = new PerfectMatching(node_num, E);
+ for (p=0; p<node_num; p++)
+ {
+ for (q=p+1; q<node_num; q++)
+ {
+ pm->AddEdge(p, q, Dist(p, q));
+ }
+ }
+ pm->options = options;
+ pm->Solve();
+ for (p=0; p<node_num; p++)
+ {
+ for (q=p+1; q<node_num; q++)
+ {
+ if (pm->GetSolution(e++))
+ {
+ matching[p] = q;
+ matching[q] = p;
+ }
+ }
+ }
+ delete pm;
+ return ComputeCost(matching);
+}
+
+GeomPerfectMatching::REAL GeomPerfectMatching::Solve()
+{
+ double start_time = get_time();
+ double perfect_matching_time = 0;
+ double negative_edges_time = 0;
+ if (options.verbose) { printf("starting geometric matching with %d points\n", node_num); fflush(stdout); }
+ PointId p, q;
+ Edge* e;
+ int _e;
+ int iter;
+ bool success = false;
+ PerfectMatching* pm = NULL;
+ GPMKDTree* kd_tree;
+
+ double init_matching_time = get_time();
+
+ if (gpm_options.init_Delaunay) InitDelaunay();
+ if (gpm_options.init_KNN > 0) InitKNN(gpm_options.init_KNN);
+ if (gpm_options.init_greedy) CompleteInitialMatching();
+
+ init_matching_time = get_time() - init_matching_time;
+
+ graph_update_time = 0;
+
+ int iter_max = gpm_options.iter_max;
+ for (iter=0; iter_max<=0 || iter<iter_max; iter++)
+ {
+ if (pm)
+ {
+ double negative_edges_start_time = get_time();
+ int edge_num0 = edge_num;
+ pm->StartUpdate();
+ for (p=0; p<node_num; p++)
+ {
+ PerfectMatching::REAL s = pm->GetTwiceSum(p);
+ if ( ((REAL)1 / 2) == 0 && ((PerfectMatching::REAL)1 / 2) != 0 ) sums[p] = (REAL)ceil((double)s);
+ else sums[p] = (REAL)s;
+ }
+ if (options.verbose) { printf("building kd_tree..."); fflush(stdout); }
+ {
+ kd_tree = new GPMKDTree(DIM+1, node_num, coords, this);
+ }
+ if (options.verbose) { printf(" done. Now adding negative edges:\n "); fflush(stdout); }
+
+ for (p=0; p<node_num; p++)
+ {
+ if (options.verbose && (p%(node_num/72)==0)) { printf("+"); fflush(stdout); }
+ for (e=nodes[p].first[0]; e; e=e->next[0]) nodes[e->head[0]].is_marked = 1;
+ for (e=nodes[p].first[1]; e; e=e->next[1]) nodes[e->head[1]].is_marked = 1;
+ kd_tree->AddNegativeEdges(p, pm);
+ for (e=nodes[p].first[0]; e; e=e->next[0]) nodes[e->head[0]].is_marked = 0;
+ for (e=nodes[p].first[1]; e; e=e->next[1]) nodes[e->head[1]].is_marked = 0;
+ }
+ delete kd_tree;
+ //if (edge_num - edge_num0 > node_num / 32)
+ if ( 0 ) // always reuse previous computation
+ {
+ delete pm;
+ pm = NULL;
+ }
+ else
+ {
+ pm->FinishUpdate();
+ if (edge_num0 == edge_num) success = true;
+ }
+ if (options.verbose) { printf("\ndone (%d edges added)\n", edge_num-edge_num0); fflush(stdout); }
+ negative_edges_time += get_time() - negative_edges_start_time;
+ }
+ if (!pm)
+ {
+ int E = 5*node_num;
+ if (E < 5*edge_num/4) E = 5*edge_num/4;
+ pm = new PerfectMatching(node_num, E);
+ for (e=edges->ScanFirst(); e; e=edges->ScanNext())
+ {
+ p = e->head[1]; q = e->head[0];
+ pm->AddEdge(p, q, Dist(p, q));
+ }
+ }
+ if (options.verbose) printf("iter %d: ", iter+1);
+ pm->options = options;
+ double perfect_matching_start = get_time();
+ pm->Solve();
+ perfect_matching_time += get_time() - perfect_matching_start;
+ if (success) break;
+ }
+
+ for (_e=0, e=edges->ScanFirst(); e; _e++, e=edges->ScanNext())
+ {
+ if (pm->GetSolution(_e))
+ {
+ p = e->head[1]; q = e->head[0];
+ matching[p] = q;
+ matching[q] = p;
+ }
+ }
+ delete pm;
+ REAL cost = ComputeCost(matching);
+ if (options.verbose)
+ {
+ printf("geometric matching finished [%.3f secs]. cost=%.1f \n", get_time()-start_time, (double)cost);
+ printf(" selecting initial edges: [%.3f secs], perfect matching: [%.3f secs]\n", init_matching_time, perfect_matching_time);
+ printf(" pricing: [%.3f secs] including graph updates: [%.3f secs]\n", negative_edges_time, graph_update_time);
+ fflush(stdout);
+ }
+ return cost;
+}
+
|