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