From a06ff64534815cbf702a3847a19443612d307b80 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 30 Sep 2022 10:55:55 +0200
Subject: Changed rbmp to use blossom algorithm.

---
 blossom5-v2.05.src/GEOM/GPMmain.cpp | 148 ++++++++++++++++++++++++++++++++++++
 1 file changed, 148 insertions(+)
 create mode 100644 blossom5-v2.05.src/GEOM/GPMmain.cpp

(limited to 'blossom5-v2.05.src/GEOM/GPMmain.cpp')

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;
+}
+
-- 
cgit v1.2.3-54-g00ecf