Unverified Commit 776fe75b authored by dwuggh's avatar dwuggh
Browse files

add blossomV source code

parent f0c9f433
Changes from version 1.0 to 2.0:
- Replaced Fibonacci heaps with pairing heaps
(M. Fredman, R. Sedgewick, D. Sleator, R. Tarjan, Algorithmica 1(1):111-129, 1986).
Pairing heaps take less memory (namely, 2 pointers per edge less) compared to Fibonacci heaps,
and seem to be marginally faster.
Finonacci heaps are still available - replace PQ.h with PQ-Fibonacci.h .
- Changed data structures so that the time in SHRINK operations
is O(m log n) per augmentation (I believe). This was not the case
in version 1.0 (see footnote 5 in the MPC paper).
I re-ran experiments corresponding to tables 1,3,4,5,9 in the paper.
The new version was marginally faster (e.g. up to 10% faster) on all examples
except for lrb744710, where it was about 3 times faster.
Changes from version 2.0 to 2.01 (thanks to Nic Schraudolph and Dmitry Kamenetsky for useful suggestions):
- Fixed bug in block.h (replaced "new char[] ... delete ..." with "new char[] ... delete[] ...";
in the first case the behavious is not specified, though most compilers would compile it correctly.)
- Removed PQ-Fibonacci.h
- Added disclaimer about using floating point numbers
Changes from version 2.01 to 2.02:
- Tweaks to stop compiler warnings,
changed "delete rev_mapping" to "delete [] rev_mapping" in misc.cpp (thanks to Nic Schraudolph for suggestions)
- Added a statement to license conditions
Changes from version 2.02 to 2.03:
- Fixed a bug: if the number of expands was too large (more than node_num/4) then the previous version could have crashed.
This was more likely to happen when dynamic updates were used (with multiple calls to PerfectMatching::Solve()).
Changes from version 2.03 to 2.04:
- Changes to make it compile under OS X and with latest gcc compilers, such as gcc 4.7.0 or 4.8.0
(thanks to Adam Whiteside and Michael Cook for suggesting these changes).
#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
#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;
}
#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;