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/GPMinit.cpp           | 200 ++++++++
 blossom5-v2.05.src/GEOM/GPMinterface.cpp      |  84 ++++
 blossom5-v2.05.src/GEOM/GPMkdtree.cpp         | 688 ++++++++++++++++++++++++++
 blossom5-v2.05.src/GEOM/GPMkdtree.h           | 105 ++++
 blossom5-v2.05.src/GEOM/GPMmain.cpp           | 148 ++++++
 blossom5-v2.05.src/GEOM/GeomPerfectMatching.h | 177 +++++++
 6 files changed, 1402 insertions(+)
 create mode 100644 blossom5-v2.05.src/GEOM/GPMinit.cpp
 create mode 100644 blossom5-v2.05.src/GEOM/GPMinterface.cpp
 create mode 100644 blossom5-v2.05.src/GEOM/GPMkdtree.cpp
 create mode 100644 blossom5-v2.05.src/GEOM/GPMkdtree.h
 create mode 100644 blossom5-v2.05.src/GEOM/GPMmain.cpp
 create mode 100644 blossom5-v2.05.src/GEOM/GeomPerfectMatching.h

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

diff --git a/blossom5-v2.05.src/GEOM/GPMinit.cpp b/blossom5-v2.05.src/GEOM/GPMinit.cpp
new file mode 100644
index 0000000..0e31231
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMinit.cpp
@@ -0,0 +1,200 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GeomPerfectMatching.h"
+#include "GPMkdtree.h"
+
+
+
+
+// greedy procedure to make sure that a perfect matching exists:
+// 1. construct a matching among existing edges
+//    (greedy procedure: pick a node, check whether there are edges leading
+//     to unmatched nodes, if there are pick the edge with the smallest length).
+// 2. take remaining unmatched nodes, construct kd-tree for them,
+//    assign an ordering to nodes (last visited time during left-most depth-first search),
+//    add edges between consecutive nodes (2*i,2*i+1)
+void GeomPerfectMatching::CompleteInitialMatching()
+{
+	if (options.verbose) printf("adding edges to make sure that a perfect matching exists...");
+	PointId p, q;
+	Edge* e;
+	double len, len_min;
+	int unmatched_num = 0, edge_num0 = edge_num;
+
+	// construct greedy matching
+	for (p=0; p<node_num; p++)
+	{
+		if (nodes[p].is_marked) continue;
+		q = -1;
+		for (e=nodes[p].first[0]; e; e=e->next[0])
+		{
+			if (nodes[e->head[0]].is_marked) continue;
+			len = Dist2(p, e->head[0]);
+			if (q < 0 || len_min > len)
+			{
+				q = e->head[0];
+				len_min = len;
+			}
+		}
+		if (q >= 0)
+		{
+			nodes[p].is_marked = nodes[q].is_marked = 1;
+		}
+		else unmatched_num ++;
+	}
+
+	if (unmatched_num == 0)
+	{
+		for (p=0; p<node_num; p++) nodes[p].is_marked = 0;
+		return;
+	}
+
+	//printf("%d unmatched\n", unmatched_num);
+
+	REAL* unmatched_coords = new REAL[unmatched_num*DIM];
+	int* rev_mapping = new int[unmatched_num];
+	unmatched_num = 0;
+	for (p=0; p<node_num; p++)
+	{
+		if (nodes[p].is_marked) nodes[p].is_marked = 0;
+		else
+		{
+			memcpy(unmatched_coords+unmatched_num*DIM, coords+p*DIM, DIM*sizeof(REAL));
+			rev_mapping[unmatched_num ++] = p;
+		}
+	}
+
+	GPMKDTree* kd_tree = new GPMKDTree(DIM, unmatched_num, unmatched_coords, this);
+	kd_tree->AddPerfectMatching(rev_mapping);
+
+	delete kd_tree;
+	delete [] unmatched_coords;
+	delete [] rev_mapping;
+
+	if (options.verbose) printf("done (%d edges)\n", edge_num-edge_num0);
+}
+
+void GeomPerfectMatching::InitKNN(int K)
+{
+	if (node_num != node_num_max) { printf("InitKNN() cannot be called before all points have been added!\n"); exit(1); }
+	if (options.verbose) printf("adding K nearest neighbors (K=%d)\n", K);
+
+	int dir, k;
+	PointId p;
+	Edge* e;
+
+	if (K > node_num - 1) K = node_num - 1;
+
+	GPMKDTree* kd_tree = new GPMKDTree(DIM, node_num, coords, this);
+	PointId* neighbors = new PointId[K];
+
+	for (p=0; p<node_num; p++)
+	{
+		for (dir=0; dir<2; dir++)
+		for (e=nodes[p].first[dir]; e; e=e->next[dir])
+		{
+			nodes[e->head[dir]].is_marked = 1;
+		}
+
+		kd_tree->ComputeKNN(p, K, neighbors);
+		for (k=0; k<K; k++)
+		{
+			if (nodes[neighbors[k]].is_marked) continue;
+			AddInitialEdge(p, neighbors[k]);
+			nodes[neighbors[k]].is_marked = 1;
+		}
+
+		for (dir=0; dir<2; dir++)
+		for (e=nodes[p].first[dir]; e; e=e->next[dir])
+		{
+			nodes[e->head[dir]].is_marked = 0;
+		}
+	}
+
+	delete kd_tree;
+	delete [] neighbors;
+}
+
+#ifdef DELAUNAY_TRIANGLE
+
+#ifdef _MSC_VER
+#pragma warning(disable: 4311)
+#pragma warning(disable: 4312)
+#endif
+
+extern "C" {
+#define ANSI_DECLARATORS
+#define TRILIBRARY
+#define NO_TIMER
+#define main NO_MAIN_FUNCTION
+#include "../triangle/triangle.c"
+}
+
+void GeomPerfectMatching::InitDelaunay()
+{
+	if (node_num < 16) return;
+	if (options.verbose) printf("adding edges in Delaunay triangulation\n");
+
+	int k;
+
+	struct triangulateio in, out, vorout;
+	in.numberofpoints = node_num;
+	in.numberofpointattributes = 0;
+	in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));
+	for (k=0; k<2*node_num; k++) in.pointlist[k] = coords[k];
+	in.pointattributelist = NULL;
+	in.pointmarkerlist = NULL;
+	in.numberofsegments = 0;
+	in.numberofholes = 0;
+	in.numberofregions = 0;
+	in.regionlist = 0;
+
+	out.pointlist = (REAL *) NULL;
+	out.pointattributelist = (REAL *) NULL;
+	out.pointmarkerlist = (int *) NULL;
+	out.trianglelist = (int *) NULL;
+	out.triangleattributelist = (REAL *) NULL;
+	out.neighborlist = (int *) NULL;
+	out.segmentlist = (int *) NULL;
+	out.segmentmarkerlist = (int *) NULL;
+	out.edgelist = (int *) NULL;
+	out.edgemarkerlist = (int *) NULL;
+
+	vorout.pointlist = (REAL *) NULL;
+	vorout.pointattributelist = (REAL *) NULL;
+	vorout.edgelist = (int *) NULL;
+	vorout.normlist = (REAL *) NULL;
+
+	triangulate("pczAevn", &in, &out, &vorout);
+
+	free(in.pointlist);
+	free(out.pointlist);
+	free(out.pointmarkerlist);
+	free(out.trianglelist);
+	free(out.neighborlist);
+	free(out.segmentlist);
+	free(out.segmentmarkerlist);
+	free(out.edgemarkerlist);
+	free(vorout.pointlist);
+	free(vorout.pointattributelist);
+	free(vorout.edgelist);
+	free(vorout.normlist);
+
+	for (k=0; k<out.numberofedges; k++) AddInitialEdge(out.edgelist[2*k], out.edgelist[2*k+1]);
+
+	free(out.edgelist);
+}
+
+
+
+#else
+
+void GeomPerfectMatching::InitDelaunay()
+{
+	printf("You need to download the 'Triangle' software from \n\thttp://www.cs.cmu.edu/~quake/triangle.html ,\nextract it to the directory GeomPerfectMatching and define DELAUNARY_TRIANG in GeomPerfectMatching.h\n");
+	exit(1);
+}
+
+#endif
+
diff --git a/blossom5-v2.05.src/GEOM/GPMinterface.cpp b/blossom5-v2.05.src/GEOM/GPMinterface.cpp
new file mode 100644
index 0000000..97a5834
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMinterface.cpp
@@ -0,0 +1,84 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GeomPerfectMatching.h"
+#include "GPMkdtree.h"
+
+
+GeomPerfectMatching::GeomPerfectMatching(int nodeNum, int _DIM)
+	: DIM(_DIM),
+	  node_num(0),
+	  node_num_max(nodeNum),
+	  edge_num(0)
+{
+	if (node_num_max < 1) { printf("too few nodes\n"); exit(1); }
+	if (node_num_max & 1) { printf("# of points is odd: perfect matching cannot exist\n"); exit(1); }
+	nodes = (Node*) malloc(node_num_max*sizeof(Node));
+	memset(nodes, 0, node_num_max*sizeof(Node));
+	edges = new Block<Edge>(512);
+	coords = (REAL*) malloc((DIM+1)*node_num_max*sizeof(REAL));
+	sums = coords + DIM*node_num_max;
+	matching = (int*) malloc(node_num_max*sizeof(int));
+	int i;
+	for (i=0; i<node_num_max; i++) matching[i] = i;
+}
+
+GeomPerfectMatching::~GeomPerfectMatching()
+{
+	free(nodes);
+	delete edges;
+	free(coords);
+	free(matching);
+}
+
+GeomPerfectMatching::PointId GeomPerfectMatching::AddPoint(REAL* coord)
+{
+	if (node_num >= node_num_max)
+	{
+		printf("Error: you are trying to add too many points!\n");
+		exit(1);
+	}
+	memcpy(coords+DIM*node_num, coord, DIM*sizeof(REAL));
+	return node_num ++;
+}
+
+void GeomPerfectMatching::AddInitialEdge(PointId _i, PointId _j)
+{
+	assert(_i>=0 && _i<node_num_max && _j>=0 && _j<node_num_max && _i!=_j);
+	if (_j < _i) { int _k = _i; _i = _j; _j = _k; }
+	Node* i = nodes + _i;
+	Node* j = nodes + _j;
+	Edge* e = edges->New();
+	edge_num ++;
+
+	e->head[1] = _i;
+	e->head[0] = _j;
+	e->next[0] = i->first[0];
+	e->next[1] = j->first[1];
+	i->first[0] = e;
+	j->first[1] = e;
+}
+
+
+GeomPerfectMatching::REAL GeomPerfectMatching::ComputeCost(PointId* matching)
+{
+	if (node_num != node_num_max) { printf("ComputeCost() cannot be called before all points have been added!\n"); exit(1); }
+
+	REAL cost = 0;
+	int i;
+	for (i=0; i<node_num; i++)
+	{
+		if (matching[i]==i || matching[i]<0 || matching[i]>=node_num || matching[matching[i]]!=i)
+		{
+			printf("ComputeCost(): not a valid matching!\n");
+			exit(1);
+		}
+		if (matching[i] > i)
+		{
+			cost += Dist(i, matching[i]);
+		}
+	}
+	return cost;
+}
+
+
diff --git a/blossom5-v2.05.src/GEOM/GPMkdtree.cpp b/blossom5-v2.05.src/GEOM/GPMkdtree.cpp
new file mode 100644
index 0000000..f8cd17b
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMkdtree.cpp
@@ -0,0 +1,688 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GPMkdtree.h"
+#include "../timer.h"
+
+
+// 'mapping' must be of size 2*N. (array[N], ... array[2*N-1] is used as a temporary buffer).
+// After the call  array[mapping[0]] <= array[mapping[1]] <= ... <= array[mapping[N-1]].
+// array is not modified.
+template <typename Type> inline void sort(Type* array, int array_skip, int* mapping, int N)
+{
+	// mergesort
+	int i;
+	int* mappingSrc = mapping;
+	int* mappingDst = mapping + N;
+	int* pSrc1;
+	int* pSrc2;
+	int* pSrc1End;
+	int* pSrc2End;
+	int* pDst;
+
+	for (i=0; i<(N&(~1)); i+=2)
+	{
+		if (array[array_skip*i] < array[array_skip*(i+1)])
+		{
+			mappingSrc[i]   = i;
+			mappingSrc[i+1] = i+1;
+		}
+		else
+		{
+			mappingSrc[i]   = i+1;
+			mappingSrc[i+1] = i;
+		}
+	}
+	if (i != N) mappingSrc[i] = i;
+
+	int step;
+	for (step=2; step<N; step*=2)
+	{
+		pSrc2End = mappingSrc;
+		pDst = mappingDst;
+		while ( 1 )
+		{
+			pSrc1 = pSrc2End;
+			pSrc1End = pSrc1 + step;
+			if (pSrc1End >= mappingSrc + N)
+			{
+				memcpy(pDst, pSrc1, (int)((char*)(mappingSrc + N) - (char*)pSrc1));
+				break;
+			}
+			pSrc2 = pSrc1End;
+			pSrc2End = pSrc2 + step;
+			if (pSrc2End > mappingSrc + N) pSrc2End = mappingSrc + N;
+			while ( 1 )
+			{
+				if (array[(array_skip)*(*pSrc1)] < array[array_skip*(*pSrc2)])
+				{
+					*pDst ++ = *pSrc1 ++;
+					if (pSrc1 == pSrc1End)
+					{
+						memcpy(pDst, pSrc2, (int)((char*)pSrc2End - (char*)pSrc2));
+						pDst = (int*) ((char*)pDst + (int)((char*)pSrc2End - (char*)pSrc2));
+						break;
+					}
+				}
+				else
+				{
+					*pDst ++ = *pSrc2 ++;
+					if (pSrc2 == pSrc2End)
+					{
+						memcpy(pDst, pSrc1, (int)((char*)pSrc1End - (char*)pSrc1));
+						pDst = (int*) ((char*)pDst + (int)((char*)pSrc1End - (char*)pSrc1));
+						break;
+					}
+				}
+			}
+		}
+		pDst = mappingDst;
+		mappingDst = mappingSrc;
+		mappingSrc = pDst;
+	}
+	if (mappingSrc != mapping) memcpy(mapping, mappingSrc, N*sizeof(int));
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+#define NEIGHBOR_PARENT(k) (((k)-1)>>1)
+#define NEIGHBOR_FIRST_CHILD(k) (((k)<<1)+1)
+
+Neighbors::Neighbors()
+{
+	K_max = 0;
+	dist_array = NULL;
+}
+
+Neighbors::~Neighbors()
+{
+	if (dist_array) delete [] dist_array;
+}
+
+void Neighbors::Init(int _K, PointId* _array)
+{
+	K = _K;
+	array = _array;
+	num = 0;
+	if (K > K_max)
+	{
+		if (dist_array) delete [] dist_array;
+		K_max = K;
+		dist_array = new double[K_max];
+	}
+}
+
+inline void Neighbors::Swap(int k1, int k2)
+{
+	PointId p = array[k1]; array[k1] = array[k2]; array[k2] = p;
+	double d = dist_array[k1]; dist_array[k1] = dist_array[k2]; dist_array[k2] = d;
+}
+
+inline void Neighbors::Add(PointId p, double dist)
+{
+	int k;
+	if (num < K)
+	{
+		k = num ++;
+		array[k] = p;
+		dist_array[k] = dist;
+		while ( k > 0 )
+		{
+			int k_parent = NEIGHBOR_PARENT(k);
+			if (dist_array[k] <= dist_array[k_parent]) break;
+			Swap(k, k_parent);
+			k = k_parent;
+		}
+	}
+	else
+	{
+		if (dist_array[0] <= dist) return;
+		array[0] = p;
+		dist_array[0] = dist;
+		k = 0;
+		while ( 1 )
+		{
+			int k_child = NEIGHBOR_FIRST_CHILD(k);
+			if (k_child >= K) break;
+			if (k_child+1 < K && dist_array[k_child+1] > dist_array[k_child]) k_child ++;
+			if (dist_array[k] >= dist_array[k_child]) break;
+			Swap(k, k_child);
+			k = k_child;
+		}
+	}
+	//for (k=1; k<num; k++) assert(dist_array[k] <= dist_array[NEIGHBOR_PARENT(k)]);
+}
+
+inline double Neighbors::GetMax()
+{
+	assert(num > 0);
+	return dist_array[0];
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+GPMKDTree::GPMKDTree(int _D, int _point_num, REAL* coords, GeomPerfectMatching* _GPM)
+	: D(_D), DIM(_GPM->DIM), point_num(_point_num), GPM(_GPM)
+{
+	Node* i;
+	Node* j;
+	int d, d0, k;
+	int* mapping = new int[(D+2)*point_num];
+	int* buf = mapping + D*point_num;
+	int* marking = buf + point_num;
+	memset(marking, 0, point_num*sizeof(int));
+	int* ptr = mapping;
+
+	int node_num_max = 4*point_num/3+2;
+	nodes = (Node*)malloc(node_num_max*sizeof(Node));
+	rev_mapping = (Node**)malloc(point_num*sizeof(Node*));
+	memset(rev_mapping, 0, point_num*sizeof(Node*));
+
+	REAL** coords_array = new REAL*[D];
+	int* skip_array = new int[D];
+	for (d=0; d<DIM; d++) { coords_array[d] = coords + d; skip_array[d] = DIM; }
+	if (d < D)            { coords_array[d] = GPM->sums; skip_array[d] = 1; }
+
+	for (d=0; d<D; d++)
+	{
+		sort<REAL>(coords_array[d], skip_array[d], ptr, point_num);
+		if (d == DIM) sum_max = GPM->sums[ptr[point_num-1]];
+		ptr += point_num;
+	}
+
+	nodes[0].parent = NULL;
+	nodes[0].order = 0;
+	nodes[0].d = point_num;
+	nodes[0].first_child = NULL;
+	node_num = 1;
+	Node* first_unprocessed = &nodes[0];
+	while ( (i=first_unprocessed) )
+	{
+		first_unprocessed = i->first_child;
+
+		int start = i->order;
+		int num0 = i->d, num;
+		if ((DIM==D && num0<=2) || (DIM<D && num0 <= CHILDREN_MAX))
+		{
+			// leaf
+			i->d = -num0;
+			for (k=0; k<num0; k++)
+			{
+				i->points[k] = mapping[start+k];
+				rev_mapping[mapping[start+k]] = i;
+			}
+			continue;
+		}
+
+		// not a leaf.
+		if (node_num + 2 > node_num_max)
+		{
+			node_num_max = 3*node_num_max/2 + 16;
+			Node* nodes_old = nodes;
+			nodes = (Node*)realloc(nodes, node_num_max*sizeof(Node));
+#define UPDATE_NODE_PTR(ptr) ptr = (Node*)((char*)ptr + ((char*)nodes-(char*)nodes_old))
+			UPDATE_NODE_PTR(i);
+			if (first_unprocessed) UPDATE_NODE_PTR(first_unprocessed);
+			for (k=0; k<node_num; k++)
+			{
+				if (nodes[k].parent)                       UPDATE_NODE_PTR(nodes[k].parent);
+				if (nodes[k].d>=0 && nodes[k].first_child) UPDATE_NODE_PTR(nodes[k].first_child);
+			}
+			for (k=0; k<point_num; k++) if (rev_mapping[k]) UPDATE_NODE_PTR(rev_mapping[k]);
+		}
+
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		// choose split dimension d0 and number of items 'num' in the first group.
+
+		//d0 = (i->parent) ? ((i->parent->d + 1) % D) : 0;
+		//num = num0/2;
+
+		const int FRACTION = 20;
+		int num_min = 1;      if (num_min < num0/FRACTION) num_min = num0/FRACTION;
+		int num_max = num0-1; if (num_max > (FRACTION-1)*num0/FRACTION) num_max = (FRACTION-1)*num0/FRACTION;
+		int num_max_DIM = num0-1;
+
+		if (D>DIM && (!i->parent || !i->parent->parent || !i->parent->parent->parent)) d0 = D-1;
+		else
+		{
+			d0 = -1;
+			REAL diff_max = 0;
+			for (d=0; d<D; d++)
+			{
+				ptr = mapping + d*point_num;
+
+				REAL c_min, c_max;
+				c_min = coords_array[d][ptr[start +          num_min             ]*skip_array[d]];
+				c_max = coords_array[d][ptr[start + ((d<DIM)?num_max:num_max_DIM)]*skip_array[d]];
+				REAL diff = c_max - c_min;
+				//if (d == DIM) diff *= 2;
+				if (d0<0 || diff_max < diff) { d0 = d; diff_max = diff; }
+			}
+		}
+
+		ptr = mapping + d0*point_num;
+		REAL c_min, c_max;
+		if (d < DIM)
+		{
+			c_min = coords_array[d0][ptr[start+num_min]*skip_array[d0]];
+			c_max = coords_array[d0][ptr[start+num_max]*skip_array[d0]];
+		}
+		else
+		{
+			c_min = coords_array[d0][ptr[start+1]*skip_array[d0]];
+			c_max = coords_array[d0][ptr[start+num0-1]*skip_array[d0]]; // for this dimension upper tail is important, so don't discard it!
+		}
+		REAL split = (c_min + c_max) / 2;
+		for (num=num_min; num<num_max; num++)
+		{
+			if (coords_array[d0][ptr[start+num]*skip_array[d0]] > split) break;
+		}
+
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+
+		ptr = mapping + d0*point_num;
+		for (k=start; k<start+num; k++) marking[ptr[k]] = 1;
+		for (d=0; d<D; d++)
+		{
+			if (d == d0) continue;
+			ptr = mapping + d*point_num;
+			int* bufa = buf;
+			int* bufb = buf + num;
+			for (k=start; k<start+num0; k++)
+			{
+				if (marking[ptr[k]]) *bufa ++ = ptr[k];
+				else                 *bufb ++ = ptr[k];
+			}
+			memcpy(ptr+start, buf, num0*sizeof(int));
+		}
+		ptr = mapping + d0*point_num;
+		for (k=start; k<start+num; k++) marking[ptr[k]] = 0;
+
+		i->d = d0;
+		PointId p = ptr[start + num - ((num >= num0-num) ? 1 : 0)];
+		i->coord = coords_array[d0][p*skip_array[d0]];
+		i->first_child = j = &nodes[node_num];
+		node_num += 2;
+		j->parent = (j+1)->parent = i;
+		j->order = start;
+		j->d = num;
+		(j+1)->order = start+num;
+		(j+1)->d = num0-num;
+		(j+1)->first_child = first_unprocessed;
+		j->first_child = j+1;
+		first_unprocessed = j;
+	}
+
+	delete [] coords_array;
+	delete [] skip_array;
+	delete [] mapping;
+
+	////////////////////////////////////////////////////////////////////////////////////
+	// set ordering and depth_max
+	int depth = 0, depth_max = 0;
+
+	// set ordering
+	i = &nodes[0];
+	k = 0;
+	if (D > DIM)
+	{
+		// ordering for AddNegativeEdges() - tree preorder
+		while ( 1 )
+		{
+			if (!IS_LEAF(i))
+			{
+				i = i->first_child;
+				depth ++;
+				if (depth_max < depth) depth_max = depth;
+			}
+			else
+			{
+				while ( 1 )
+				{
+					i->order = k ++;
+
+					if (!i->parent) break;
+					if (i->parent->first_child == i) { i ++; break; }
+					i = i->parent;
+					depth --;
+				}
+				if (!i->parent) break;
+			}
+		}
+	}
+	else
+	{
+		// compute tree inorder - useful for nearest neighbor search so that that branches close to the input node are explored first
+		while ( 1 )
+		{
+			if (!IS_LEAF(i))
+			{
+				i = i->first_child;
+				depth ++;
+				if (depth_max < depth) depth_max = depth;
+			}
+			else
+			{
+				i->order = k ++;
+				while ( 1 )
+				{
+					if (!i->parent) break;
+					if (i->parent->first_child == i) { i->parent->order = k ++; i ++; break; }
+					i = i->parent;
+					depth --;
+				}
+				if (!i->parent) break;
+			}
+		}
+	}
+
+	traversing_buf = (REAL*) malloc((D + depth_max + 2)*sizeof(REAL));
+}
+
+GPMKDTree::~GPMKDTree()
+{
+	free(nodes);
+	free(rev_mapping);
+	free(traversing_buf);
+}
+
+void GPMKDTree::AddPerfectMatching(PointId* rev_mapping)
+{
+	Node* i;
+	int k;
+	PointId p, q = -1;
+	i = &nodes[0];
+	do
+	{
+		if (IS_LEAF(i))
+		{
+			for (k=0; k<-i->d; k++)
+			{
+				p = i->points[k];
+				if (q < 0) q = p;
+				else { GPM->AddInitialEdge(rev_mapping[p], rev_mapping[q]); q = -1; }
+			}
+		}
+		else
+		{
+			i = i->first_child;
+			continue;
+		}
+
+		while ( i->parent )
+		{
+			if (i->parent->first_child == i) { i ++; break; }
+			i = i->parent;
+		}
+	} while (i->parent);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+#define MOVE_DOWN_LEFT(i)\
+	{\
+		*stack ++ = current_diff[i->d];\
+		if (current_diff[i->d] <= 0 && (diff = i->coord - coord0[i->d]) < 0) current_diff[i->d] = diff;\
+		i = i->first_child;\
+	}
+
+#define MOVE_DOWN_RIGHT(i)\
+	{\
+		*stack ++ = current_diff[i->d];\
+		if (current_diff[i->d] >= 0 && (diff = i->coord - coord0[i->d]) > 0) current_diff[i->d] = diff;\
+		i = i->first_child+1;\
+	}
+
+#define MOVE_LEFT(i)\
+	{\
+		int d_prev = i->parent->d;\
+		current_diff[d_prev] = stack[-1];\
+		if (current_diff[d_prev] <= 0 && (diff = i->parent->coord - coord0[d_prev]) < 0) current_diff[d_prev] = diff;\
+		i --;\
+	}
+
+#define MOVE_RIGHT(i)\
+	{\
+		int d_prev = i->parent->d;\
+		current_diff[d_prev] = stack[-1];\
+		if (current_diff[d_prev] >= 0 && (diff = i->parent->coord - coord0[d_prev]) > 0) current_diff[d_prev] = diff;\
+		i ++;\
+	}
+
+#define MOVE_UP(i)\
+	{\
+		i = i->parent;\
+		current_diff[i->d] = *(-- stack);\
+	}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+#define MOVE_DOWN_LEFT_X(i)\
+	{\
+		*stack ++ = current_diff[i->d];\
+		if (i->d < D-1) { if (current_diff[i->d] <= 0 && (diff = i->coord - coord0[i->d]) < 0) current_diff[i->d] = diff; }\
+		else current_diff[i->d] = i->coord;\
+		i = i->first_child;\
+	}
+
+#define MOVE_DOWN_RIGHT_X(i)\
+	{\
+		*stack ++ = current_diff[i->d];\
+		if (i->d < D-1) { if (current_diff[i->d] >= 0 && (diff = i->coord - coord0[i->d]) > 0) current_diff[i->d] = diff; }\
+		i = i->first_child+1;\
+	}
+
+#define MOVE_LEFT_X(i)\
+	{\
+		int d_prev = i->parent->d;\
+		current_diff[d_prev] = stack[-1];\
+		if (d_prev < D-1) { if (current_diff[d_prev] <= 0 && (diff = i->parent->coord - coord0[d_prev]) < 0) current_diff[d_prev] = diff; }\
+		else current_diff[d_prev] = i->parent->coord;\
+		i --;\
+	}
+
+#define MOVE_RIGHT_X(i)\
+	{\
+		int d_prev = i->parent->d;\
+		current_diff[d_prev] = stack[-1];\
+		if (d_prev < D-1) { if (current_diff[d_prev] >= 0 && (diff = i->parent->coord - coord0[d_prev]) > 0) current_diff[d_prev] = diff; }\
+		i ++;\
+	}
+
+#define MOVE_UP_X(i)\
+	{\
+		i = i->parent;\
+		current_diff[i->d] = *(-- stack);\
+	}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+/*
+void GPMKDTree::ComputeKNN(PointId p0, int K, PointId* neighbors_array)
+{
+	int p;
+	REAL* coord0 = GPM->coords + p0*DIM;
+
+	neighbors.Init(K, neighbors_array);
+
+	for (p=0; p<point_num; p++)
+	if (p != p0)
+	{
+		neighbors.Add(p, GPM->Dist(coord0, GPM->coords+p*DIM));
+	}
+}
+*/
+
+void GPMKDTree::ComputeKNN(PointId p0, int K, PointId* neighbors_array)
+{
+	int neighbor_num = 0;
+	Node* i = rev_mapping[p0];
+	int k, order0 = i->order;
+	REAL diff;
+	REAL* coords = GPM->coords;
+	REAL* coord0 = GPM->coords + p0*DIM;
+	REAL* current_diff = traversing_buf;
+	REAL* stack = traversing_buf + D;
+
+	neighbors.Init(K, neighbors_array);
+
+	for (k=0; k<D; k++) current_diff[k] = 0;
+	i = &nodes[0];
+	do
+	{
+		if (IS_LEAF(i))
+		{
+			for (k=0; k<-i->d; k++)
+			{
+				PointId p = i->points[k];
+				if (p == p0) continue;
+				double dist2;
+				REAL* coord2 = coords+p*DIM;
+				GPM_GET_DIST2(dist2, coord0, coord2);
+				neighbors.Add(p, dist2);
+			}
+		}
+		else
+		{
+			if (neighbors.GetNum() < K || GPM->Norm2(current_diff) < neighbors.GetMax())
+			{
+				if (i->order > order0)
+				{
+					MOVE_DOWN_LEFT(i);
+				}
+				else
+				{
+					MOVE_DOWN_RIGHT(i);
+				}
+				continue;
+			}
+		}
+
+		while ( i->parent )
+		{
+			if (i->parent->order > order0)
+			{
+				if (i->parent->first_child == i)
+				{
+					MOVE_RIGHT(i);
+					break;
+				}
+			}
+			else
+			{
+				if (i->parent->first_child != i)
+				{
+					MOVE_LEFT(i);
+					break;
+				}
+			}
+			MOVE_UP(i);
+		}
+	} while ( i->parent );
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+/*
+void GPMKDTree::AddNegativeEdges(PointId p, PerfectMatching* pm)
+{
+	PointId q;
+	for (q=p+1; q<node_num; q++)
+	{
+		if (GPM->nodes[q].is_marked) continue;
+		REAL len = GPM->Dist(p, q);
+		if (2*len - GPM->sums[p] - GPM->sums[q] < 0)
+		{
+			if (pm->AddNewEdge(p, q, len, true)>=0)
+			{
+				GPM->AddInitialEdge(p, q);
+			}
+		}
+	}
+}
+*/
+
+
+void GPMKDTree::AddNegativeEdges(PointId p0, PerfectMatching* pm)
+{
+	Node* i = rev_mapping[p0];
+	int k, order0 = i->order;
+	bool check;
+	REAL diff;
+	REAL* coords = GPM->coords;
+	REAL* coord0 = coords + p0*DIM;
+	REAL* sums = GPM->sums;
+	REAL sum0 = sums[p0];
+	REAL* current_diff = traversing_buf;
+	REAL* stack = traversing_buf + D;
+
+	for (k=0; i->points[k]!=p0; k++) {}
+	for (k++; k<-i->d; k++)
+	{
+		PointId p = i->points[k];
+		
+		REAL len = GPM->Dist(coord0, GPM->coords+p*DIM);
+		if (2*len - GPM->sums[p] < GPM->sums[p0])
+		{
+			double start_time = get_time();
+			if (pm->AddNewEdge(p0, p, len, true)>=0) GPM->AddInitialEdge(p0, p);
+			GPM->graph_update_time += get_time() - start_time;
+		}
+	}
+
+	for (k=0; k<DIM; k++) current_diff[k] = 0;
+	current_diff[k] = sum_max;
+	i = &nodes[0];
+	do
+	{
+		if (i->order > order0)
+		{
+			if (IS_LEAF(i))
+			{
+				for (k=0; k<-i->d; k++)
+				{
+					PointId p = i->points[k];
+					if (!GPM->nodes[p].is_marked)
+					{
+						//REAL* coord2 = coords+p*DIM;
+						//REAL threshold = sums[p0]+sums[p];
+						//GPM_CHECK_DIST(check, coord0, coord2, threshold);
+						//if (check)
+						{
+							REAL len = GPM->Dist(coord0, GPM->coords+p*DIM);
+							if (2*len - GPM->sums[p] < GPM->sums[p0])
+							{
+								double start_time = get_time();
+								if (pm->AddNewEdge(p0, p, len, true)>=0) GPM->AddInitialEdge(p0, p);
+								GPM->graph_update_time += get_time() - start_time;
+							}
+						}
+					}
+				}
+			}
+			else
+			{
+				REAL threshold = current_diff[D-1] + sum0;
+				GPM_CHECK_NORM(check, current_diff, threshold);
+				if (check)
+				{
+					MOVE_DOWN_LEFT_X(i);
+					continue;
+				}
+			}
+		}
+
+		while ( i->parent )
+		{
+			if (i->parent->first_child == i)
+			{
+				MOVE_RIGHT_X(i);
+				break;
+			}
+			MOVE_UP_X(i);
+		}
+	} while (i->parent);
+}
diff --git a/blossom5-v2.05.src/GEOM/GPMkdtree.h b/blossom5-v2.05.src/GEOM/GPMkdtree.h
new file mode 100644
index 0000000..67e3d91
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMkdtree.h
@@ -0,0 +1,105 @@
+/*
+    GPMkdtree.h - kd-tree data structure for pricing in complete geometric instances
+
+    Copyright 2008 Vladimir Kolmogorov (vnk@adastral.ucl.ac.uk)
+
+    This software can be used for research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
+#ifndef NJKASKJTASJNBAJSNRBAJS
+#define NJKASKJTASJNBAJSNRBAJS
+
+#include "GeomPerfectMatching.h"
+
+struct Neighbors
+{
+	typedef GeomPerfectMatching::REAL REAL;
+	typedef GeomPerfectMatching::PointId PointId;
+
+	Neighbors();
+	~Neighbors();
+
+	int GetNum() { return num; }
+
+	void Init(int K, PointId* array);
+	void Add(PointId p, double dist);
+	double GetMax();
+
+private:
+	void Swap(int k1, int k2);
+	PointId* array;
+	double* dist_array;
+	int num, K, K_max;
+};
+
+
+
+struct GPMKDTree
+{
+	typedef GeomPerfectMatching::REAL REAL;
+	typedef GeomPerfectMatching::PointId PointId;
+
+	GPMKDTree(int D, int point_num, REAL* coords, GeomPerfectMatching* GPM);
+	~GPMKDTree();
+
+	// if D == DIM
+	void AddPerfectMatching(PointId* rev_mapping);
+	void ComputeKNN(PointId p, int K, PointId* neighbors);
+
+	// if D == DIM+1
+	void AddNegativeEdges(PointId p, PerfectMatching* pm);
+
+	//////////////////////////////////////////////////////////////////////////
+
+public:
+#define CHILDREN_MAX 2
+	struct Node
+	{
+		Node* parent;
+		int d; // split dimension. d<0 indicates a leaf.
+#define IS_LEAF(i) ((i)->d < 0)
+		union
+		{
+			struct // for non-leaves
+			{
+				REAL coord;
+				Node* first_child; // the second child is first_child+1
+			};
+			struct // for leaves
+			{
+				PointId points[CHILDREN_MAX]; // the number of points is -d
+			};
+		};
+		int order;
+	}* nodes;
+
+
+	//////////////////////////////////////////////////////////////////////////
+	Node** rev_mapping;
+
+	int D, DIM, point_num, node_num;
+	REAL sum_max;
+	REAL* traversing_buf;
+	GeomPerfectMatching* GPM;
+
+	Neighbors neighbors;
+};
+
+
+
+
+#endif
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;
+}
+
diff --git a/blossom5-v2.05.src/GEOM/GeomPerfectMatching.h b/blossom5-v2.05.src/GEOM/GeomPerfectMatching.h
new file mode 100644
index 0000000..b54d471
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GeomPerfectMatching.h
@@ -0,0 +1,177 @@
+/*
+    GeomPerfectMatching.h - computing min cost perfect matching in complete geometric instances (with edge weights equal to Euclidean distances)
+
+    Copyright 2008 Vladimir Kolmogorov (vnk@adastral.ucl.ac.uk)
+
+    This software can be used for research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef ASKHALKSJBRASMNABFAJSTAS
+#define ASKHALKSJBRASMNABFAJSTAS
+
+#include <assert.h>
+#include <math.h>
+#include "../PerfectMatching.h"
+
+//#define DELAUNAY_TRIANGLE
+
+struct GPMKDTree;
+
+class GeomPerfectMatching
+{
+public:
+	typedef int REAL; // if you change it to double, you should also change PerfectMatching::REAL to double!
+	typedef int PointId;
+
+	// pointNum must be a positive even number.
+	GeomPerfectMatching(int pointNum, int DIM);
+	~GeomPerfectMatching();
+
+	// Must be called exactly point_num times (before anything else can be called).
+	// First call returns 0, second returns 1, and so on.
+	// coord must be an array of size DIM. (This array is read into internal memory.)
+	PointId AddPoint(REAL* coord);
+
+
+	////////////// SOLVING //////////////
+
+	// solves perfect matching in the complete graph. Inefficient, just for testing.
+	REAL SolveComplete(); 
+
+
+	// options for Solve()
+	struct GPMOptions
+	{
+		GPMOptions() : init_Delaunay(true), init_KNN(0), init_greedy(true), iter_max(0) {}
+
+		// three variables below determine the initial subset of edges. 
+		// Delaunay initialization seems to be more robust than K nearest neighbors
+		// (but you need to download the "Triangle" package of Shewchuk
+		// from http://www.cs.cmu.edu/~quake/triangle.html , extract it to the directory 'triangle'
+		// and define DELAUNAY_TRIANGLE above). 
+		bool	init_Delaunay;  // add Delaunay triangulation edges.
+		int		init_KNN;       // use init_KNN nearest neighbors for each point.
+		bool	init_greedy;    // add edges greedily to make sure that a perfect matching exists (see comments before CompleteInitialMatching() in GPMinit.cpp for details)
+
+		int		iter_max;   // If iter_max <= 0 then adds subsets of edges until an optimal solution is found.
+				            // Otherwise runs at most iter_max iterations, so the solution may be suboptimal. 
+				            // (iter_max=1 runs perfect matching just for the initial subset).
+	};
+	struct PerfectMatching::Options options;
+	struct GPMOptions gpm_options;
+
+	REAL Solve(); 
+
+	// You can also specify the initial subset manually
+	void AddInitialEdge(PointId i, PointId j);
+
+	//////////// READING RESULTS //////////
+	// Can be called after Solve() or SolveComplete(). Returns pointer to an array 'matching' 
+	// of size node_num. (matching[i] is the point corresponding to point i).
+	// User must not modify this array.
+	PointId GetMatch(PointId p) { return matching[p]; }
+
+
+
+	REAL Dist(REAL* coord1, REAL* coord2);
+	REAL Dist(PointId p, PointId q);
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+private:
+
+friend struct GPMKDTree;
+
+	struct Edge
+	{
+		PointId		head[2];
+		Edge*		next[2];
+	};
+	struct Node
+	{
+		Edge*		first[2];
+		int			is_marked;
+	};
+
+	Node*			nodes;
+	Block<Edge>*	edges;
+	REAL*			coords; // array of size DIM*node_num_max
+	REAL*			sums; // array of size node_num_max
+	PointId*		matching; // array of size node_num_max
+
+	int DIM;
+	int node_num, node_num_max;
+	int edge_num;
+
+	double graph_update_time;
+
+// below x denotes REAL, X denotes double, c,c1,c2 are coords (REAL*)
+#define GPM_ROUND(X) (REAL)( ( ((REAL)1 / 2) == 0 ) ? ((X)+0.5) : (X) )
+#define GPM_GET_NORM2(X, c) { int d; X = 0; for (d=0; d<DIM; d++) X += ((double)(c)[d])*(c)[d]; }
+#define GPM_GET_NORM(x, c) { double X; GPM_GET_NORM2(X, c); X = sqrt(X); x = GPM_ROUND(X); }
+#define GPM_GET_DIST2(X, c1, c2) { int d; X = 0; for (d=0; d<DIM; d++) X += ((double)((c1)[d]-(c2)[d]))*((c1)[d]-(c2)[d]); }
+#define GPM_GET_DIST(x, c1, c2) { double X; GPM_GET_DIST2(X, c1, c2); X = sqrt(X); x = GPM_ROUND(X); }
+// if returns false, then 2*||coord||-treshold >= 0
+#define GPM_CHECK_NORM(result, c, x)\
+	{\
+		if (threshold <= 0) result = false;\
+		else\
+		{\
+			double X;\
+			GPM_GET_NORM2(X, c);\
+			result = (4*X < (double)x*threshold);\
+		}\
+	}
+#define GPM_CHECK_DIST(result, c1, c2, x)\
+	{\
+		if (threshold <= 0) result = false;\
+		else\
+		{\
+			double X;\
+			GPM_GET_DIST2(X, c1, c2);\
+			result = (4*X < (double)x*threshold);\
+		}\
+	}
+
+	double Norm2(REAL* coord)                { double norm2; GPM_GET_NORM2(norm2, coord); return norm2; }
+	REAL Norm(REAL* coord)                   { REAL norm;    GPM_GET_NORM (norm,  coord); return norm;  }
+	double Dist2(REAL* coord1, REAL* coord2) { double dist2; GPM_GET_DIST2(dist2, coord1, coord2); return dist2; }
+	double Dist2(PointId p, PointId q) { return Dist2(coords+DIM*p, coords+DIM*q); }
+
+	void CompleteInitialMatching(); // add edges so that a perfect matching is guaranteed to exist
+	void InitKNN(int K);
+	void InitDelaunay();
+
+	REAL ComputeCost(PointId* matching);
+};
+
+inline GeomPerfectMatching::REAL GeomPerfectMatching::Dist(REAL* coord1, REAL* coord2)
+{
+	REAL dist;
+	GPM_GET_DIST (dist,  coord1, coord2);
+	return dist; 
+}
+
+inline GeomPerfectMatching::REAL GeomPerfectMatching::Dist(PointId p, PointId q)
+{
+	return Dist(coords+DIM*p,  coords+DIM*q);
+}
+
+#endif
-- 
cgit v1.2.3-70-g09d2