From a06ff64534815cbf702a3847a19443612d307b80 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 30 Sep 2022 10:55:55 +0200 Subject: Changed rbmp to use blossom algorithm. --- .gitignore | 2 + Makefile | 24 + blossom5-v2.05.src/CHANGES.TXT | 45 ++ 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 +++++++ blossom5-v2.05.src/GRAPH1.TXT | 10 + blossom5-v2.05.src/GRAPH2.TXT | 11 + blossom5-v2.05.src/LCA.h | 297 +++++++++++ blossom5-v2.05.src/LICENSE.TXT | 25 + blossom5-v2.05.src/Makefile | 26 + blossom5-v2.05.src/MinCost/MinCost.cpp | 287 +++++++++++ blossom5-v2.05.src/MinCost/MinCost.cpp.orig | 287 +++++++++++ blossom5-v2.05.src/MinCost/MinCost.cpp.rej | 25 + blossom5-v2.05.src/MinCost/MinCost.h | 504 +++++++++++++++++++ blossom5-v2.05.src/MinCost/instances.inc | 12 + blossom5-v2.05.src/MinCost/p | 26 + blossom5-v2.05.src/PMduals.cpp | 441 ++++++++++++++++ blossom5-v2.05.src/PMexpand.cpp | 304 +++++++++++ blossom5-v2.05.src/PMimplementation.h | 356 +++++++++++++ blossom5-v2.05.src/PMinit.cpp | 544 ++++++++++++++++++++ blossom5-v2.05.src/PMinterface.cpp | 270 ++++++++++ blossom5-v2.05.src/PMmain.cpp | 694 ++++++++++++++++++++++++++ blossom5-v2.05.src/PMrepair.cpp | 431 ++++++++++++++++ blossom5-v2.05.src/PMshrink.cpp | 318 ++++++++++++ blossom5-v2.05.src/PQ.h | 396 +++++++++++++++ blossom5-v2.05.src/PerfectMatching.h | 246 +++++++++ blossom5-v2.05.src/README.TXT | 80 +++ blossom5-v2.05.src/USAGE.TXT | 31 ++ blossom5-v2.05.src/block.h | 268 ++++++++++ blossom5-v2.05.src/example.cpp | 285 +++++++++++ blossom5-v2.05.src/misc.cpp | 172 +++++++ blossom5-v2.05.src/timer.h | 98 ++++ hungarian.cpp | 14 +- rbmp.cpp | 106 ++-- 38 files changed, 7972 insertions(+), 65 deletions(-) create mode 100644 .gitignore create mode 100644 Makefile create mode 100644 blossom5-v2.05.src/CHANGES.TXT 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 create mode 100644 blossom5-v2.05.src/GRAPH1.TXT create mode 100644 blossom5-v2.05.src/GRAPH2.TXT create mode 100644 blossom5-v2.05.src/LCA.h create mode 100644 blossom5-v2.05.src/LICENSE.TXT create mode 100644 blossom5-v2.05.src/Makefile create mode 100644 blossom5-v2.05.src/MinCost/MinCost.cpp create mode 100644 blossom5-v2.05.src/MinCost/MinCost.cpp.orig create mode 100644 blossom5-v2.05.src/MinCost/MinCost.cpp.rej create mode 100644 blossom5-v2.05.src/MinCost/MinCost.h create mode 100644 blossom5-v2.05.src/MinCost/instances.inc create mode 100644 blossom5-v2.05.src/MinCost/p create mode 100644 blossom5-v2.05.src/PMduals.cpp create mode 100644 blossom5-v2.05.src/PMexpand.cpp create mode 100644 blossom5-v2.05.src/PMimplementation.h create mode 100644 blossom5-v2.05.src/PMinit.cpp create mode 100644 blossom5-v2.05.src/PMinterface.cpp create mode 100644 blossom5-v2.05.src/PMmain.cpp create mode 100644 blossom5-v2.05.src/PMrepair.cpp create mode 100644 blossom5-v2.05.src/PMshrink.cpp create mode 100644 blossom5-v2.05.src/PQ.h create mode 100644 blossom5-v2.05.src/PerfectMatching.h create mode 100644 blossom5-v2.05.src/README.TXT create mode 100644 blossom5-v2.05.src/USAGE.TXT create mode 100644 blossom5-v2.05.src/block.h create mode 100644 blossom5-v2.05.src/example.cpp create mode 100644 blossom5-v2.05.src/misc.cpp create mode 100644 blossom5-v2.05.src/timer.h diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..dae9299 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.o +*.dat diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..a64e9f3 --- /dev/null +++ b/Makefile @@ -0,0 +1,24 @@ +BLOSSOM_DIR := blossom5-v2.05.src +BLOSSOM_DIRS := $(BLOSSOM_DIR) $(BLOSSOM_DIR)/MinCost $(BLOSSOM_DIR)/GEOM + +BLOSSOM_SOURCES := $(filter-out $(BLOSSOM_DIR)/example.cpp, $(foreach dir, $(BLOSSOM_DIRS), $(wildcard $(dir)/*.cpp))) +BLOSSOM_OBJS := $(patsubst %.cpp, %.o, $(BLOSSOM_SOURCES)) + +CFLAGS := -flto -O3 -D_NDEBUG -DPERFECT_MATCHING_DOUBLE +CXX = clang++ +LD = ld.lld +LIBS := -lrt + +all: rbmp + +rbmp: rbmp.cpp $(BLOSSOM_DIR)/blossom5.o + $(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o rbmp.cpp -o $@ + +$(BLOSSOM_DIR)/blossom5.o: ${BLOSSOM_OBJS} + $(LD) -r -o $@ ${BLOSSOM_OBJS} + +.cpp.o: + $(CXX) $(CFLAGS) $< -c -o $@ + +clean: + rm -f ${BLOSSOM_OBJS} $(BLOSSOM_DIR)/blossom5.o rbmp diff --git a/blossom5-v2.05.src/CHANGES.TXT b/blossom5-v2.05.src/CHANGES.TXT new file mode 100644 index 0000000..61320a6 --- /dev/null +++ b/blossom5-v2.05.src/CHANGES.TXT @@ -0,0 +1,45 @@ +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). + + + 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 +#include +#include +#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; pnext[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; pAddPerfectMatching(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; pnext[dir]) + { + nodes[e->head[dir]].is_marked = 1; + } + + kd_tree->ComputeKNN(p, K, neighbors); + for (k=0; knext[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 +#include +#include +#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(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) + { + 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=0 && _jNew(); + 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 || 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 +#include +#include +#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 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= 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 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; dsums; skip_array[d] = 1; } + + for (d=0; d(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) || (DIMd = -num0; + for (k=0; kpoints[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=0 && nodes[k].first_child) UPDATE_NODE_PTR(nodes[k].first_child); + } + for (k=0; kparent) ? ((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 split) break; + } + + //////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////// + + ptr = mapping + d0*point_num; + for (k=start; kd = 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; pDist(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; kd; 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; qnodes[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; korder > 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 +#include +#include +#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; pAddEdge(p, q, Dist(p, q)); + } + } + pm->options = options; + pm->Solve(); + for (p=0; pGetSolution(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 || iterStartUpdate(); + for (p=0; pGetTwiceSum(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; pnext[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 +#include +#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* 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= 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 diff --git a/blossom5-v2.05.src/GRAPH1.TXT b/blossom5-v2.05.src/GRAPH1.TXT new file mode 100644 index 0000000..13ec532 --- /dev/null +++ b/blossom5-v2.05.src/GRAPH1.TXT @@ -0,0 +1,10 @@ +6 9 +0 1 3 +0 3 10 +0 4 7 +1 2 -1 +1 4 5 +1 5 4 +2 5 -7 +3 4 0 +4 5 4 diff --git a/blossom5-v2.05.src/GRAPH2.TXT b/blossom5-v2.05.src/GRAPH2.TXT new file mode 100644 index 0000000..b3e487f --- /dev/null +++ b/blossom5-v2.05.src/GRAPH2.TXT @@ -0,0 +1,11 @@ +c example of a graph in DIMACS format +p edge 6 9 +e 1 2 3 +e 1 4 10 +e 1 5 7 +e 2 3 -1 +e 2 5 5 +e 2 6 4 +e 3 6 -7 +e 4 5 0 +e 5 6 4 diff --git a/blossom5-v2.05.src/LCA.h b/blossom5-v2.05.src/LCA.h new file mode 100644 index 0000000..3426f7e --- /dev/null +++ b/blossom5-v2.05.src/LCA.h @@ -0,0 +1,297 @@ +/* + +LCA.h - Imlementation of LCA data structure described in + +O. Berkman, U. Vishkin: Recursive Star-Tree Parallel Data Structure. SIAM J. Comput. 22(2): 221-242 (1993) + +It takes O(n \log n) space (and O(n \log n) time for construction). The memory requirement can be +reduced to O(n) by defining LCA_BLOCKS. The implementation then follows the approach described in + +J. Fischer, V. Heun: Theoretical and Practical Improvements on the RMQ-Problem, with Applications to LCA and LCE. + Proceedings of the 17th Annual Symposium on Combinatorial Pattern Matching + (CPM'06), Lecture Notes in Computer Science 4009, 36-48, Springer-Verlag, 2006. + +(although in-block queries are computed naively in O(K) worst-case time, where K is the size of the block, +rather than as described by Fischer and Heun). + +Written by Vladimir Kolmogorov, vnk@adastral.ucl.ac.uk. + +// example for +// 4 +// / \ +// 2 3 +// / \ +// 0 1 +// (Note: the ordering of nodes is in the tree preorder!!!) + + LCATree* lca = new LCATree(5); + char A0, A1, A2, A3, A4; + + lca->Add(&A0, &A2); + lca->Add(&A1, &A2); + lca->Add(&A2, &A4); + lca->Add(&A3, &A4); + lca->AddRoot(&A4); + + int result = lca->GetLCA(1, 3); // should be 4 + delete lca; + + +*/ + +#ifndef GNAKDLATHJSTHAJSRNAKSJDA +#define GNAKDLATHJSTHAJSRNAKSJDA + +#include +#include +#include +#include + +//#define LCA_BLOCKS + +class LCATree +{ +public: + typedef void* NodeId; // can be any type, e.g. int. (The code checks NodeId's only for equalities.) + typedef int PreorderId; + + LCATree(int node_num_max); + ~LCATree(); + + // construct tree. Nodes must be added in the tree preorder!!! + // First call returns 0, second returns 1, and so on. + PreorderId Add(NodeId i, NodeId i_parent); + PreorderId AddRoot(NodeId i); // completes tree construction + + PreorderId GetLCA(PreorderId i, PreorderId j); + // Let i0=i, j0=j be the input nodes, and let r = LCA(i0,j0). + // This function sets i and j to be the immediate children of r + // such that i is a descendant of i0 and j is a descendant of j0. + // There must hold i0!=r, j0!=r. + void GetPenultimateNodes(PreorderId& i, PreorderId& j); + + +////////////////////////////////////////////////////////////////////////// +private: + int n, n_max, K, k_max; + int** array; + + NodeId* buf0; + int* buf1; + NodeId* parent_current; + int* child_current; + + int* parents; + + int _GetLCA(int i, int j); // same as GetLCA, but assumes that i 0) { K ++; n /= 2; } + if (K < 1) K = 1; +#else + K = 1; + n = 0; +#endif + parents = new int[n_max]; + buf0 = new NodeId[n_max]; + buf1 = new int[n_max]; + parent_current = buf0; + child_current = buf1; +} + +inline LCATree::~LCATree() +{ + int k; + delete [] parents; + if (buf0) delete [] buf0; + if (buf1) delete [] buf1; + if (array) + { + for (k=1; k<=k_max; k++) delete [] array[k]; + delete [] array; + } +} + +inline LCATree::PreorderId LCATree::Add(NodeId i, NodeId i_parent) +{ + assert(n < n_max); + + if (n == 0) + { + *parent_current = i; + *(++ parent_current) = i_parent; + parents[0] = -1; + } + else + { + if (i == *parent_current) + { + int c = *child_current --; + while ( 1 ) + { + int c_next = parents[c]; + parents[c] = n; + if (c_next < 0) break; + c = c_next; + } + parent_current --; + } + if (i_parent == *parent_current) parents[n] = *child_current; + else + { + *(++ parent_current) = i_parent; + parents[n] = -1; + child_current ++; + } + } + *child_current = n; + return n ++; +} + + +inline LCATree::PreorderId LCATree::AddRoot(NodeId i) +{ + assert(n < n_max); + + if (n > 0) + { + if (i != *parent_current || parent_current != buf0+1) + { + printf("Error in LCATree construction: wrong sequence of calls!\n"); + exit(1); + } + int c = *child_current --; + while ( 1 ) + { + int c_next = parents[c]; + parents[c] = n; + if (c_next < 0) break; + c = c_next; + } + child_current ++; + } + parents[n++] = -1; + + delete [] buf0; + buf0 = NULL; + delete [] buf1; + buf1 = NULL; + + // initialize array + int b, k = 1, block_num = (n-1)/K+1; + if (block_num < 3) return n-1; + int d = (block_num-1)/4; + while (d) { k ++; d >>= 1; } + k_max = k; + + array = new int*[k_max+1]; + array[0] = parents; + for (k=1, d=2; k<=k_max; k++, d*=2) + { + array[k] = new int[block_num-d]; + if (k == 1) + { + for (b=0; b j) ? i : j; + } + } + } + + return n-1; +} + +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// + +inline int LCATree::GetLCADirect(int i, int j) +{ + while (i < j) i = parents[i]; + return i; +} + + +inline int LCATree::_GetLCA(int i, int j) +{ +#ifdef LCA_BLOCKS + + int bi = i/K, bj = j/K; + if (bi == bj) return GetLCADirect(i, j); + int i_last = (bi+1)*K-1, j_first = bj*K; + i = GetLCADirect(i, i_last); + j = GetLCADirect(j_first, j); + if (i < j) i = j; + // set j = LCA(i_last, j_first) + if (j_first - i_last == 1) j = parents[i_last]; + else + { + int k = 1, d = (bj-bi)/4; + while (d) { k ++; d >>= 1; } + int diff = 1<bj-diff); + j = (array[k][bi] > array[k][bj-diff]) ? array[k][bi] : array[k][bj-diff]; + } + return (i > j) ? i : j; + +#else + + if (j == i) return i; + + int k = 0, d = (j-i)/2; + while (d) { k ++; d >>= 1; } + int diff = 1<j-diff); + return (array[k][i] > array[k][j-diff]) ? array[k][i] : array[k][j-diff]; + +#endif +} + +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// + +inline LCATree::PreorderId LCATree::GetLCA(PreorderId i, PreorderId j) +{ + if (i > j) { PreorderId k = i; i = j; j = k; } + return _GetLCA(i, j); +} + +inline void LCATree::GetPenultimateNodes(PreorderId& _i, PreorderId& _j) +{ + int i, j, d, swap; + if (_i < _j) { i = _i; j = _j; swap = 0; } + else { i = _j; j = _i; swap = 1; } + int r = _GetLCA(i, j); + assert(i!=r && j!=r); + while (parents[i] != r) + { + int i0 = parents[i]; + d = (j - i0)/2; + while ( (i=_GetLCA(i0, i0+d)) == r ) d /= 2; + } + while (parents[j] != r) + { + int j0 = parents[j]; + d = (r - j0)/2; + while ( (j=_GetLCA(j0, j0+d)) == r ) d /= 2; + } + if (swap == 0) { _i = i; _j = j; } + else { _j = i; _i = j; } +} + +#endif diff --git a/blossom5-v2.05.src/LICENSE.TXT b/blossom5-v2.05.src/LICENSE.TXT new file mode 100644 index 0000000..346e2db --- /dev/null +++ b/blossom5-v2.05.src/LICENSE.TXT @@ -0,0 +1,25 @@ + Copyright 2008-2009 UCL Business PLC, Author Vladimir Kolmogorov (vnk@ist.ac.at) + + This software can be used for evaluation and non-commercial research purposes only. Commercial use is prohibited. + Public redistribution of the code or its derivatives is prohibited. + If you use this software for research purposes, you should cite the following paper in any resulting publication: + + Vladimir Kolmogorov. "Blossom V: A new implementation of a minimum cost perfect matching algorithm." + In Mathematical Programming Computation (MPC), July 2009, 1(1):43-67. + + For commercial use of the software not covered by this agreement, you may obtain a licence from + the copyright holders UCL Business via their licensing site: www.e-lucid.com/i/software/optimisation_software/BlossomV.html. + + 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. + + diff --git a/blossom5-v2.05.src/Makefile b/blossom5-v2.05.src/Makefile new file mode 100644 index 0000000..da30258 --- /dev/null +++ b/blossom5-v2.05.src/Makefile @@ -0,0 +1,26 @@ +DIRS := . MinCost GEOM + +SOURCES := $(foreach dir, $(DIRS), $(wildcard $(dir)/*.cpp)) +OBJS := $(patsubst %.cpp, %.o, $(SOURCES)) + +CFLAGS := -O3 -D_NDEBUG +CXX ?= c++ +LIBS := +INCLUDES := +LIBDIR := + +# Add librt if the target platform is not Darwin (OS X) +ifneq ($(shell uname -s),Darwin) + LIBS += -lrt +endif + +all: blossom5 + +blossom5: ${OBJS} + $(CXX) $(CFLAGS) ${LIBDIR} -o $@ ${OBJS} ${LIBS} + +.cpp.o: + $(CXX) $(CFLAGS) ${INCLUDES} $< -c -o $@ + +clean: + rm -f ${OBJS} blossom5 diff --git a/blossom5-v2.05.src/MinCost/MinCost.cpp b/blossom5-v2.05.src/MinCost/MinCost.cpp new file mode 100644 index 0000000..a0a9759 --- /dev/null +++ b/blossom5-v2.05.src/MinCost/MinCost.cpp @@ -0,0 +1,287 @@ +#include +#include +#include +#include "MinCost.h" + + +template + MinCost::MinCost(int _nodeNum, int _edgeNumMax, void (*err_function)(const char *)) + : nodeNum(_nodeNum), + edgeNum(0), + edgeNumMax(_edgeNumMax), + counter(0), + cost(0), + error_function(err_function) +{ + nodes = (Node*) malloc(nodeNum*sizeof(Node)); + arcs = (Arc*) malloc(2*edgeNumMax*sizeof(Arc)); + if (!nodes || !arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + + memset(nodes, 0, nodeNum*sizeof(Node)); + memset(arcs, 0, 2*edgeNumMax*sizeof(Arc)); + firstActive = &nodes[nodeNum]; +#ifdef MINCOST_DEBUG + for (int i=0; i + MinCost::~MinCost() +{ + free(nodes); + free(arcs); +} + +template + void MinCost::Init() +{ + Node* i; + Arc* a; + + for (a=arcs; ar_cap > 0 && a->GetRCost() < 0) PushFlow(a, a->r_cap); + } + + Node** lastActivePtr = &firstActive; + for (i=nodes; iexcess > 0) + { + *lastActivePtr = i; + lastActivePtr = &i->next; + } + else i->next = NULL; + } + *lastActivePtr = &nodes[nodeNum]; +} + + +template + FlowType MinCost::Augment(Node* start, Node* end) +{ + FlowType delta = (start->excess < -end->excess) ? start->excess : -end->excess; + Arc* a; + + for (a=end->parent; a; a=a->sister->head->parent) + { + if (delta > a->r_cap) delta = a->r_cap; + } + assert(delta > 0); + + end->excess += delta; + for (a=end->parent; a; a=a->head->parent) + { + DecreaseRCap(a, delta); + a = a->sister; + IncreaseRCap(a, delta); + } + start->excess -= delta; + + return delta; +} + +template + void MinCost::Dijkstra(Node* start) +{ + assert(start->excess > 0); + + Node* i; + Node* j; + Arc* a; + CostType d; + Node* permanentNodes; + + int FLAG0 = ++ counter; // permanently labeled nodes + int FLAG1 = ++ counter; // temporarily labeled nodes + + start->parent = NULL; + start->flag = FLAG1; + queue.Reset(); + queue.Add(start, 0); + + permanentNodes = NULL; + + while ( (i=queue.RemoveMin(d)) ) + { + if (i->excess < 0) + { + FlowType delta = Augment(start, i); + cost += delta*(d - i->pi + start->pi); + for (i=permanentNodes; i; i=i->next_permanent) i->pi += d; + break; + } + + i->pi -= d; + i->flag = FLAG0; + i->next_permanent = permanentNodes; + permanentNodes = i; + + for (a=i->firstNonsaturated; a; a=a->next) + { + j = a->head; + if (j->flag == FLAG0) continue; + d = a->GetRCost(); + if (j->flag == FLAG1) + { + if (d >= queue.GetKey(j)) continue; + queue.DecreaseKey(j, d); + } + else + { + queue.Add(j, d); + j->flag = FLAG1; + } + j->parent = a; + } + + } +} + + +template + CostType MinCost::Solve() +{ + Node* i; + //Init(); + while ( 1 ) + { + i = firstActive; + if (i == &nodes[nodeNum]) break; + firstActive = i->next; + i->next = NULL; + if (i->excess > 0) + { + Dijkstra(i); + if (i->excess > 0 && !i->next) + { + i->next = firstActive; + firstActive = i; + } + } + } +#ifdef MINCOST_DEBUG + TestOptimality(); + TestCosts(); +#endif + + return cost; +} + + +template + void MinCost::TestOptimality() +{ + Node* i; + Arc* a; + + for (i=nodes; iexcess != 0) + { + assert(0); + } + for (a=i->firstSaturated; a; a=a->next) + { + if (a->r_cap != 0) + { + assert(0); + } + } + for (a=i->firstNonsaturated; a; a=a->next) + { + CostType c = a->GetRCost(); + if (a->r_cap <= 0 || a->GetRCost() < -1e-5) + { + assert(0); + } + } + } +} + +#ifdef MINCOST_DEBUG + +template + void MinCost::TestCosts() +{ + Arc* a; + + CostType _cost = 0; + + for (a=arcs; ar_cap + a->sister->r_cap == a->cap_orig + a->sister->cap_orig); + _cost += a->cost*(a->cap_orig - a->r_cap); + } + + CostType delta = cost - _cost; + if (delta < 0) delta = -delta; + if (delta >= 1e-5) + { + assert(0); + } +} + +#endif + + + + +/////////////////////////////////////////////////////////////////////////////////////// + +#define FLOW_INFTY ((int)0x00fffffff) + +template + DualMinCost::DualMinCost(int _nodeNum, int _edgeNumMax) + : MinCost(_nodeNum+1, _edgeNumMax+2*_nodeNum) +{ + source = _nodeNum; +} + +template + DualMinCost::~DualMinCost() +{ +} + +template + void DualMinCost::AddUnaryTerm(NodeId i, int objective_coef) +{ + MinCost::AddNodeExcess(i, objective_coef); + MinCost::AddNodeExcess(source, -objective_coef); +} + +template + void DualMinCost::SetLowerBound(NodeId i, CostType cmin) +{ + DualMinCost::AddEdge(i, source, FLOW_INFTY, 0, -cmin); +} + +template + void DualMinCost::SetUpperBound(NodeId i, CostType cmax) +{ + DualMinCost::AddEdge(source, i, FLOW_INFTY, 0, cmax); +} + +template + void DualMinCost::AddConstraint(NodeId i, NodeId j, CostType cmax) +{ + DualMinCost::AddEdge(i, j, FLOW_INFTY, 0, cmax); +} + +template + void DualMinCost::Solve() +{ + MinCost::Solve(); +} + +template + CostType DualMinCost::GetSolution(NodeId i) +{ + return MinCost::nodes[source].pi - MinCost::nodes[i].pi; +} + + + +#include "instances.inc" + + diff --git a/blossom5-v2.05.src/MinCost/MinCost.cpp.orig b/blossom5-v2.05.src/MinCost/MinCost.cpp.orig new file mode 100644 index 0000000..43d94bd --- /dev/null +++ b/blossom5-v2.05.src/MinCost/MinCost.cpp.orig @@ -0,0 +1,287 @@ +#include +#include +#include +#include "MinCost.h" + + +template + MinCost::MinCost(int _nodeNum, int _edgeNumMax, void (*err_function)(const char *)) + : nodeNum(_nodeNum), + edgeNum(0), + edgeNumMax(_edgeNumMax), + counter(0), + cost(0), + error_function(err_function) +{ + nodes = (Node*) malloc(nodeNum*sizeof(Node)); + arcs = (Arc*) malloc(2*edgeNumMax*sizeof(Arc)); + if (!nodes || !arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + + memset(nodes, 0, nodeNum*sizeof(Node)); + memset(arcs, 0, 2*edgeNumMax*sizeof(Arc)); + firstActive = &nodes[nodeNum]; +#ifdef MINCOST_DEBUG + for (int i=0; i + MinCost::~MinCost() +{ + free(nodes); + free(arcs); +} + +template + void MinCost::Init() +{ + Node* i; + Arc* a; + + for (a=arcs; ar_cap > 0 && a->GetRCost() < 0) PushFlow(a, a->r_cap); + } + + Node** lastActivePtr = &firstActive; + for (i=nodes; iexcess > 0) + { + *lastActivePtr = i; + lastActivePtr = &i->next; + } + else i->next = NULL; + } + *lastActivePtr = &nodes[nodeNum]; +} + + +template + FlowType MinCost::Augment(Node* start, Node* end) +{ + FlowType delta = (start->excess < -end->excess) ? start->excess : -end->excess; + Arc* a; + + for (a=end->parent; a; a=a->sister->head->parent) + { + if (delta > a->r_cap) delta = a->r_cap; + } + assert(delta > 0); + + end->excess += delta; + for (a=end->parent; a; a=a->head->parent) + { + DecreaseRCap(a, delta); + a = a->sister; + IncreaseRCap(a, delta); + } + start->excess -= delta; + + return delta; +} + +template + void MinCost::Dijkstra(Node* start) +{ + assert(start->excess > 0); + + Node* i; + Node* j; + Arc* a; + CostType d; + Node* permanentNodes; + + int FLAG0 = ++ counter; // permanently labeled nodes + int FLAG1 = ++ counter; // temporarily labeled nodes + + start->parent = NULL; + start->flag = FLAG1; + queue.Reset(); + queue.Add(start, 0); + + permanentNodes = NULL; + + while ( (i=queue.RemoveMin(d)) ) + { + if (i->excess < 0) + { + FlowType delta = Augment(start, i); + cost += delta*(d - i->pi + start->pi); + for (i=permanentNodes; i; i=i->next_permanent) i->pi += d; + break; + } + + i->pi -= d; + i->flag = FLAG0; + i->next_permanent = permanentNodes; + permanentNodes = i; + + for (a=i->firstNonsaturated; a; a=a->next) + { + j = a->head; + if (j->flag == FLAG0) continue; + d = a->GetRCost(); + if (j->flag == FLAG1) + { + if (d >= queue.GetKey(j)) continue; + queue.DecreaseKey(j, d); + } + else + { + queue.Add(j, d); + j->flag = FLAG1; + } + j->parent = a; + } + + } +} + + +template + CostType MinCost::Solve() +{ + Node* i; + //Init(); + while ( 1 ) + { + i = firstActive; + if (i == &nodes[nodeNum]) break; + firstActive = i->next; + i->next = NULL; + if (i->excess > 0) + { + Dijkstra(i); + if (i->excess > 0 && !i->next) + { + i->next = firstActive; + firstActive = i; + } + } + } +#ifdef MINCOST_DEBUG + TestOptimality(); + TestCosts(); +#endif + + return cost; +} + + +template + void MinCost::TestOptimality() +{ + Node* i; + Arc* a; + + for (i=nodes; iexcess != 0) + { + assert(0); + } + for (a=i->firstSaturated; a; a=a->next) + { + if (a->r_cap != 0) + { + assert(0); + } + } + for (a=i->firstNonsaturated; a; a=a->next) + { + CostType c = a->GetRCost(); + if (a->r_cap <= 0 || a->GetRCost() < -1e-5) + { + assert(0); + } + } + } +} + +#ifdef MINCOST_DEBUG + +template + void MinCost::TestCosts() +{ + Arc* a; + + CostType _cost = 0; + + for (a=arcs; ar_cap + a->sister->r_cap == a->cap_orig + a->sister->cap_orig); + _cost += a->cost*(a->cap_orig - a->r_cap); + } + + CostType delta = cost - _cost; + if (delta < 0) delta = -delta; + if (delta >= 1e-5) + { + assert(0); + } +} + +#endif + + + + +/////////////////////////////////////////////////////////////////////////////////////// + +#define FLOW_INFTY ((int)0x00fffffff) + +template + DualMinCost::DualMinCost(int _nodeNum, int _edgeNumMax) + : MinCost(_nodeNum+1, _edgeNumMax+2*_nodeNum) +{ + source = _nodeNum; +} + +template + DualMinCost::~DualMinCost() +{ +} + +template + void DualMinCost::AddUnaryTerm(NodeId i, int objective_coef) +{ + MinCost::AddNodeExcess(i, objective_coef); + MinCost::AddNodeExcess(source, -objective_coef); +} + +template + void DualMinCost::SetLowerBound(NodeId i, CostType cmin) +{ + AddEdge(i, source, FLOW_INFTY, 0, -cmin); +} + +template + void DualMinCost::SetUpperBound(NodeId i, CostType cmax) +{ + AddEdge(source, i, FLOW_INFTY, 0, cmax); +} + +template + void DualMinCost::AddConstraint(NodeId i, NodeId j, CostType cmax) +{ + AddEdge(i, j, FLOW_INFTY, 0, cmax); +} + +template + void DualMinCost::Solve() +{ + MinCost::Solve(); +} + +template + CostType DualMinCost::GetSolution(NodeId i) +{ + return MinCost::nodes[source].pi - MinCost::nodes[i].pi; +} + + + +#include "instances.inc" + + diff --git a/blossom5-v2.05.src/MinCost/MinCost.cpp.rej b/blossom5-v2.05.src/MinCost/MinCost.cpp.rej new file mode 100644 index 0000000..1597c49 --- /dev/null +++ b/blossom5-v2.05.src/MinCost/MinCost.cpp.rej @@ -0,0 +1,25 @@ +--- MinCost.cpp ++++ MinCost.cpp +@@ -253,19 +253,19 @@ + template + void DualMinCost::SetLowerBound(NodeId i, CostType cmin) + { +- AddEdge(i, source, FLOW_INFTY, 0, -cmin); ++ DualMinCost::AddEdge(i, source, FLOW_INFTY, 0, -cmin); + } + + template + void DualMinCost::SetUpperBound(NodeId i, CostType cmax) + { +- AddEdge(source, i, FLOW_INFTY, 0, cmax); ++ DualMinCost::AddEdge(source, i, FLOW_INFTY, 0, cmax); + } + + template + void DualMinCost::AddConstraint(NodeId i, NodeId j, CostType cmax) + { +- AddEdge(i, j, FLOW_INFTY, 0, cmax); ++ DualMinCost::AddEdge(i, j, FLOW_INFTY, 0, cmax); + } + + template diff --git a/blossom5-v2.05.src/MinCost/MinCost.h b/blossom5-v2.05.src/MinCost/MinCost.h new file mode 100644 index 0000000..fa3bc91 --- /dev/null +++ b/blossom5-v2.05.src/MinCost/MinCost.h @@ -0,0 +1,504 @@ +/* + MinCost.h - successive shortest path algorithm of Ford and Fulkerson for solving a minimum cost flow problem + + 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 __MINCOST_H__ +#define __MINCOST_H__ + +#include +#include + +// if GRAPH_ASSERT is defined then all calls to graph construction functions are assert'ed for correctness +// (e.g. that node_id's are valid id's and edge capacities are non-negative). +//#define GRAPH_ASSERT + +//#define MINCOST_DEBUG + + + +template class MinCost +{ +public: + typedef int NodeId; + typedef int EdgeId; + + MinCost(int NodeNum, int edgeNumMax, void (*err_function)(const char *) = NULL); + + // Destructor + ~MinCost(); + + void AddNodeExcess(NodeId i, FlowType excess); + + // first call returns 0, second 1, and so on. + // cap, rev_cap must be non-negative. + // cost can be negative. + EdgeId AddEdge(NodeId i, NodeId j, FlowType cap, FlowType rev_cap, CostType cost); + + CostType Solve(); + + /////////////////////////////////////////////////// + + FlowType GetRCap(EdgeId e); + void SetRCap(EdgeId e, FlowType new_rcap); + FlowType GetReverseRCap(EdgeId e); + void SetReverseRCap(EdgeId e, FlowType new_rcap); + void PushFlow(EdgeId e, FlowType delta); + void UpdateCost(EdgeId e, FlowType cap_orig, CostType delta); + + CostType GetDual(NodeId i) { return nodes[i].pi; } + +///////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////// + +protected: + // internal variables and functions + + struct Node; + struct Arc; + + struct Node + { + Arc *firstNonsaturated; + Arc *firstSaturated; + + Arc *parent; + Node *next; // list of nodes with positive excesses + + FlowType excess; + CostType pi; + int flag; + union + { + int heap_ptr; + Node* next_permanent; + }; +#ifdef MINCOST_DEBUG + int id; +#endif + }; + + struct Arc + { + Node *head; + Arc *prev; + Arc *next; + Arc *sister; // reverse arc + + FlowType r_cap; // residual capacity +#ifdef MINCOST_DEBUG + FlowType cap_orig; +#endif + CostType cost; + CostType GetRCost() { return cost + head->pi - sister->head->pi; } + }; + + int nodeNum, edgeNum, edgeNumMax; + Node *nodes; + Arc *arcs; + Node* firstActive; + int counter; + CostType cost; + + + void (*error_function)(const char *); // this function is called if a error occurs, + // with a corresponding error message + // (or exit(1) is called if it's NULL) + + ///////////////////////////////////////////////////////////////////////// + + struct PriorityQueue + { + PriorityQueue(); + ~PriorityQueue(); + void Reset(); + CostType GetKey(Node* i); + void Add(Node* i, CostType key); + void DecreaseKey(Node* i, CostType key); + Node* RemoveMin(CostType& key); + + private: + struct Item + { + Node* i; + CostType key; + }* array; + int N, arraySize; + void Swap(int k1, int k2); + }; + + PriorityQueue queue; + + ///////////////////////////////////////////////////////////////////////// + + void SetRCap(Arc* a, FlowType new_rcap); + void PushFlow(Arc* a, FlowType delta); + + void Init(); + void DecreaseRCap(Arc* a, FlowType delta); + void IncreaseRCap(Arc* a, FlowType delta); + FlowType Augment(Node* start, Node* end); + void Dijkstra(Node* start); + + void TestOptimality(); +#ifdef MINCOST_DEBUG + void TestCosts(); +#endif +}; + + +template class DualMinCost : private MinCost +{ +public: + typedef int NodeId; + DualMinCost(int node_num, int constraint_num_max); + ~DualMinCost(); + + void AddUnaryTerm(NodeId i, int objective_coef); + void SetLowerBound(NodeId, CostType cmin); + void SetUpperBound(NodeId, CostType cmax); + void AddConstraint(NodeId i, NodeId j, CostType cmax); // xj - xi <= cmax + + void Solve(); + CostType GetSolution(NodeId i); + +private: + NodeId source; +}; + + + + + + + + + +/////////////////////////////////////// +// Implementation - inline functions // +/////////////////////////////////////// + + + +template + inline void MinCost::AddNodeExcess(NodeId _i, FlowType excess) +{ + assert(_i>=0 && _i 0 && !nodes[_i].next) + { + nodes[_i].next = firstActive; + firstActive = &nodes[_i]; + } +} + +template + inline typename MinCost::EdgeId MinCost::AddEdge(NodeId _i, NodeId _j, FlowType cap, FlowType rev_cap, CostType cost) +{ + assert(_i>=0 && _i=0 && _j= 0); + assert(rev_cap >= 0); + + Arc *a = &arcs[2*edgeNum]; + Arc *a_rev = a+1; + edgeNum ++; + + Node* i = nodes + _i; + Node* j = nodes + _j; + + a -> sister = a_rev; + a_rev -> sister = a; + if (cap > 0) + { + if (i->firstNonsaturated) i->firstNonsaturated->prev = a; + a -> next = i -> firstNonsaturated; + i -> firstNonsaturated = a; + } + else + { + if (i->firstSaturated) i->firstSaturated->prev = a; + a -> next = i -> firstSaturated; + i -> firstSaturated = a; + } + a->prev = NULL; + if (rev_cap > 0) + { + if (j->firstNonsaturated) j->firstNonsaturated->prev = a_rev; + a_rev -> next = j -> firstNonsaturated; + j -> firstNonsaturated = a_rev; + } + else + { + if (j->firstSaturated) j->firstSaturated->prev = a_rev; + a_rev -> next = j -> firstSaturated; + j -> firstSaturated = a_rev; + } + a_rev->prev = NULL; + + a -> head = j; + a_rev -> head = i; + a -> r_cap = cap; + a_rev -> r_cap = rev_cap; + a -> cost = cost; + a_rev -> cost = -cost; +#ifdef MINCOST_DEBUG + a->cap_orig = cap; + a_rev->cap_orig = rev_cap; +#endif + + if (a->r_cap > 0 && a->GetRCost() < 0) PushFlow(a, a->r_cap); + if (a_rev->r_cap > 0 && a_rev->GetRCost() < 0) PushFlow(a_rev, a_rev->r_cap); + + return edgeNum-1; +} + +/////////////////////////////////////// +/////////////////////////////////////// +/////////////////////////////////////// + +template + inline void MinCost::DecreaseRCap(Arc* a, FlowType delta) +{ + a->r_cap -= delta; + if (a->r_cap == 0) + { + Node* i = a->sister->head; + if (a->next) a->next->prev = a->prev; + if (a->prev) a->prev->next = a->next; + else i->firstNonsaturated = a->next; + a->next = i->firstSaturated; + if (a->next) a->next->prev = a; + a->prev = NULL; + i->firstSaturated = a; + } +} + +template + inline void MinCost::IncreaseRCap(Arc* a, FlowType delta) +{ + if (a->r_cap == 0) + { + Node* i = a->sister->head; + if (a->next) a->next->prev = a->prev; + if (a->prev) a->prev->next = a->next; + else i->firstSaturated = a->next; + a->next = i->firstNonsaturated; + if (a->next) a->next->prev = a; + a->prev = NULL; + i->firstNonsaturated = a; + } + a->r_cap += delta; +} + +template + inline FlowType MinCost::GetRCap(EdgeId e) +{ + Arc* a = &arcs[2*e]; + return a->r_cap; +} + +template + inline void MinCost::SetRCap(Arc* a, FlowType new_rcap) +{ + assert(new_rcap >= 0); +#ifdef MINCOST_DEBUG + a->cap_orig += new_rcap - a->r_cap; +#endif + if (a->r_cap == 0) + { + Node* i = a->sister->head; + if (a->next) a->next->prev = a->prev; + if (a->prev) a->prev->next = a->next; + else i->firstSaturated = a->next; + a->next = i->firstNonsaturated; + if (a->next) a->next->prev = a; + a->prev = NULL; + i->firstNonsaturated = a; + } + a->r_cap = new_rcap; + if (a->r_cap == 0) + { + Node* i = a->sister->head; + if (a->next) a->next->prev = a->prev; + if (a->prev) a->prev->next = a->next; + else i->firstNonsaturated = a->next; + a->next = i->firstSaturated; + if (a->next) a->next->prev = a; + a->prev = NULL; + i->firstSaturated = a; + } +} + +template + inline void MinCost::SetRCap(EdgeId e, FlowType new_rcap) +{ + SetRCap(&arcs[2*e], new_rcap); +} + +template + inline FlowType MinCost::GetReverseRCap(EdgeId e) +{ + Arc* a = &arcs[2*e+1]; + return a->r_cap; +} + +template + inline void MinCost::SetReverseRCap(EdgeId e, FlowType new_rcap) +{ + SetRCap(&arcs[2*e+1], new_rcap); +} + +template + inline void MinCost::PushFlow(Arc* a, FlowType delta) +{ + if (delta < 0) { a = a->sister; delta = -delta; } + DecreaseRCap(a, delta); + IncreaseRCap(a->sister, delta); + a->head->excess += delta; + a->sister->head->excess -= delta; + cost += delta*a->cost; + if (a->head->excess > 0 && !a->head->next) + { + a->head->next = firstActive; + firstActive = a->head; + } +} + +template + inline void MinCost::PushFlow(EdgeId e, FlowType delta) +{ + PushFlow(&arcs[2*e], delta); +} + +template + inline void MinCost::UpdateCost(EdgeId e, FlowType cap_orig, CostType delta) +{ + Arc* a = &arcs[2*e]; + cost += delta*(cap_orig-a->r_cap); + a->cost += delta; + a->sister->cost = -a->cost; + + if (a->GetRCost() > 0) a = a->sister; + if (a->r_cap > 0 && a->GetRCost() < 0) PushFlow(a, a->r_cap); +} + +/////////////////////////////////////// +/////////////////////////////////////// +/////////////////////////////////////// + +template + inline MinCost::PriorityQueue::PriorityQueue() +{ + N = 0; + arraySize = 16; + array = (Item*) malloc(arraySize*sizeof(Item)); +} + +template + inline MinCost::PriorityQueue::~PriorityQueue() +{ + free(array); +} + +template + inline void MinCost::PriorityQueue::Reset() +{ + N = 0; +} + +template + inline CostType MinCost::PriorityQueue::GetKey(Node* i) +{ + return array[i->heap_ptr].key; +} + +template + inline void MinCost::PriorityQueue::Swap(int k1, int k2) +{ + Item* a = array+k1; + Item* b = array+k2; + a->i->heap_ptr = k2; + b->i->heap_ptr = k1; + Node* i = a->i; a->i = b->i; b->i = i; + CostType key = a->key; a->key = b->key; b->key = key; +} + +template + inline void MinCost::PriorityQueue::Add(Node* i, CostType key) +{ + if (N == arraySize) + { + arraySize *= 2; + array = (Item*) realloc(array, arraySize*sizeof(Item)); + } + int k = i->heap_ptr = N ++; + array[k].i = i; + array[k].key = key; + while (k > 0) + { + int k2 = (k-1)/2; + if (array[k2].key <= array[k].key) break; + Swap(k, k2); + k = k2; + } +} + +template + inline void MinCost::PriorityQueue::DecreaseKey(Node* i, CostType key) +{ + int k = i->heap_ptr; + array[k].key = key; + while (k > 0) + { + int k2 = (k-1)/2; + if (array[k2].key <= array[k].key) break; + Swap(k, k2); + k = k2; + } +} + +template + inline typename MinCost::Node* MinCost::PriorityQueue::RemoveMin(CostType& key) +{ + if (N == 0) return NULL; + + Swap(0, N-1); + N --; + + int k = 0; + while ( 1 ) + { + int k1 = 2*k + 1, k2 = k1 + 1; + if (k1 >= N) break; + int k_min = (k2 >= N || array[k1].key <= array[k2].key) ? k1 : k2; + if (array[k].key <= array[k_min].key) break; + Swap(k, k_min); + k = k_min; + } + + key = array[N].key; + return array[N].i; +} + + +#endif diff --git a/blossom5-v2.05.src/MinCost/instances.inc b/blossom5-v2.05.src/MinCost/instances.inc new file mode 100644 index 0000000..a6d099e --- /dev/null +++ b/blossom5-v2.05.src/MinCost/instances.inc @@ -0,0 +1,12 @@ +#include "MinCost.h" + +#ifdef _MSC_VER +#pragma warning(disable: 4661) +#endif + + +template class MinCost; +template class MinCost; + +template class DualMinCost; +template class DualMinCost; diff --git a/blossom5-v2.05.src/MinCost/p b/blossom5-v2.05.src/MinCost/p new file mode 100644 index 0000000..84855b0 --- /dev/null +++ b/blossom5-v2.05.src/MinCost/p @@ -0,0 +1,26 @@ +index 3b6bedd..44338f7 100644 +--- Autotune/blossomv/MinCost/MinCost.cpp ++++ Autotune/blossomv/MinCost/MinCost.cpp +@@ -253,19 +253,19 @@ template + template + void DualMinCost::SetLowerBound(NodeId i, CostType cmin) + { +- AddEdge(i, source, FLOW_INFTY, 0, -cmin); ++ DualMinCost::AddEdge(i, source, FLOW_INFTY, 0, -cmin); + } + + template + void DualMinCost::SetUpperBound(NodeId i, CostType cmax) + { +- AddEdge(source, i, FLOW_INFTY, 0, cmax); ++ DualMinCost::AddEdge(source, i, FLOW_INFTY, 0, cmax); + } + + template + void DualMinCost::AddConstraint(NodeId i, NodeId j, CostType cmax) + { +- AddEdge(i, j, FLOW_INFTY, 0, cmax); ++ DualMinCost::AddEdge(i, j, FLOW_INFTY, 0, cmax); + } + + template diff --git a/blossom5-v2.05.src/PMduals.cpp b/blossom5-v2.05.src/PMduals.cpp new file mode 100644 index 0000000..ac3fa66 --- /dev/null +++ b/blossom5-v2.05.src/PMduals.cpp @@ -0,0 +1,441 @@ +#include +#include +#include +#include "PMimplementation.h" +#include "MinCost/MinCost.h" + + +void PerfectMatching::ComputeEpsGlobal() +{ + Node* r; + PriorityQueue::Item* q; + Tree* t; + Tree* t2; + TreeEdge* e; + int i, j, k, N = 0, E = 0; + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t = r->tree; + t->id = N; + N += 2; + for (k=0; k<2; k++) + for (e=t->first[k]; e; e=e->next[k]) E += 6; + } + DualMinCost* m = new DualMinCost(N, E); + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t = r->tree; + i = t->id; + m->AddUnaryTerm(i, -1); + m->SetLowerBound(i, 0); + m->AddUnaryTerm(i+1, 1); + m->SetUpperBound(i+1, 0); + + if (t->eps_delta < PM_INFTY) + { + m->SetUpperBound(i, t->eps_delta); + m->SetLowerBound(i+1, -t->eps_delta); + } + for (e=t->first[0]; e; e=e->next[0]) + { + t2 = e->head[0]; + if (t2 == NULL) continue; + j = e->head[0]->id; + if ((q=e->pq01[0].GetMin())) + { + m->AddConstraint(j, i, q->slack - t->eps + t2->eps); + m->AddConstraint(i+1, j+1, q->slack - t->eps + t2->eps); + } + if ((q=e->pq01[1].GetMin())) + { + m->AddConstraint(i, j, q->slack - t2->eps + t->eps); + m->AddConstraint(j+1, i+1, q->slack - t2->eps + t->eps); + } + if ((q=e->pq00.GetMin())) + { + m->AddConstraint(i+1, j, q->slack - t->eps - t2->eps); + m->AddConstraint(j+1, i, q->slack - t->eps - t2->eps); + } + } + } + m->Solve(); + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t = r->tree; + i = t->id; + t->eps_delta = (m->GetSolution(i) - m->GetSolution(i+1))/2; + } + delete m; +} + +void PerfectMatching::ComputeEpsSingle() +{ + Node* r; + PriorityQueue::Item* q; + Tree* t; + Tree* t2; + TreeEdge* e; + REAL eps = PM_INFTY; + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t = r->tree; + if (eps > t->eps_delta) eps = t->eps_delta; + for (e=t->first[0]; e; e=e->next[0]) + { + t2 = e->head[0]; + if ((q=e->pq00.GetMin()) && 2*eps > q->slack-t->eps-t2->eps) + { + eps = (q->slack-t->eps-t2->eps)/2; + } + } + } + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + r->tree->eps_delta = eps; + } +} + + +void PerfectMatching::ComputeEpsCC() +{ + Node* r; + PriorityQueue::Item* q; + Tree* t0; + Tree* t; + Tree* t2; + Tree* t_next; + TreeEdge* e; + REAL eps, eps2; + Tree* queue_last; + int dir; + Tree* FIXED_TREE = trees-1; + int component_num = 0; + TreeEdge** e_ptr; + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t0 = r->tree; + t0->next = NULL; + } + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t0 = r->tree; + if (t0->next) continue; + eps = t0->eps_delta; + + t0->next = queue_last = t = t0; + while ( 1 ) + { + for (dir=0; dir<2; dir++) + for (e_ptr=&t->first[dir], e=*e_ptr; e; e=*e_ptr) + { + t2 = e->head[dir]; + if (t2 == NULL) { *e_ptr = e->next[dir]; tree_edges->Delete(e); continue; } + e_ptr = &e->next[dir]; + + REAL eps00 = ((q=e->pq00.GetMin())) ? (q->slack - t->eps - t2->eps) : PM_INFTY; + if (t2->next && t2->next != FIXED_TREE) + { + if (2*eps > eps00) eps = eps00/2; + continue; + } + + REAL eps01[2]; + eps01[dir] = ((q=e->pq01[dir].GetMin())) ? (q->slack - t->eps + t2->eps) : PM_INFTY; + eps01[1-dir] = ((q=e->pq01[1-dir].GetMin())) ? (q->slack - t2->eps + t->eps) : PM_INFTY; + + if (t2->next == FIXED_TREE) eps2 = t2->eps_delta; + else if (eps01[0] > 0 && eps01[1] > 0) eps2 = 0; + else + { + queue_last->next = t2; + queue_last = t2; + t2->next = t2; + if (eps > eps00) eps = eps00; + if (eps > t2->eps_delta) eps = t2->eps_delta; + continue; + } + if (eps > eps00 - eps2) eps = eps00 - eps2; + if (eps > eps2 + eps01[dir]) eps = eps2 + eps01[dir]; + } + + if (t->next == t) break; + t = t->next; + } + for (t=t0; ; t=t_next) + { + t->eps_delta = eps; + t_next = t->next; + t->next = FIXED_TREE; + if (t_next == t) break; + } + component_num ++; + } + //printf("%d CCs ", component_num); +} + + +void PerfectMatching::ComputeEpsSCC() +{ + PriorityQueue::Item* q; + Node* r; + Tree* t0; + Tree* t; + Tree* t2; + TreeEdge* e; + TreeEdge** e_ptr; + REAL eps; + int c, dir; + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t0 = r->tree; + t0->dfs_parent = NULL; + + for (dir=0; dir<2; dir++) + for (e_ptr=&t0->first[dir], e=*e_ptr; e; e=*e_ptr) + { + t2 = e->head[dir]; + if (t2 == NULL) { *e_ptr = e->next[dir]; tree_edges->Delete(e); continue; } + e_ptr = &e->next[dir]; + } + } + Tree* stack = NULL; + + // first DFS + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t0 = r->tree; + if (t0->dfs_parent) continue; + t = t0; + e = (t->first[0]) ? t->first[0] : t->first[1]; + t->dfs_parent = (TreeEdge*)trees; + while ( 1 ) + { + if (e == NULL) + { + t->next = stack; + stack = t; + + if (t == t0) break; + + e = t->dfs_parent; + if (t == e->head[0]) { t = e->head[1]; e = (e->next[0]) ? e->next[0] : t->first[1]; } + else { t = e->head[0]; e = e->next[1]; } + continue; + } + + if (e->head[1] == t) + { + if (e->head[0]->dfs_parent || !(q=e->pq01[0].GetMin()) || q->slack - t->eps + e->head[0]->eps > 0) { e = (e->next[0]) ? e->next[0] : t->first[1]; continue; } + t = e->head[0]; + } + else + { + if (e->head[1]->dfs_parent || !(q=e->pq01[1].GetMin()) || q->slack - t->eps + e->head[1]->eps > 0) { e = e->next[1]; continue; } + t = e->head[1]; + } + t->dfs_parent = e; + e = (t->first[0]) ? t->first[0] : t->first[1]; + } + } + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) r->tree->dfs_parent = NULL; + + int component_num = 0; + while (stack) + { + t0 = stack; + stack = t0->next; + if (t0->dfs_parent) continue; + t = t0; + e = (t->first[0]) ? t->first[0] : t->first[1]; + t->dfs_parent = (TreeEdge*)trees; + while ( 1 ) + { + if (e == NULL) + { + e = t->dfs_parent; + t->dfs_parent = (TreeEdge*)((char*)trees + component_num); + if (t == t0) break; + if (t == e->head[0]) { t = e->head[1]; e = (e->next[0]) ? e->next[0] : t->first[1]; } + else { t = e->head[0]; e = e->next[1]; } + continue; + } + + if (e->head[1] == t) + { + if (e->head[0]->dfs_parent || !(q=e->pq01[1].GetMin()) || q->slack - e->head[0]->eps + t->eps > 0) { e = (e->next[0]) ? e->next[0] : t->first[1]; continue; } + t = e->head[0]; + } + else + { + if (e->head[1]->dfs_parent || !(q=e->pq01[0].GetMin()) || q->slack - e->head[1]->eps + t->eps > 0) { e = e->next[1]; continue; } + t = e->head[1]; + } + t->dfs_parent = e; + e = (t->first[0]) ? t->first[0] : t->first[1]; + } + component_num ++; + } + + Tree** array = new Tree*[component_num]; + memset(array, 0, component_num*sizeof(Tree*)); + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t = r->tree; + t->id = (int)((char*)t->dfs_parent - (char*)trees); + t->next = array[t->id]; + array[t->id] = t; + } + + for (c=component_num-1; c>=0; c--) + { + eps = PM_INFTY; + for (t=array[c]; t; t=t->next) + { + if (eps > t->eps_delta) eps = t->eps_delta; + FOR_ALL_TREE_EDGES(t, e, dir) + { + t2 = e->head[dir]; + REAL eps00 = (q=e->pq00.GetMin()) ? (q->slack-t->eps-t2->eps) : PM_INFTY; + REAL eps01[2]; + eps01[dir] = ((q=e->pq01[dir].GetMin())) ? (q->slack - t->eps + t2->eps) : PM_INFTY; + eps01[1-dir] = ((q=e->pq01[1-dir].GetMin())) ? (q->slack - t2->eps + t->eps) : PM_INFTY; + if (t2->id < c) + { + if (eps > eps01[dir]) eps = eps01[dir]; + if (eps > eps00) eps = eps00; + } + else if (t2->id == c) + { + if (2*eps > eps00) eps = eps00 / 2; + } + else + { + if (eps > eps01[dir] + t2->eps_delta) eps = eps01[dir] + t2->eps_delta; + if (eps > eps00 - t2->eps_delta) eps = eps00 - t2->eps_delta; + } + } + } + for (t=array[c]; t; t=t->next) t->eps_delta = eps; + } + + delete [] array; + //printf("%d SCCs ", component_num); +} + +void PerfectMatching::CommitEps() +{ + printf("CommitEps()\n"); + Node* i; + Node* j; + Node* r; + int dir; + Edge* a; + EdgeIterator I; + Tree* t; + TreeEdge* e; + TreeEdge** e_ptr; + REAL eps, eps2; + PriorityQueue::Item* q; + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + t = r->tree; + eps = t->eps; + + i = r; + while ( 1 ) + { + i->y += eps; + if (!i->is_tree_root) + { + Node* i0 = i; + i = ARC_HEAD(i0->match); + if (i->is_blossom) ARC_TO_EDGE_PTR(i0->match)->slack -= eps; + else i->y -= eps; + FOR_ALL_EDGES(i, a, dir, I) + { + GET_OUTER_HEAD(a, dir, j); + + a->slack += eps; + if (j->flag == 0) a->slack -= j->tree->eps; + } + i = i0; + } + + MOVE_NODE_IN_TREE(i); + } + + t->pq0.Update(-eps); + + PriorityQueue pq00 = t->pq00; + t->pq00.Reset(); + for (q=pq00.GetAndResetFirst(); q; q=pq00.GetAndResetNext()) + { + a = (Edge*)q; + if (ProcessEdge00(a)) t->pq00.Add(a); + } + + for (e_ptr=&t->first[0], e=*e_ptr; e; e=*e_ptr) + { + if (e->head[0] == NULL) { *e_ptr = e->next[0]; tree_edges->Delete(e); continue; } + e_ptr = &e->next[0]; + + eps2 = e->head[0]->eps; + e->pq00.Update( - eps - eps2 ); + } + } + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) r->tree->eps = 0; +} + +bool PerfectMatching::UpdateDuals() +{ + Node* r; + + double start_time = get_time(); + + //////////////////////////////////////////////////////////////////////////////////// + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + Tree* t = r->tree; + PriorityQueue::Item* q; + REAL eps = PM_INFTY; + if ((q=t->pq0.GetMin())) eps = q->slack; + if ((q=t->pq_blossoms.GetMin()) && eps > q->slack) eps = q->slack; + while ((q=t->pq00.GetMin())) + { + if (ProcessEdge00((Edge*)q, false)) break; + t->pq00.Remove(q, pq_buf); + } + if (q && 2*eps > q->slack) eps = q->slack/2; + t->eps_delta = eps - t->eps; + } + + if (tree_num >= options.dual_LP_threshold*node_num) + { + if (options.dual_greedy_update_option == 0) ComputeEpsCC(); + else if (options.dual_greedy_update_option == 1) ComputeEpsSCC(); + else ComputeEpsSingle(); + } + else ComputeEpsGlobal(); + + REAL delta = 0; + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + if (r->tree->eps_delta > 0) + { + delta += r->tree->eps_delta; + r->tree->eps += r->tree->eps_delta; + } + } + + stat.dual_time += get_time() - start_time; + + return (delta > PM_THRESHOLD); +} diff --git a/blossom5-v2.05.src/PMexpand.cpp b/blossom5-v2.05.src/PMexpand.cpp new file mode 100644 index 0000000..79ba684 --- /dev/null +++ b/blossom5-v2.05.src/PMexpand.cpp @@ -0,0 +1,304 @@ +#include +#include +#include +#include "PMimplementation.h" + + + + +inline void PerfectMatching::ProcessSelfloop(Node* b, Edge* a) +{ + int dir; + Node* j; + Node* prev[2]; + for (dir=0; dir<2; dir++) + { + j = a->head[dir]; + GET_PENULTIMATE_BLOSSOM(j); + prev[dir] = j; + } + if (prev[0] != prev[1]) + { + ADD_EDGE(prev[0], a, 1); + ADD_EDGE(prev[1], a, 0); + a->slack -= 2*prev[0]->blossom_eps; + } + else + { + a->next[0] = prev[0]->blossom_selfloops; + prev[0]->blossom_selfloops = a; + } + +} + + + +void PerfectMatching::Expand(Node* b) +{ + assert(b->is_blossom); + assert(b->is_outer); + assert(b->flag == 1); + + double start_time = get_time(); + + Node* i; + Node* j; + Node* k; + Edge* a; + EdgeIterator I; + int dir; + ExpandTmpItem* tmp_item; + Tree* t = b->tree; + REAL eps = t->eps; + Edge* a_augment = NULL; + + GET_TREE_PARENT(b, i); + a = ARC_TO_EDGE_PTR(b->tree_parent); + dir = ARC_TO_EDGE_DIR(b->tree_parent); + + j = a->head0[1-dir]; + GET_PENULTIMATE_BLOSSOM(j); + MOVE_EDGE(b, j, a, dir); + + a = ARC_TO_EDGE_PTR(b->match); + dir = ARC_TO_EDGE_DIR(b->match); + k = a->head0[1-dir]; + GET_PENULTIMATE_BLOSSOM(k); + MOVE_EDGE(b, k, a, dir); + + i = ARC_HEAD(k->blossom_sibling); + while ( 1 ) + { + tmp_item = expand_tmp_list->New(); + tmp_item->i = i; tmp_item->blossom_parent = i->blossom_parent; tmp_item->blossom_grandparent = i->blossom_grandparent; + i->flag = 2; + + // blossom_selfloops + i->is_outer = 1; + while ((a=i->blossom_selfloops)) + { + i->blossom_selfloops = a->next[0]; + ProcessSelfloop(i, a); + } + i->is_outer = 0; + + if (i == k) break; + i->match = i->blossom_sibling; + j = ARC_HEAD(i->match); + tmp_item = expand_tmp_list->New(); + tmp_item->i = j; tmp_item->blossom_parent = j->blossom_parent; tmp_item->blossom_grandparent = j->blossom_grandparent; + j->flag = 2; + + // blossom_selfloops + j->is_outer = 1; + while ((a=j->blossom_selfloops)) + { + j->blossom_selfloops = a->next[0]; + ProcessSelfloop(j, a); + } + j->is_outer = 0; + + j->match = ARC_REV(i->match); + i = ARC_HEAD(j->blossom_sibling); + } + k->match = b->match; + i = ARC_TAIL(b->tree_parent); + Arc* aa = i->blossom_sibling; + i->flag = 1; i->tree = b->tree; i->y += b->tree->eps; + i->tree_parent = b->tree_parent; + if (i != k) + { + Node** i_ptr; + if (i->match == aa) + { + i = ARC_HEAD(i->match); + i_ptr = &j; + while ( 1 ) + { + aa = i->blossom_sibling; + i->flag = 0; i->tree = b->tree; i->y -= t->eps; + *i_ptr = i; + i_ptr = &i->first_tree_child; + i->tree_sibling_prev = i; + i->tree_sibling_next = NULL; + i = ARC_HEAD(aa); + i->flag = 1; i->tree = b->tree; i->y += t->eps; + i->tree_parent = ARC_REV(aa); + if (i == k) break; + i = ARC_HEAD(i->match); + } + *i_ptr = ARC_HEAD(k->match); + } + else + { + i = k; + j = ARC_HEAD(k->match); + do + { + i->tree_parent = i->blossom_sibling; + i->flag = 1; i->tree = b->tree; i->y += b->tree->eps; + i = ARC_HEAD(i->tree_parent); + i->flag = 0; i->tree = b->tree; i->y -= b->tree->eps; + i->first_tree_child = j; + j = i; + i->tree_sibling_prev = i; + i->tree_sibling_next = NULL; + i = ARC_HEAD(i->match); + } while ( i->flag != 1 ); + } + i = ARC_HEAD(k->match); + + j->tree_sibling_prev = i->tree_sibling_prev; + j->tree_sibling_next = i->tree_sibling_next; + if (i->tree_sibling_prev->tree_sibling_next) i->tree_sibling_prev->tree_sibling_next = j; + else ARC_HEAD(b->tree_parent)->first_tree_child = j; + if (i->tree_sibling_next) i->tree_sibling_next->tree_sibling_prev = j; + else ARC_HEAD(b->tree_parent)->first_tree_child->tree_sibling_prev = j; + + i->tree_sibling_prev = i; + i->tree_sibling_next = NULL; + } + + // go through inner arcs + i = k; + while ( 1 ) + { + // "-" node + if (i->is_blossom) + { + a = ARC_TO_EDGE_PTR(i->match); + REAL tmp = a->slack; a->slack = i->y; i->y = tmp; + t->pq_blossoms.Add(a); + } + FOR_ALL_EDGES(i, a, dir, I) + { + j = a->head[dir]; + if (j->flag != 0) a->slack -= eps; + } + i->is_processed = 1; + if (i->tree_parent == b->tree_parent) break; + i = ARC_HEAD(i->tree_parent); + // "+" node + FOR_ALL_EDGES(i, a, dir, I) + { + j = a->head[dir]; + if (j->flag == 2) + { + a->slack += eps; + t->pq0.Add(a); + } + else if (j->flag == 0 && i < j) + { + a->slack += 2*eps; + t->pq00.Add(a); + } + } + i->is_processed = 1; + i = ARC_HEAD(i->match); + } + + // go through boundary arcs + for (tmp_item=expand_tmp_list->ScanFirst(); tmp_item; tmp_item=expand_tmp_list->ScanNext()) + { + i = tmp_item->i; + j = tmp_item->blossom_parent; tmp_item->blossom_parent = i->blossom_parent; i->blossom_parent = j; + j = tmp_item->blossom_grandparent; tmp_item->blossom_grandparent = i->blossom_grandparent; i->blossom_grandparent = j; + } + for (dir=0; dir<2; dir++) + { + if (!b->first[dir]) continue; + b->first[dir]->prev[dir]->next[dir] = NULL; + + Edge* a_next; + for (a=b->first[dir]; a; a=a_next) + { + a_next = a->next[dir]; + i = a->head0[1-dir]; + GET_PENULTIMATE_BLOSSOM2(i); + ADD_EDGE(i, a, dir); + GET_OUTER_HEAD(a, dir, j); + + if (i->flag == 1) continue; + + if (j->flag == 0 && j->tree != t) j->tree->pq_current->pq01[1-j->tree->dir_current].Remove(a, pq_buf); + + if (i->flag == 2) + { + a->slack += eps; + if (j->flag == 0) j->tree->pq0.Add(a); + } + else + { + a->slack += 2*eps; + if (j->flag == 2) t->pq0.Add(a); + else if (j->flag == 0) + { + if (j->tree != t) + { + if (!j->tree->pq_current) AddTreeEdge(t, j->tree); + if (a->slack <= j->tree->eps + eps) a_augment = a; + } + j->tree->pq_current->pq00.Add(a); + } + else if (j->tree != t) + { + if (!j->tree->pq_current) AddTreeEdge(t, j->tree); + j->tree->pq_current->pq01[j->tree->dir_current].Add(a); + } + + } + } + } + for (tmp_item=expand_tmp_list->ScanFirst(); tmp_item; tmp_item=expand_tmp_list->ScanNext()) + { + i = tmp_item->i; + i->blossom_parent = tmp_item->blossom_parent; + i->blossom_grandparent = tmp_item->blossom_grandparent; + i->is_outer = 1; + } + expand_tmp_list->Reset(); + + b->is_removed = 1; + b->tree_sibling_next = removed_first; + removed_first = b; + removed_num ++; + if (4*removed_num > node_num) FreeRemoved(); + + blossom_num --; + stat.expand_count ++; + + stat.expand_time += get_time() - start_time; + + if (a_augment) Augment(a_augment); +} + +void PerfectMatching::FreeRemoved() +{ + Node* i0; + Node* i; + for (i0=nodes; i0is_outer && !i->is_marked; i=i->blossom_parent) + { + i->is_marked = 1; + if (i->blossom_grandparent->is_removed) i->blossom_grandparent = i->blossom_parent; + } + } + for (i0=nodes; i0is_outer && i->is_marked; i=i->blossom_parent) + { + i->is_marked = 0; + } + } + + while ((i=removed_first)) + { + removed_first = i->tree_sibling_next; + blossoms->Delete(i); + removed_num --; + } + + assert(removed_num == 0); +} + diff --git a/blossom5-v2.05.src/PMimplementation.h b/blossom5-v2.05.src/PMimplementation.h new file mode 100644 index 0000000..b949c14 --- /dev/null +++ b/blossom5-v2.05.src/PMimplementation.h @@ -0,0 +1,356 @@ +/* + PMimplementation.h + + Copyright 2008 Vladimir Kolmogorov (vnk@ist.ac.at) + + 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 ASKHAKJSNTJAKSNBBAVASRA +#define ASKHAKJSNTJAKSNBBAVASRA + +#ifdef _MSC_VER +#pragma warning(disable: 4311) +#pragma warning(disable: 4312) +#endif + +#include "PerfectMatching.h" +#include "timer.h" +#include "PQ.h" +#include "LCA.h" + + +#define LCA_REPAIRS + + + + +#define IS_INT ( ((REAL)1 / 2) == 0 ) +#define COST_FACTOR 2 + +#define PM_THRESHOLD ((REAL)1e-12) + +struct PerfectMatching::Node +{ + unsigned int is_outer : 1; // 0 - the node is contained in another blossom, 1 - it's the outermost node + unsigned int flag : 2; // 0 corresponds to +, 1 corresponds to -, 2 corresponds to a free node + unsigned int is_tree_root : 1; + unsigned int is_processed : 1; + unsigned int is_blossom : 1; + unsigned int is_marked : 1; + unsigned int is_removed : 1; + + Edge* first[2]; + union + { + Arc* match; // used if not a tree root (is_tree_root = 0) or it's an inner node (is_outer == 0) + Node* blossom_grandparent; + }; + REAL y; + + union + { + struct // used when is_outer = 0 + { + Arc* blossom_sibling; + Node* blossom_parent; + union + { + Edge* blossom_selfloops; + Node* blossom_ptr; // used in repairs +#ifdef LCA_REPAIRS + int lca_preorder; // used in repairs +#endif + }; + REAL blossom_eps; // stores 'eps' of the tree at the moment when the node was shrunk into a blossom. + // (it is used for determining slacks of self-loops) + }; + struct // used when is_outer = 1 + { + union + { + struct // used for "+" nodes (flag = 0) + { + Node* first_tree_child; + Node* tree_sibling_prev; // circular list (with one exception: parent->first_tree_child->tree_sibling_prev->tree_sibling_next is NULL) + Node* tree_sibling_next; + }; + Arc* tree_parent; // used for "-" nodes (flag = 1) + }; + union + { + Tree* tree; + Edge* best_edge; // used during InitGlobal() for non-tree nodes +#ifdef LCA_REPAIRS + int lca_size; // used in repairs + LCATreeX* lca; // used in repairs +#endif + }; + }; + }; +}; + +struct PerfectMatching::Edge : PriorityQueue::Item +{ + Node* head[2]; + Node* head0[2]; + Edge* next[2]; + Edge* prev[2]; +}; + +typedef unsigned long POINTER_TYPE; +// if the declaration below fails, set POINTER_TYPE to be the appropriate integer type of the same length as (void*) +extern char dummy_array[2*(sizeof(void*)==sizeof(POINTER_TYPE))-1]; + +#define ARC_TO_EDGE_PTR(a) ( (Edge*) ( ((POINTER_TYPE)(a)) & (~1) ) ) +#define ARC_TO_EDGE_DIR(a) ( (int) ( ((POINTER_TYPE)(a)) & 1 ) ) +#define EDGE_DIR_TO_ARC(a, dir) ( (Arc*) ( (char*)(a) + (dir)) ) + +#define ARC_REV(a) ( (Arc*) ( ((POINTER_TYPE)(a)) ^ 1 ) ) + +#define ARC_TAIL(a) (ARC_TO_EDGE_PTR(a)->head [1-ARC_TO_EDGE_DIR(a)]) +#define ARC_TAIL0(a) (ARC_TO_EDGE_PTR(a)->head0[1-ARC_TO_EDGE_DIR(a)]) +#define ARC_HEAD(a) (ARC_TO_EDGE_PTR(a)->head [ARC_TO_EDGE_DIR(a)]) +#define ARC_HEAD0(a) (ARC_TO_EDGE_PTR(a)->head0[ARC_TO_EDGE_DIR(a)]) + +struct PerfectMatching::PQPointers +{ + PriorityQueue pq00; // plus-plus edges + union + { + PriorityQueue pq01[2]; // plus-minus, minus-plus edges. Used for tree edges. + struct // used for trees. + { + PriorityQueue pq0; // plus-free edges + PriorityQueue pq_blossoms; + }; + }; +}; + +struct PerfectMatching::Tree : PQPointers +{ + REAL eps; + TreeEdge* first[2]; + Node* root; + + PQPointers* pq_current; + int dir_current; + + ///////////////////////////////////////// + // used while computing dual updates + REAL eps_delta; + Tree* next; + union + { + int id; + TreeEdge* dfs_parent; + }; +}; + +struct PerfectMatching::TreeEdge : PQPointers +{ + Tree* head[2]; + TreeEdge* next[2]; +}; + + + +/////////////////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////////////// + + +#define GET_PENULTIMATE_BLOSSOM(j)\ + {\ + Node* jtmp1 = j;\ + while ( 1 )\ + {\ + if (!j->blossom_grandparent->is_outer) j = j->blossom_grandparent;\ + else if (j->blossom_grandparent != j->blossom_parent) j->blossom_grandparent = j->blossom_parent;\ + else break;\ + }\ + Node* jtmp2;\ + for ( ; jtmp1!=j; jtmp1=jtmp2)\ + {\ + jtmp2 = jtmp1->blossom_grandparent;\ + jtmp1->blossom_grandparent = j;\ + }\ + } +#define GET_PENULTIMATE_BLOSSOM2(j)\ + {\ + Node* jtmp1 = j;\ + Node* jtmp_prev = NULL;\ + while ( 1 )\ + {\ + if (!j->blossom_grandparent->is_outer) { jtmp_prev = j; j = j->blossom_grandparent; }\ + else if (j->blossom_grandparent != j->blossom_parent) j->blossom_grandparent = j->blossom_parent;\ + else break;\ + }\ + if (jtmp_prev)\ + {\ + Node* jtmp2;\ + for ( ; jtmp1!=jtmp_prev; jtmp1=jtmp2)\ + {\ + jtmp2 = jtmp1->blossom_grandparent;\ + jtmp1->blossom_grandparent = jtmp_prev;\ + }\ + }\ + } + + +/////////////////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////////////// + + +struct PerfectMatching::EdgeIterator +{ + Edge* a_last; + int start_flag; +}; + +#define FOR_ALL_EDGES(i, a, dir, I)\ + for ( dir = (i->first[0]) ? 0 : 1, I.a_last = a = i->first[dir], I.start_flag = (a) ? 0 : 1;\ + a != I.a_last || (I.start_flag ++ == 0) || (dir ++ == 0 && (I.a_last = a = i->first[1]));\ + a = a->next[dir] ) + +#define CONTINUE_FOR_ALL_EDGES(i, a, dir, I)\ + for ( a = a->next[dir];\ + a != I.a_last || (I.start_flag ++ == 0) || (dir ++ == 0 && (I.a_last = a = i->first[1]));\ + a = a->next[dir] ) + +#define REMOVE_EDGE(i, a, dir)\ + {\ + if ((a)->prev[dir]==(a)) (i)->first[dir] = NULL;\ + else\ + {\ + (a)->prev[dir]->next[dir] = (a)->next[dir];\ + (a)->next[dir]->prev[dir] = (a)->prev[dir];\ + (i)->first[dir] = (a)->next[dir];\ + }\ + } + +#define ADD_EDGE(i, a, dir)\ + {\ + if ((i)->first[dir])\ + {\ + (a)->prev[dir] = (i)->first[dir]->prev[dir];\ + (a)->next[dir] = (i)->first[dir];\ + (i)->first[dir]->prev[dir]->next[dir] = (a);\ + (i)->first[dir]->prev[dir] = (a);\ + }\ + else (i)->first[dir] = (a)->prev[dir] = (a)->next[dir] = (a);\ + (a)->head[1-(dir)] = (i);\ + } + +#define MOVE_EDGE(i_old, i_new, a, dir)\ + {\ + REMOVE_EDGE(i_old, a, dir);\ + ADD_EDGE(i_new, a, dir);\ + } + +#define GET_OUTER_HEAD(a, dir, j)\ + {\ + j = (a)->head[dir];\ + if (!j->is_outer)\ + {\ + Node* j_orig = j;\ + GET_PENULTIMATE_BLOSSOM(j);\ + j = j->blossom_parent;\ + int dir_rev = 1 - (dir);\ + MOVE_EDGE(j_orig, j, a, dir_rev);\ + }\ + } + +#define GET_TREE_PARENT(child, parent)\ + {\ + Arc* a = (child)->tree_parent;\ + Edge* e = ARC_TO_EDGE_PTR(a);\ + int dir = ARC_TO_EDGE_DIR(a);\ + GET_OUTER_HEAD(e, dir, parent);\ + } + +/////////////////////////////////////////////////////////////////////////////////////////////// + +struct PerfectMatching::TreeEdgeIterator +{ + TreeEdge** e_ptr; +}; + +#define FOR_ALL_TREE_EDGES(t, e, dir)\ + for ( dir = (t->first[0]) ? 0 : 1, e = t->first[dir];\ + e || (dir ++ == 0 && (e = t->first[1]));\ + e = e->next[dir] ) + +#define FOR_ALL_TREE_EDGES_X(t, e, dir, T)\ + for ( dir = (t->first[0]) ? 0 : 1, T.e_ptr = &t->first[dir], e = *T.e_ptr;\ + e || (dir ++ == 0 && (e = *(T.e_ptr = &t->first[1])));\ + e = *T.e_ptr )\ + if (e->head[dir] == NULL) { *T.e_ptr = e->next[dir]; tree_edges->Delete(e); }\ + else if ((T.e_ptr = &e->next[dir])) + +/////////////////////////////////////////////////////////////////////////////////////////////// + +#define MOVE_NODE_IN_TREE(i)\ + {\ + if ((i)->first_tree_child) (i) = (i)->first_tree_child;\ + else\ + {\ + while (!(i)->is_tree_root && !(i)->tree_sibling_next) { (i) = ARC_HEAD((i)->match); GET_TREE_PARENT(i, i); }\ + if ((i)->is_tree_root) break;\ + (i) = (i)->tree_sibling_next;\ + }\ + } + +// i=parent, j=child +#define ADD_TREE_CHILD(i, j)\ + {\ + (j)->flag = 0;\ + (j)->tree = (i)->tree;\ + (j)->first_tree_child = NULL;\ + (j)->tree_sibling_next = (i)->first_tree_child;\ + if ((i)->first_tree_child)\ + {\ + (j)->tree_sibling_prev = (i)->first_tree_child->tree_sibling_prev;\ + (i)->first_tree_child->tree_sibling_prev = j;\ + }\ + else\ + {\ + (j)->tree_sibling_prev = j;\ + }\ + (i)->first_tree_child = j;\ + } + + +#define REMOVE_FROM_TREE(i)\ + {\ + if ((i)->tree_sibling_next) (i)->tree_sibling_next->tree_sibling_prev = (i)->tree_sibling_prev;\ + else\ + {\ + Node* i_NEXT = ARC_HEAD((i)->match); i_NEXT = ARC_HEAD(i_NEXT->tree_parent); i_NEXT = i_NEXT->first_tree_child;\ + i_NEXT->tree_sibling_prev = (i)->tree_sibling_prev;\ + }\ + if ((i)->tree_sibling_prev->tree_sibling_next) (i)->tree_sibling_prev->tree_sibling_next = (i)->tree_sibling_next;\ + else\ + {\ + Node* i_PARENT = ARC_HEAD((i)->match); i_PARENT = ARC_HEAD(i_PARENT->tree_parent);\ + i_PARENT->first_tree_child = (i)->tree_sibling_next;\ + }\ + } + + +#endif + diff --git a/blossom5-v2.05.src/PMinit.cpp b/blossom5-v2.05.src/PMinit.cpp new file mode 100644 index 0000000..cbcbb99 --- /dev/null +++ b/blossom5-v2.05.src/PMinit.cpp @@ -0,0 +1,544 @@ +#include +#include +#include +#include "PMimplementation.h" + + +void PerfectMatching::InitGreedy(bool allocate_trees) +{ + Node* i; + int dir; + Edge* a; + EdgeIterator I; + Tree* t = NULL; + Node* last_root = &nodes[node_num]; + REAL slack_min; + + for (i=nodes; iy = PM_INFTY; + for (a=edges; ahead[0]->y > a->slack) a->head[0]->y = a->slack; + if (a->head[1]->y > a->slack) a->head[1]->y = a->slack; + } + for (a=edges; ahead[0]; + if (!i->is_outer) + { + i->is_outer = 1; + i->y /= 2; + } + a->slack -= i->y; + i = a->head[1]; + if (!i->is_outer) + { + i->is_outer = 1; + i->y /= 2; + } + a->slack -= i->y; + } + + tree_num = node_num; + for (i=nodes; iflag == 2) continue; + slack_min = PM_INFTY; + FOR_ALL_EDGES(i, a, dir, I) if (slack_min > a->slack) slack_min = a->slack; + i->y += slack_min; + FOR_ALL_EDGES(i, a, dir, I) + { + if (a->slack <= slack_min && i->flag == 0 && a->head[dir]->flag == 0) + { + i->flag = 2; + a->head[dir]->flag = 2; + i->match = EDGE_DIR_TO_ARC(a, dir); + a->head[dir]->match = EDGE_DIR_TO_ARC(a, 1-dir); + tree_num -= 2; + } + a->slack -= slack_min; + } + } + if (allocate_trees) + { + if (tree_num > tree_num_max) + { + if (trees) free(trees); + tree_num_max = tree_num; + trees = (Tree*) malloc(tree_num_max*sizeof(Tree)); + } + t = trees; + } + for (i=nodes; iflag != 0) continue; + i->is_tree_root = 1; + i->first_tree_child = NULL; + i->tree_sibling_prev = last_root; + last_root->tree_sibling_next = i; + last_root = i; + if (allocate_trees) + { + i->tree = t; + t->root = i; + t->eps = 0; + t->first[0] = t->first[1] = NULL; + t->pq_current = NULL; + t->pq00.Reset(); + t->pq0.Reset(); + t->pq_blossoms.Reset(); + t ++; + } + } + last_root->tree_sibling_next = NULL; +} + +//////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////// + +PerfectMatching::Node* PerfectMatching::FindBlossomRootInit(Edge* a0) +{ + Node* i; + Node* j; + Node* _i[2]; + Node* r; + int branch; + + _i[0] = ARC_HEAD(a0); + _i[1] = ARC_TAIL(a0); + branch = 0; + while ( 1 ) + { + if (!_i[branch]->is_outer) + { + r = _i[branch]; + j = _i[1-branch]; + break; + } + _i[branch]->is_outer = 0; + if (_i[branch]->is_tree_root) + { + j = _i[branch]; + i = _i[1-branch]; + while (i->is_outer) + { + i->is_outer = 0; + i = ARC_HEAD(i->match); + i->is_outer = 0; + i = ARC_HEAD(i->tree_parent); + } + r = i; + break; + } + i = ARC_HEAD(_i[branch]->match); + i->is_outer = 0; + _i[branch] = ARC_HEAD(i->tree_parent); + branch = 1 - branch; + } + i = r; + while ( i != j ) + { + i = ARC_HEAD(i->match); + i->is_outer = 1; + i = ARC_HEAD(i->tree_parent); + i->is_outer = 1; + } + return r; +} + +void PerfectMatching::ShrinkInit(Edge* a0, Node* tree_root) +{ + int branch, flag; + Node* i; + Node* j; + Node* r; + Arc* a_prev; + Arc* aa; + + tree_root->flag = 2; + i = tree_root->first_tree_child; + if ( i ) + while ( 1 ) + { + ARC_HEAD(i->match)->flag = 2; + i->flag = 2; + + MOVE_NODE_IN_TREE(i); + } + + r = FindBlossomRootInit(a0); + + if ( !r->is_tree_root ) + { + j = ARC_HEAD(r->match); + j->match = aa = j->tree_parent; + i = ARC_HEAD(aa); + while ( !i->is_tree_root ) + { + j = ARC_HEAD(i->match); + i->match = ARC_REV(aa); + j->match = aa = j->tree_parent; + i = ARC_HEAD(aa); + } + i->match = ARC_REV(aa); + } + + tree_root->is_tree_root = 0; + + branch = 0; + flag = 0; + a_prev = EDGE_DIR_TO_ARC(a0, 0); + i = ARC_HEAD(a_prev); + while ( 1 ) + { + Arc* a_next = (flag == 0) ? i->match : i->tree_parent; + flag = 1 - flag; + i->flag = 0; + i->match = NULL; + if (branch == 0) + { + i->blossom_sibling = a_next; + if (i == r) + { + branch = 1; + flag = 0; + a_prev = ARC_REV(a0); + i = ARC_HEAD(a_prev); + if (i == r) break; + } + else + { + a_prev = i->blossom_sibling; + i = ARC_HEAD(a_prev); + } + } + else + { + i->blossom_sibling = ARC_REV(a_prev); + a_prev = a_next; + i = ARC_HEAD(a_prev); + if (i == r) break; + } + } + i->blossom_sibling = ARC_REV(a_prev); +} + +void PerfectMatching::ExpandInit(Node* k) +{ + Node* i = ARC_HEAD(k->blossom_sibling); + Node* j; + + while ( 1 ) + { + i->flag = 2; i->is_outer = 1; + if (i == k) break; + i->match = i->blossom_sibling; + j = ARC_HEAD(i->match); + j->flag = 2; j->is_outer = 1; + j->match = ARC_REV(i->match); + i = ARC_HEAD(j->blossom_sibling); + } +} + +void PerfectMatching::AugmentBranchInit(Node* i0, Node* r) +{ + Node* tree_root_prev = r->tree_sibling_prev; + Node* i; + Node* j; + Arc* aa; + + r->flag = 2; + i = r->first_tree_child; + if ( i ) + while ( 1 ) + { + ARC_HEAD(i->match)->flag = 2; + i->flag = 2; + + MOVE_NODE_IN_TREE(i); + } + i = i0; + if ( !i0->is_tree_root ) + { + j = ARC_HEAD(i0->match); + j->match = aa = j->tree_parent; + i = ARC_HEAD(aa); + while ( !i->is_tree_root ) + { + j = ARC_HEAD(i->match); + i->match = ARC_REV(aa); + j->match = aa = j->tree_parent; + i = ARC_HEAD(aa); + } + i->match = ARC_REV(aa); + } + r->is_tree_root = 0; + tree_root_prev->tree_sibling_next = r->tree_sibling_next; + if (r->tree_sibling_next) r->tree_sibling_next->tree_sibling_prev = tree_root_prev; + tree_num --; +} + +//////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////// + +// true_slack(a) = slack(a) + ... + +// i->flag=0, i->is_processed=1: true_slack -= eps +// i->flag=1, i->match->head->is_processed=1: true_slack += eps - slack(i->match) + +void PerfectMatching::InitGlobal() +{ + Node* i; + Node* j; + Node* r; + Node* r2; + Node* r3 = NULL; // initialize to prevent compiler warning + Edge* a; + EdgeIterator I; + int dir; + Tree TREE; + enum { NONE, AUGMENT, SHRINK } flag; + + InitGreedy(); + + for (i=nodes; ibest_edge = NULL; + + PriorityQueue pq; + + for (r=nodes[node_num].tree_sibling_next; r; ) + { + r2 = r->tree_sibling_next; + if (r2) r3 = r2->tree_sibling_next; + i = r; + + pq.Reset(); + + r->tree = &TREE; + + REAL eps = 0; + Arc* critical_arc = NULL; + REAL critical_eps = PM_INFTY; + flag = NONE; + Node* branch_root = i; + + while ( 1 ) + { + i->is_processed = 1; + i->y -= eps; + if (!i->is_tree_root) ARC_HEAD(i->match)->y += eps; + + FOR_ALL_EDGES(i, a, dir, I) + { + a->slack += eps; + j = a->head[dir]; + + if (j->tree == &TREE) + { + // same tree + if (j->flag == 0) + { + REAL slack = a->slack; + if (!j->is_processed) slack += eps; + if (2*critical_eps > slack || critical_arc == NULL) + { + flag = SHRINK; + critical_eps = slack/2; + critical_arc = EDGE_DIR_TO_ARC(a, dir); + if (critical_eps <= eps) break; + //pq.DecreaseUpperBound(critical_eps); + } + } + } + else if (j->flag == 0) + { + // different tree + if (critical_eps >= a->slack || critical_arc == NULL) + { + flag = AUGMENT; + critical_eps = a->slack; + critical_arc = EDGE_DIR_TO_ARC(a, dir); + if (critical_eps <= eps) break; + //pq.DecreaseUpperBound(critical_eps); + } + } + else + { + // free node + if (a->slack > eps) + { + if (a->slack < critical_eps) + { + if (j->best_edge == NULL) + { + j->best_edge = a; + pq.Add(a); + } + else + { + if (a->slack < j->best_edge->slack) + { + pq.Decrease(j->best_edge, a, pq_buf); + j->best_edge = a; + } + } + } + } + else + { + assert(j->flag == 2 && !j->is_blossom && !ARC_HEAD(j->match)->is_blossom); + if (j->best_edge) pq.Remove(j->best_edge, pq_buf); + j->flag = 1; + j->tree = i->tree; + j->tree_parent = EDGE_DIR_TO_ARC(a, 1-dir); + j = ARC_HEAD(j->match); + if (j->best_edge) pq.Remove(j->best_edge, pq_buf); + ADD_TREE_CHILD(i, j); + } + } + } + + if (dir < 2 && a) + { + Edge* atmp = a; + int dirtmp = dir; + CONTINUE_FOR_ALL_EDGES(i, atmp, dirtmp, I) atmp->slack += eps; + break; + } + + // move i + if (i->first_tree_child) i = i->first_tree_child; + else + { + while (i != branch_root && !i->tree_sibling_next) { i = ARC_HEAD(i->match); i = ARC_HEAD(i->tree_parent); } + if (i == branch_root) + { + PriorityQueue::Item* q = pq.GetMin(); + if (q == NULL || q->slack >= critical_eps) + { + eps = critical_eps; + break; + } + pq.Remove(q, pq_buf); + a = (Edge*)q; + dir = (a->head[0]->flag == 2) ? 0 : 1; + j = a->head[0]; + Arc* aa = EDGE_DIR_TO_ARC(a, dir); + eps = a->slack; + assert(eps < critical_eps); + + // continue growth + i = ARC_TAIL(aa); + j = ARC_HEAD(aa); + + assert(j->flag == 2 && !j->is_blossom && !ARC_HEAD(j->match)->is_blossom); + j->flag = 1; + j->tree = i->tree; + j->tree_parent = ARC_REV(aa); + j = ARC_HEAD(j->match); + if (j->best_edge) pq.Remove(j->best_edge, pq_buf); + ADD_TREE_CHILD(i, j); + i = branch_root = j; + continue; + } + i = i->tree_sibling_next; + } + } + + // update slacks + i = r; + while ( 1 ) + { + if (i->is_processed) + { + i->y += eps; + if (!i->is_tree_root) + { + j = ARC_HEAD(i->match); + j->y -= eps; + REAL delta = eps - ARC_TO_EDGE_PTR(i->match)->slack; + FOR_ALL_EDGES(j, a, dir, I) a->slack += delta; + j->best_edge = NULL; + } + FOR_ALL_EDGES(i, a, dir, I) + { + if (!PriorityQueue::isReset(a)) + { + assert(a->head[dir]->flag == 2 && a->head[dir]->best_edge == a); + a->head[dir]->best_edge = NULL; + PriorityQueue::ResetItem(a); + } + a->slack -= eps; + } + + i->is_processed = 0; + } + else + { + if (!i->is_tree_root) ARC_HEAD(i->match)->best_edge = NULL; + } + i->best_edge = NULL; + + MOVE_NODE_IN_TREE(i); + } + + i = ARC_TAIL(critical_arc); + j = ARC_HEAD(critical_arc); + if (flag == SHRINK) + { + // shrink + ShrinkInit(ARC_TO_EDGE_PTR(critical_arc), r); + } + else + { + // augment + AugmentBranchInit(i, r); + if (j->is_outer) + { + AugmentBranchInit(j, j); + } + else + { + ExpandInit(j); + tree_num --; + } + i->match = critical_arc; + j->match = ARC_REV(critical_arc); + } + + r = r2; + if (r && !r->is_tree_root) r = r3; + } + + if (tree_num > tree_num_max) + { + if (trees) free(trees); + tree_num_max = tree_num; + trees = (Tree*) malloc(tree_num_max*sizeof(Tree)); + } + Tree* t = trees; + for (r=nodes; ris_outer) + { + ExpandInit(r); + r->is_tree_root = 1; + r->flag = 0; + r->first_tree_child = NULL; + if (t == trees) { nodes[node_num].tree_sibling_next = r; r->tree_sibling_prev = &nodes[node_num]; } + else { (t-1)->root->tree_sibling_next = r; r->tree_sibling_prev = (t-1)->root; } + r->tree = t; + t->root = r; + t->eps = 0; + t->first[0] = t->first[1] = NULL; + t->pq_current = NULL; + t->pq00.Reset(); + t->pq0.Reset(); + t->pq_blossoms.Reset(); + t ++; + } + } + assert(t == trees+tree_num); + if (t == trees) nodes[node_num].tree_sibling_next = NULL; + else (t-1)->root->tree_sibling_next = NULL; +} diff --git a/blossom5-v2.05.src/PMinterface.cpp b/blossom5-v2.05.src/PMinterface.cpp new file mode 100644 index 0000000..baab2d6 --- /dev/null +++ b/blossom5-v2.05.src/PMinterface.cpp @@ -0,0 +1,270 @@ +#include +#include +#include +#include "PMimplementation.h" + + +PerfectMatching::PerfectMatching(int nodeNum, int edgeNumMax) + : node_num(nodeNum), + edge_num(0), + edge_num_max(edgeNumMax), + trees(NULL), + tree_num_max(0), + removed_first(NULL), + blossom_num(0), + removed_num(0), + first_solve(true) +{ + if (node_num & 1) { printf("# of nodes is odd: perfect matching cannot exist\n"); exit(1); } + nodes = (Node*) malloc((node_num+1)*sizeof(Node)); + edges_orig = (char*) malloc(edge_num_max*sizeof(Edge)+1); + edges = (Edge*) ( ( ((POINTER_TYPE)edges_orig) & 1 ) ? (edges_orig + 1) : edges_orig ); + memset(nodes, 0, (node_num+1)*sizeof(Node)); + + blossoms = new DBlock(256); + tree_edges = new DBlock(256); + expand_tmp_list = new Block(256); + pq_buf = PriorityQueue::AllocateBuf(); +} + + +void PerfectMatching::Save(char* filename, int format) +{ + if (!first_solve) { printf("Save() cannot be called after Solve()!\n"); exit(1); } + int e; + FILE* fp = fopen(filename, "w"); + if (!fp) { printf("Can't open %s\n", filename); exit(1); } + if (format == 0) + { + fprintf(fp, "p edge %d %d\n", node_num, edge_num); + for (e=0; e::DeallocateBuf(pq_buf); +} + + +PerfectMatching::EdgeId PerfectMatching::AddEdge(NodeId _i, NodeId _j, REAL cost) +{ + if (_i<0 || _i>=node_num || _j<0 || _j>node_num || _i==_j) + { + printf("wrong node id's! (%d,%d)\n", _i, _j); exit(1); + } + if (edge_num >= edge_num_max) ReallocateEdges(); + Node* i = nodes + _i; + Node* j = nodes + _j; + Edge* a = edges + edge_num; + + ADD_EDGE(i, a, 0); + ADD_EDGE(j, a, 1); + a->head0[0] = j; + a->head0[1] = i; + + a->slack = cost*COST_FACTOR; + PriorityQueue::ResetItem(a); + + return edge_num ++; +} + +int PerfectMatching::GetSolution(EdgeId e) +{ + assert(e>=0 && ehead0[1]->match == EDGE_DIR_TO_ARC(a, 0)) ? 1 : 0; +} + +PerfectMatching::NodeId PerfectMatching::GetMatch(NodeId i) +{ + assert(i>=0 && ihead0[1]; !i->is_outer; i=i->blossom_parent, delta--) {} + for (j=a->head0[0]; !j->is_outer; j=j->blossom_parent, delta++) {} + if ( i == j ) + { + i = a->head0[1]; + j = a->head0[0]; + while ( delta < 0 ) { i = i->blossom_parent; delta ++; } + while ( delta > 0 ) { j = j->blossom_parent; delta --; } + while ( i->blossom_parent != j->blossom_parent ) + { + i = i->blossom_parent; + j = j->blossom_parent; + } + } + tail = i; + head = j; + assert((i->is_outer && j->is_outer) || (i->blossom_parent==j->blossom_parent && !i->is_outer && !j->is_outer)); +} + +void PerfectMatching::ReallocateEdges() +{ + printf("Warning: reallocating edges. Increasing edge_num_max in the constructor may improve memory efficiency!\n"); + edge_num_max = edge_num_max*3/2 + 16; + char* edges_orig_old = edges_orig; + Edge* edges_old = edges; + edges_orig = (char*) realloc(edges_orig_old, edge_num_max*sizeof(Edge)+1); + edges = (Edge*) ( ( ((POINTER_TYPE)edges_orig_old) & 1 ) ? (edges_orig + 1) : edges_orig ); + if ( ((POINTER_TYPE)edges) & 1 ) + { + char* edges_orig_old2 = edges_orig; + Edge* edges_old2 = edges; + + edges_orig = (char*) malloc(edge_num_max*sizeof(Edge)+1); + edges = (Edge*) ( ( ((POINTER_TYPE)edges_orig_old) & 1 ) ? (edges_orig + 1) : edges_orig ); + memcpy(edges, edges_old2, edge_num*sizeof(Edge)); + free(edges_orig_old2); + } + +#define UPDATE_EDGE_PTR(ptr) ptr = (Edge*)((char*)(ptr) + ((char*)edges - (char*)edges_old)) +#define UPDATE_ARC_PTR(ptr) ptr = (Arc*)((char*)(ptr) + ((char*)edges - (char*)edges_old)) + + Node* i; + Edge* a; + for (a=edges; anext[0]) UPDATE_EDGE_PTR(a->next[0]); + if (a->next[1]) UPDATE_EDGE_PTR(a->next[1]); + if (a->prev[0]) UPDATE_EDGE_PTR(a->prev[0]); + if (a->prev[1]) UPDATE_EDGE_PTR(a->prev[1]); + } + if (first_solve) + { + for (i=nodes; ifirst[0]) UPDATE_EDGE_PTR(i->first[0]); + if (i->first[1]) UPDATE_EDGE_PTR(i->first[1]); + } + } + else + { + Node* i0; + for (i0=nodes; i0is_outer) + { + UPDATE_ARC_PTR(i->match); + if (i->first[0]) UPDATE_EDGE_PTR(i->first[0]); + if (i->first[1]) UPDATE_EDGE_PTR(i->first[1]); + break; + } + UPDATE_ARC_PTR(i->blossom_sibling); + if (i->first[0]) UPDATE_EDGE_PTR(i->first[0]); + if (i->first[1]) UPDATE_EDGE_PTR(i->first[1]); + + i = i->blossom_parent; + if (i->is_outer) { if ( i->is_marked) break; i->is_marked = 1; } + else { if (!i->is_marked) break; i->is_marked = 0; } + } + } + for (i0=nodes; i0is_outer) break; + + i = i->blossom_parent; + if (i->is_outer) { if (!i->is_marked) break; i->is_marked = 0; } + else { if ( i->is_marked) break; i->is_marked = 1; } + } + } + } +} + +int PerfectMatching::GetBlossomNum() +{ + return blossom_num; +} + +void PerfectMatching::GetDualSolution(int* blossom_parents, REAL* twice_y) +{ + int _i0, id = node_num; + int* child_ptr; + Node* i0; + Node* i; + int* tmp_array = new int[blossom_num]; + + int* tmp_array_ptr = tmp_array; + for (_i0=0, i0=nodes; _i0y; + if (i0->is_outer) + { + blossom_parents[_i0] = -1; + continue; + } + child_ptr = &blossom_parents[_i0]; + i = i0->blossom_parent; + while ( 1 ) + { + if (i->is_marked) + { + *child_ptr = i->lca_preorder; + break; + } + i->is_marked = 1; + *tmp_array_ptr ++ = i->lca_preorder; + *child_ptr = i->lca_preorder = id ++; + child_ptr = &blossom_parents[i->lca_preorder]; + twice_y[i->lca_preorder] = i->y; + if (i->is_outer) + { + *child_ptr = -1; + break; + } + i = i->blossom_parent; + } + } + + assert(id == node_num+blossom_num && tmp_array_ptr == tmp_array + blossom_num); + + tmp_array_ptr = tmp_array; + for (_i0=0, i0=nodes; _i0is_outer) continue; + i = i0->blossom_parent; + while ( 1 ) + { + if (!i->is_marked) break; + i->is_marked = 0; + i->lca_preorder = *tmp_array_ptr ++; + if (i->is_outer) break; + i = i->blossom_parent; + } + } + + delete [] tmp_array; +} diff --git a/blossom5-v2.05.src/PMmain.cpp b/blossom5-v2.05.src/PMmain.cpp new file mode 100644 index 0000000..0914fc6 --- /dev/null +++ b/blossom5-v2.05.src/PMmain.cpp @@ -0,0 +1,694 @@ +#include +#include +#include +#include "PMimplementation.h" + + +void PerfectMatching::Finish() +{ + +#define IS_VALID_MATCH(i) ((Edge*)(i->match) >= edges && (Edge*)(i->match) < edges + edge_num) + + Node* i0; + Node* i; + Node* j; + Node* k; + Node* b; + Node* b_prev; + Node* b_prev_prev; + + for (i0=nodes; i0blossom_grandparent = b_prev; + b_prev = b; + b = b->blossom_parent; + } while (!IS_VALID_MATCH(b)); + + b_prev_prev = b_prev->blossom_grandparent; + while ( 1 ) + { + for (k=ARC_TAIL0(b->match); k->blossom_parent!=b; k=k->blossom_parent) {} + k->match = b->match; + i = ARC_HEAD(k->blossom_sibling); + while ( i != k ) + { + i->match = i->blossom_sibling; + j = ARC_HEAD(i->match); + j->match = ARC_REV(i->match); + i = ARC_HEAD(j->blossom_sibling); + } + + b = b_prev; + if (!b->is_blossom) break; + b_prev = b_prev_prev; + b_prev_prev = b_prev->blossom_grandparent; + } + } +} + + +void PerfectMatching::AddTreeEdge(Tree* t0, Tree* t1) +{ + TreeEdge* e = tree_edges->New(); + e->head[0] = t1; + e->head[1] = t0; + e->next[0] = t0->first[0]; + t0->first[0] = e; + e->next[1] = t1->first[1]; + t1->first[1] = e; + + e->pq00.Reset(); + e->pq01[0].Reset(); + e->pq01[1].Reset(); + + t1->pq_current = e; + t1->dir_current = 0; +} + +bool PerfectMatching::ProcessEdge00(Edge* a, bool update_boundary_edge) +{ + int dir; + Node* j; + Node* prev[2]; + Node* last[2]; + for (dir=0; dir<2; dir++) + { + if (a->head[dir]->is_outer) { prev[dir] = NULL; last[dir] = a->head[dir]; } + else + { + j = a->head[dir]; + GET_PENULTIMATE_BLOSSOM(j); + prev[dir] = j; + last[dir] = prev[dir]->blossom_parent; + //assert(last[dir]->is_outer); + } + } + + if (last[0] != last[1]) + { + for (dir=0; dir<2; dir++) + { + j = a->head[dir]; + if (j != last[dir]) { int dir_rev = 1 - dir; MOVE_EDGE(j, last[dir], a, dir_rev); } + } + if (update_boundary_edge) a->slack -= 2*a->head[0]->tree->eps; + return true; + } + + if (prev[0] != prev[1]) + { + for (dir=0; dir<2; dir++) + { + j = a->head[dir]; + if (j != prev[dir]) { int dir_rev = 1 - dir; MOVE_EDGE(j, prev[dir], a, dir_rev); } + } + a->slack -= 2*prev[0]->blossom_eps; + return false; + } + + for (dir=0; dir<2; dir++) + { + j = a->head[1-dir]; + REMOVE_EDGE(j, a, dir); + } + a->next[0] = prev[0]->blossom_selfloops; + prev[0]->blossom_selfloops = a; + return false; +} + + +inline void PerfectMatching::AugmentBranch(Node* i0) +{ + int dir; + Tree* t = i0->tree; + Node* r = t->root; + Node* tree_root_prev = r->tree_sibling_prev; + Node* i; + Node* j; + Edge* a; + EdgeIterator I; + Arc* aa; + REAL eps = t->eps; + PriorityQueue::Item* q; + TreeEdge* e; + TreeEdgeIterator T; + Tree* t2; + + t = r->tree; + t->pq_current = t; + + FOR_ALL_TREE_EDGES_X(t, e, dir, T) + { + t2 = e->head[dir]; + e->head[1-dir] = NULL; // mark it for deletion + + t2->pq_current = e; + t2->dir_current = dir; + } + + i = r->first_tree_child; + if ( i ) + while ( 1 ) + { + Node* i0 = i; + i = ARC_HEAD(i->match); + if (i->is_processed) + { + if (i->is_blossom) + { + a = ARC_TO_EDGE_PTR(i->match); + REAL tmp = a->slack; a->slack = i->y; i->y = tmp; + PriorityQueue::ResetItem(a); + } + FOR_ALL_EDGES(i, a, dir, I) + { + GET_OUTER_HEAD(a, dir, j); + + if (j->flag == 0 && j->is_processed) + { + if (j->tree != t) + { + a->slack += eps; + if (PriorityQueue::isReset(a)) j->tree->pq0.Add(a); + } + } + else a->slack += eps; + } + } + + i = i0; + MOVE_NODE_IN_TREE(i); + } + + /////////////////////////////////////////////////////////////////// + + FOR_ALL_TREE_EDGES(t, e, dir) + { + t2 = e->head[dir]; + t2->pq_current = NULL; + + e->pq01[1-dir].Merge(t2->pq0); + for (q=e->pq00.GetFirst(); q; q=e->pq00.GetNext(q)) + { + q->slack -= eps; + int dir2; + for (dir2=0; dir2<2; dir2++) GET_OUTER_HEAD((Edge*)q, dir2, j); + } + e->pq00.Merge(t2->pq0); + for (q=e->pq01[dir].GetAndResetFirst(); q; q=e->pq01[dir].GetAndResetNext()) + { + q->slack -= eps; + int dir2; + for (dir2=0; dir2<2; dir2++) GET_OUTER_HEAD((Edge*)q, dir2, j); + } + } + for (q=t->pq0.GetAndResetFirst(); q; q=t->pq0.GetAndResetNext()) + { + q->slack -= eps; + int dir2; + for (dir2=0; dir2<2; dir2++) GET_OUTER_HEAD((Edge*)q, dir2, j); + } + for (q=t->pq00.GetAndResetFirst(); q; q=t->pq00.GetAndResetNext()) + { + ProcessEdge00((Edge*)q); + } + + /////////////////////////////////////////////////////////////////// + + r->flag = 2; + r->is_processed = 0; + i = r->first_tree_child; + r->y += eps; + if ( i ) + while ( 1 ) + { + j = ARC_HEAD(i->match); + j->flag = 2; + i->flag = 2; + j->is_processed = 0; + i->is_processed = 0; + j->y -= eps; + i->y += eps; + + MOVE_NODE_IN_TREE(i); + } + + /////////////////////////////////////////////////////////////////// + + i = i0; + if ( !i0->is_tree_root ) + { + j = ARC_HEAD(i0->match); + GET_TREE_PARENT(j, i); + j->match = aa = j->tree_parent; + while ( !i->is_tree_root ) + { + j = ARC_HEAD(i->match); + i->match = ARC_REV(aa); + GET_TREE_PARENT(j, i); + j->match = aa = j->tree_parent; + } + i->match = ARC_REV(aa); + } + r->is_tree_root = 0; + tree_root_prev->tree_sibling_next = r->tree_sibling_next; + if (r->tree_sibling_next) r->tree_sibling_next->tree_sibling_prev = tree_root_prev; + tree_num --; +} + + +void PerfectMatching::Augment(Edge* a) +{ + Node* j; + int dir; + + for (dir=0; dir<2; dir++) + { + GET_OUTER_HEAD(a, dir, j); + AugmentBranch(j); + j->match = EDGE_DIR_TO_ARC(a, 1-dir); + } + if (options.verbose) + { + int k = 1; + while (k < tree_num) k *= 2; + if (k == tree_num || tree_num<=8 || (tree_num<=64 && (tree_num%8)==0)) { printf("%d.", tree_num); fflush(stdout); } + } +} + +inline void PerfectMatching::GrowNode(Node* i) +{ + //assert(i->is_outer); + //assert(i->flag == 0); + + Edge* a; + EdgeIterator I; + int dir; + Node* j; + Tree* t = i->tree; + REAL eps = t->eps; + Edge* a_augment = NULL; + + FOR_ALL_EDGES(i, a, dir, I) + { + GET_OUTER_HEAD(a, dir, j); + + if (j->flag == 2) + { + a->slack += eps; + if (a->slack > 0) + { + t->pq0.Add(a); + } + else + { + j->flag = 1; + j->tree = i->tree; + j->tree_parent = EDGE_DIR_TO_ARC(a, 1-dir); + j->y += eps; + j = ARC_HEAD(j->match); + j->y -= eps; + ADD_TREE_CHILD(i, j); + } + } + else + { + if (j->flag == 0 && j->is_processed) + { + if (!PriorityQueue::isReset(a)) j->tree->pq0.Remove(a, pq_buf); + if (a->slack <= j->tree->eps && j->tree != t) a_augment = a; + a->slack += eps; + if (!j->tree->pq_current) AddTreeEdge(t, j->tree); + j->tree->pq_current->pq00.Add(a); + } + else + { + a->slack += eps; + if (j->flag == 1 && j->tree != t) + { + if (!j->tree->pq_current) AddTreeEdge(t, j->tree); + j->tree->pq_current->pq01[j->tree->dir_current].Add(a); + } + } + } + } + + //assert(!i->is_processed); + i->is_processed = 1; + + if (!i->is_tree_root) + { + j = ARC_HEAD(i->match); + //assert(!j->is_processed); + j->is_processed = 1; + if (j->is_blossom) + { + a = ARC_TO_EDGE_PTR(i->match); + REAL tmp = a->slack; a->slack = j->y; j->y = tmp; + t->pq_blossoms.Add(a); + } + } + + if (a_augment) Augment(a_augment); + + stat.grow_count ++; +} + + + +void PerfectMatching::GrowTree(Node* r, bool new_subtree) +{ + //assert(r->flag == 0); + + Node* i = r; + Node* j; + Node* stop = r->tree_sibling_next; + if (new_subtree && r->first_tree_child) stop = r->first_tree_child; + Edge* a; + EdgeIterator I; + int dir; + Tree* t = r->tree; + REAL eps = t->eps; + int tree_num0 = tree_num; + + while ( 1 ) + { + if (!i->is_tree_root) + { + // process "-" node + i = ARC_HEAD(i->match); + FOR_ALL_EDGES(i, a, dir, I) + { + GET_OUTER_HEAD(a, dir, j); + + if (j->flag == 2) a->slack -= eps; + else + { + if (j->flag == 0 && j->is_processed) + { + if (!PriorityQueue::isReset(a)) j->tree->pq0.Remove(a, pq_buf); + a->slack -= eps; + if (j->tree != t) + { + if (!j->tree->pq_current) AddTreeEdge(t, j->tree); + j->tree->pq_current->pq01[1-j->tree->dir_current].Add(a); + } + } + else a->slack -= eps; + } + } + i = ARC_HEAD(i->match); + } + // process "+" node + GrowNode(i); + if (tree_num != tree_num0) break; + + if (i->first_tree_child) i = i->first_tree_child; + else + { + while (i != r && !i->tree_sibling_next) { i = ARC_HEAD(i->match); GET_TREE_PARENT(i, i); } + i = i->tree_sibling_next; + } + if (i == stop) break; + } +} + +void PerfectMatching::Solve(bool finish) +{ + Node* i; + Node* j; + Node* r; + Node* r2; + Node* r3 = NULL; // initialize to prevent compiler warning + PriorityQueue::Item* q; + Edge* a; + Tree* t; + Tree* t2; + TreeEdge* e; + TreeEdgeIterator T; + int dir; + REAL eps; + + double start_time = get_time(); + + if (IS_INT) + { + if (options.dual_greedy_update_option == 2) + { + printf("Fixed eps approach can only be used with floating point REAL!\n"); + printf("Change REAL to double in PerfectMatching.h and recompile\n"); + exit(1); + } + if (options.dual_LP_threshold > 0) + { + printf("LP approach can only be used with floating point REAL!\n"); + printf("Change REAL to double in PerfectMatching.h and recompile\n"); + exit(1); + } + } + if (options.verbose) { printf("perfect matching with %d nodes and %d edges\n", node_num, edge_num); fflush(stdout); } + + if (first_solve) + { + if (options.verbose) { printf(" starting init..."); fflush(stdout); } + if (options.fractional_jumpstart) InitGlobal(); + else InitGreedy(); + if (options.verbose) printf("done [%.3f secs]. ", get_time() - start_time); + first_solve = false; + } + else if (options.verbose) printf(" solving updated problem. "); + + if (options.verbose) { printf("%d trees\n .", tree_num); fflush(stdout); } + + memset(&stat, 0, sizeof(Stat)); + + /////////////////////////////////////////////////////// + // first pass - initialize auxiliary graph // + /////////////////////////////////////////////////////// + + for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) + { + //assert(!r->is_processed); + t = r->tree; + //assert(!t->first[0] && !t->first[1]); + + EdgeIterator I; + FOR_ALL_EDGES(r, a, dir, I) + { + j = a->head[dir]; + if (j->flag == 2) t->pq0.Add(a); + else if (j->is_processed) + { + //assert(j->flag == 0); + if (!j->tree->pq_current) AddTreeEdge(t, j->tree); + j->tree->pq_current->pq00.Add(a); + } + } + r->is_processed = 1; + FOR_ALL_TREE_EDGES(t, e, dir) e->head[dir]->pq_current = NULL; + } + + /////////////////////////////////////////////////////// + // main loop // + /////////////////////////////////////////////////////// + + while ( 1 ) + { + int tree_num0 = tree_num; + Stat stat0 = stat; + REAL delta = 0; + + for (r=nodes[node_num].tree_sibling_next; r; ) + { + r2 = r->tree_sibling_next; + if (r2) r3 = r2->tree_sibling_next; + t = r->tree; + + int tree_num1 = tree_num; + + ////////////////////////////////////////////////////////////////////// + // step 1 - traversing auxiliary graph, setting pq_current pointers // + ////////////////////////////////////////////////////////////////////// + t->pq_current = t; + if (options.update_duals_before) + { + eps = PM_INFTY; + Edge* a_augment = NULL; + REAL eps_augment = PM_INFTY; + if ((q=t->pq0.GetMin())) eps = q->slack; + if ((q=t->pq_blossoms.GetMin()) && eps > q->slack) eps = q->slack; + while ((q=t->pq00.GetMin())) + { + if (ProcessEdge00((Edge*)q, false)) break; + t->pq00.Remove(q, pq_buf); + } + if (q && 2*eps > q->slack) eps = q->slack/2; + FOR_ALL_TREE_EDGES_X(t, e, dir, T) + { + t2 = e->head[dir]; + t2->pq_current = e; + t2->dir_current = dir; + if ((q=e->pq00.GetMin()) && (!a_augment || eps_augment > q->slack-t2->eps)) { a_augment = (Edge*)q; eps_augment = q->slack-t2->eps; } + if ((q=e->pq01[dir].GetMin()) && eps > q->slack+t2->eps) eps = q->slack+t2->eps; + } + if (eps > eps_augment) eps = eps_augment; + if (eps > t->eps) + { + delta += eps - t->eps; + t->eps = eps; + } + if (a_augment && eps_augment <= t->eps) Augment(a_augment); + } + else + { + FOR_ALL_TREE_EDGES_X(t, e, dir, T) + { + t2 = e->head[dir]; + t2->pq_current = e; + t2->dir_current = dir; + + if ((q=e->pq00.GetMin()) && (q->slack - t->eps <= t2->eps)) + { + Augment((Edge*)q); + break; + } + } + } + + ///////////////////////////////// + // step 2 - growing tree // + ///////////////////////////////// + eps = t->eps; + REAL twice_eps = 2*eps; + + while ( tree_num1 == tree_num ) + { + if ((q=t->pq0.GetMin()) && q->slack <= t->eps) + { + a = (Edge*)q; + dir = (a->head[1]->flag == 2 && a->head[1]->is_outer) ? 1 : 0; + GET_OUTER_HEAD(a, 1-dir, i); + j = a->head[dir]; + //assert(i->flag==0 && j->flag==2 && i->is_outer && j->is_outer && i->tree==t); + + j->flag = 1; + j->tree = i->tree; + j->tree_parent = EDGE_DIR_TO_ARC(a, 1-dir); + j->y += eps; + j = ARC_HEAD(j->match); + j->y -= eps; + ADD_TREE_CHILD(i, j); + + GrowTree(j, true); + } + else if ((q=t->pq00.GetMin()) && q->slack <= twice_eps) + { + t->pq00.Remove(q, pq_buf); + a = (Edge*)q; + if (ProcessEdge00(a)) Shrink(a); + } + else if ((q=t->pq_blossoms.GetMin()) && q->slack <= eps) + { + t->pq_blossoms.Remove(q, pq_buf); + a = (Edge*)q; + j = (a->head[0]->flag == 1) ? a->head[0] : a->head[1]; + REAL tmp = a->slack; a->slack = j->y; j->y = tmp; + Expand(j); + } + else break; + } + + /////////////////////////////////////////////////////////////////////// + // step 3 - traversing auxiliary graph, clearing pq_current pointers // + /////////////////////////////////////////////////////////////////////// + if ( tree_num1 == tree_num ) + { + t->pq_current = NULL; + if (options.update_duals_after) + { + eps = PM_INFTY; + Edge* a_augment = NULL; + REAL eps_augment = PM_INFTY; + if ((q=t->pq0.GetMin())) eps = q->slack; + if ((q=t->pq_blossoms.GetMin()) && eps > q->slack) eps = q->slack; + while ((q=t->pq00.GetMin())) + { + if (ProcessEdge00((Edge*)q, false)) break; + t->pq00.Remove(q, pq_buf); + } + if (q && 2*eps > q->slack) eps = q->slack/2; + FOR_ALL_TREE_EDGES(t, e, dir) + { + t2 = e->head[dir]; + e->head[dir]->pq_current = NULL; + if ((q=e->pq00.GetMin()) && (!a_augment || eps_augment > q->slack-t2->eps)) { a_augment = (Edge*)q; eps_augment = q->slack-t2->eps; } + if ((q=e->pq01[dir].GetMin()) && eps > q->slack+t2->eps) eps = q->slack+t2->eps; + } + if (eps > eps_augment) eps = eps_augment; + bool progress = false; + if (eps > t->eps) + { + delta += eps - t->eps; + t->eps = eps; + progress = true; + } + if (a_augment && eps_augment <= t->eps) Augment(a_augment); + else if (progress && tree_num >= options.single_tree_threshold*node_num) + { + // continue with the same tree + r = t->root; + continue; + } + } + else + { + FOR_ALL_TREE_EDGES(t, e, dir) e->head[dir]->pq_current = NULL; + } + } + + /////////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////////////// + + r = r2; + if (r && !r->is_tree_root) r = r3; + } + + if (tree_num == 0) break; + + if ( tree_num == tree_num0 ) + //&& stat.grow_count == stat0.grow_count + //&& stat.shrink_count == stat0.shrink_count + //&& stat.expand_count == stat0.expand_count ) + { + if (!UpdateDuals()) + { + if (!IS_INT && delta <= PM_THRESHOLD) // for numerical stability + { + //CommitEps(); + int dual_greedy_update_option = options.dual_greedy_update_option; + options.dual_greedy_update_option = 2; + UpdateDuals(); + options.dual_greedy_update_option = dual_greedy_update_option; + } + } + } + } + + if (finish) Finish(); + + if (options.verbose) + { + printf("\ndone [%.3f secs]. %d grows, %d expands, %d shrinks\n", get_time()-start_time, stat.grow_count, stat.expand_count, stat.shrink_count); + printf(" expands: [%.3f secs], shrinks: [%.3f secs], dual updates: [%.3f secs]\n", stat.expand_time, stat.shrink_time, stat.dual_time); + fflush(stdout); + } +} + + + diff --git a/blossom5-v2.05.src/PMrepair.cpp b/blossom5-v2.05.src/PMrepair.cpp new file mode 100644 index 0000000..978b38f --- /dev/null +++ b/blossom5-v2.05.src/PMrepair.cpp @@ -0,0 +1,431 @@ +#include +#include +#include +#include "PMimplementation.h" + + +struct PerfectMatching::LCATreeX : LCATree +{ + LCATreeX(int size) : LCATree(size) { rev_mapping = new Node*[size]; } + ~LCATreeX() { delete [] rev_mapping; } + Node** rev_mapping; +}; + +void PerfectMatching::StartUpdate() +{ + Node* i0; + Node* i; + Node* j; + Node* b; + + while ((i=removed_first)) + { + removed_first = i->tree_sibling_next; + blossoms->Delete(i); + removed_num --; + } + + Edge* a; + Edge* selfloop_first = NULL; + Edge* selfloop_last = NULL; + + for (i0=nodes; i0is_processed = 0; + if (i0->is_outer) continue; + + i0->is_tree_root = 0; + i0->blossom_ptr = NULL; + i = i0; + while ( 1 ) + { + j = i->blossom_parent; + j->is_processed = 0; + if (j->is_outer) { j->first_tree_child = i; break; } + if (j->is_marked) break; + if ((a=j->blossom_selfloops)) + { + if (selfloop_last) selfloop_last->next[1] = a; + else selfloop_first = a; + selfloop_last = a; + a->next[1] = NULL; + } + j->blossom_ptr = i; + i = j; + } + b = (i->blossom_parent->is_outer) ? i->blossom_parent : i->blossom_parent->blossom_grandparent; +#ifdef LCA_REPAIRS + if (!b->is_marked) + { + b->lca_size = 1; + b->is_marked = 1; + } +#endif + while ( 1 ) + { +#ifdef LCA_REPAIRS + b->lca_size ++; +#endif + ARC_TO_EDGE_PTR(i->blossom_sibling)->y_saved = i->y; + i->y += i->blossom_parent->y; + if (!i->is_blossom) break; + i->is_marked = 1; + j = i; + i = i->blossom_ptr; + j->blossom_grandparent = b; + } + i->blossom_grandparent = b; + } + +#ifdef LCA_REPAIRS + for (i0=nodes; i0is_outer) continue; + b = i0->blossom_grandparent; + if (!b->is_marked) continue; + b->is_marked = 0; + LCATreeX* lca = new LCATreeX(b->lca_size); + b->blossom_ptr = b->first_tree_child; + i = b; + while ( 1 ) + { + if (i->blossom_ptr) i = i->blossom_ptr; + else + { + while ( 1 ) + { + if (i->is_outer) break; + i->lca_preorder = lca->Add(i, i->blossom_parent); + lca->rev_mapping[i->lca_preorder] = i; + i = ARC_HEAD(i->blossom_sibling); + if (i != i->blossom_parent->blossom_ptr) break; + i = i->blossom_parent; + } + if (i->is_outer) + { + lca->AddRoot(i); + break; + } + } + } + b->lca = lca; + } +#endif + + while ((a=selfloop_first)) + { + selfloop_first = a->next[1]; + do + { + Edge* a_next = a->next[0]; + +#ifdef LCA_REPAIRS + int _i = a->head0[1]->lca_preorder; + int _j = a->head0[0]->lca_preorder; + Node* b = a->head0[1]->blossom_grandparent; + b->lca->GetPenultimateNodes(_i, _j); + i = b->lca->rev_mapping[_i]; + j = b->lca->rev_mapping[_j]; +#else + GetRealEndpoints(a, i, j); +#endif + ADD_EDGE(i, a, 0); + ADD_EDGE(j, a, 1); + a->slack -= 2*i->blossom_eps; + a = a_next; + } while ( a ); + } + + /* + for (i0=nodes; i0is_outer) continue; + b = i0->blossom_grandparent; + if (b->lca) + { + delete b->lca; + b->lca = NULL; + } + } + */ + + nodes[node_num].first_tree_child = NULL; +} + +void PerfectMatching::FinishUpdate() +{ + Node* i0; + Node* i; + Node* j; + Edge* a; + EdgeIterator I; + int dir; + Tree* t; + + for (i0=nodes; i0is_outer) continue; + +#ifdef LCA_REPAIRS + if (i0->blossom_grandparent->lca) + { + delete i0->blossom_grandparent->lca; + i0->blossom_grandparent->lca = NULL; + } +#endif + + ////////////////////////////////////////////////////////////// + if (!i0->blossom_grandparent->is_removed) + { + i = i0; + do + { + i->y = ARC_TO_EDGE_PTR(i->blossom_sibling)->y_saved; + i->is_marked = 0; + i->blossom_selfloops = NULL; + i = i->blossom_parent; + } while (i->is_marked); + continue; + } + ////////////////////////////////////////////////////////////// + + i = i0->blossom_parent; + while ( 1 ) + { + if (i->is_removed && !i->is_outer) break; + REAL y_parent = (i->is_outer) ? 0 : i->blossom_parent->y; + for (dir=0; dir<2; dir++) + { + if (!i->first[dir]) continue; + i->first[dir]->prev[dir]->next[dir] = NULL; + Edge* a_next; + for (a=i->first[dir]; a; a=a_next) + { + a_next = a->next[dir]; + j = a->head0[1-dir]; + ADD_EDGE(j, a, dir); + a->slack += j->blossom_parent->y - y_parent; + } + i->first[dir] = NULL; + } + if (i->is_removed) break; + + j = i->blossom_parent; + i->is_removed = 1; + i->tree_sibling_next = removed_first; + removed_first = i; + i = j; + } + i0->y = ARC_TO_EDGE_PTR(i0->blossom_sibling)->y_saved; + i0->is_outer = 1; + i0->flag = 2; + i0->is_tree_root = 1; + } + + Node* blossom_list = nodes[node_num].first_tree_child; + + + + for (i=nodes; iis_tree_root) continue; + i->first_tree_child = nodes[node_num].first_tree_child; + nodes[node_num].first_tree_child = i; + REAL slack_min = PM_INFTY; + FOR_ALL_EDGES(i, a, dir, I) + { + if (slack_min > a->slack) slack_min = a->slack; + } + i->y += slack_min; + FOR_ALL_EDGES(i, a, dir, I) a->slack -= slack_min; + } + + tree_num = 0; + for (i=nodes[node_num].first_tree_child; i!=blossom_list; i=i->first_tree_child) + { + tree_num ++; + if (!i->is_tree_root) continue; + FOR_ALL_EDGES(i, a, dir, I) + { + j = a->head[dir]; + if (a->slack <= 0 && j->is_tree_root) + { + i->is_tree_root = j->is_tree_root = 0; + i->match = EDGE_DIR_TO_ARC(a, dir); + j->match = EDGE_DIR_TO_ARC(a, 1-dir); + tree_num -= 2; + break; + } + } + } + for ( ; i; i=i->first_tree_child) + { + if (i->is_removed) { i->is_tree_root = 0; continue; } + tree_num ++; + } + + if (tree_num > tree_num_max) + { + if (trees) free(trees); + tree_num_max = tree_num; + trees = (Tree*) malloc(tree_num_max*sizeof(Tree)); + } + t = trees; + + Node* last_root = &nodes[node_num]; + Node* i_next; + for (i=nodes; i; i=i_next) + { + if (!i->is_blossom) i_next = (ifirst_tree_child; + if (!i->is_tree_root) continue; + + i->flag = 0; + i->first_tree_child = NULL; + i->tree_sibling_prev = last_root; + last_root->tree_sibling_next = i; + last_root = i; + i->tree = t; + t->root = i; + t->eps = 0; + t->first[0] = t->first[1] = NULL; + t->pq_current = NULL; + t->pq00.Reset(); + t->pq0.Reset(); + t->pq_blossoms.Reset(); + t ++; + } + + assert(t == trees + tree_num); + last_root->tree_sibling_next = NULL; + + while ((i=removed_first)) + { + removed_first = i->tree_sibling_next; + blossoms->Delete(i); + blossom_num --; + } +} + +PerfectMatching::REAL PerfectMatching::GetTwiceSum(NodeId i) +{ + assert(i>=0 && ihead0[dir]; + if (i->is_outer) + { + if (!i->is_tree_root) + { + i->is_tree_root = 1; + i = ARC_HEAD(i->match); + assert(!i->is_tree_root && i->is_outer); + i->is_tree_root = 1; + if (i->is_blossom) + { + i->first_tree_child = nodes[node_num].first_tree_child; + nodes[node_num].first_tree_child = i; + } + } + return; + } + if (i->blossom_grandparent->is_removed) return; + } + + Node* b = i->blossom_grandparent; + assert(b->is_outer); + + if (!b->is_tree_root) + { + b->is_tree_root = 1; + i = ARC_HEAD(b->match); + assert(!i->is_tree_root && i->is_outer); + i->is_tree_root = 1; + if (i->is_blossom) + { + i->first_tree_child = nodes[node_num].first_tree_child; + nodes[node_num].first_tree_child = i; + } + } + + b->is_removed = 1; + b->tree_sibling_next = removed_first; + removed_first = b; +} + +PerfectMatching::EdgeId PerfectMatching::AddNewEdge(NodeId _i, NodeId _j, REAL cost, bool do_not_add_if_positive_slack) +{ + assert(_i>=0 && _i=0 && _j= edge_num_max) ReallocateEdges(); + Node* i = nodes + _i; + Node* j = nodes + _j; + Edge* a = edges + edge_num; + + a->slack = cost*COST_FACTOR; + a->head0[0] = j; + a->head0[1] = i; + Node* bi = (i->is_outer) ? i : i->blossom_grandparent; + Node* bj = (j->is_outer) ? j : j->blossom_grandparent; + if (bi == bj) + { +#ifdef LCA_REPAIRS + int _i = i->lca_preorder; + int _j = j->lca_preorder; + bi->lca->GetPenultimateNodes(_i, _j); + i = bi->lca->rev_mapping[_i]; + j = bi->lca->rev_mapping[_j]; +#else + GetRealEndpoints(a, i, j); +#endif + a->slack += i->blossom_parent->y + j->blossom_parent->y; + } + else + { + i = bi; + j = bj; + } + a->slack -= a->head0[0]->y + a->head0[1]->y; + + if (a->slack >= 0 && do_not_add_if_positive_slack) return -1; + + ADD_EDGE(i, a, 0); + ADD_EDGE(j, a, 1); + PriorityQueue::ResetItem(a); + + if (a->slack < 0) + { + ProcessNegativeEdge(a); + } + + return edge_num ++; +} + +void PerfectMatching::UpdateCost(EdgeId e, REAL delta_cost) +{ + assert(e>=0 && eslack += delta_cost*COST_FACTOR; + if (a->slack == 0) return; + if (a->slack > 0) + { + Node* i = a->head[1]; + Node* j = a->head[0]; + if (i->is_outer) + { + if (ARC_TO_EDGE_PTR(i->match) != a && ARC_TO_EDGE_PTR(j->match) != a) return; + } + else + { + if (ARC_TO_EDGE_PTR(i->blossom_sibling) != a && ARC_TO_EDGE_PTR(j->blossom_sibling) != a) return; + } + } + ProcessNegativeEdge(a); +} + diff --git a/blossom5-v2.05.src/PMshrink.cpp b/blossom5-v2.05.src/PMshrink.cpp new file mode 100644 index 0000000..eabad1a --- /dev/null +++ b/blossom5-v2.05.src/PMshrink.cpp @@ -0,0 +1,318 @@ +#include +#include +#include +#include "PMimplementation.h" + + +PerfectMatching::Node* PerfectMatching::FindBlossomRoot(Edge* a0) +{ + Node* i; + Node* j; + Node* _i[2]; + Node* r; + int branch; + + _i[0] = ARC_HEAD(a0); + _i[1] = ARC_TAIL(a0); + branch = 0; + while ( 1 ) + { + if (_i[branch]->is_marked) + { + r = _i[branch]; + j = _i[1-branch]; + break; + } + _i[branch]->is_marked = 1; + if (_i[branch]->is_tree_root) + { + j = _i[branch]; + i = _i[1-branch]; + while (!i->is_marked) + { + i->is_marked = 1; + i = ARC_HEAD(i->match); + GET_TREE_PARENT(i, i); + } + r = i; + break; + } + i = ARC_HEAD(_i[branch]->match); + GET_TREE_PARENT(i, _i[branch]); + branch = 1 - branch; + } + i = r; + while ( i != j ) + { + i = ARC_HEAD(i->match); + i = ARC_HEAD(i->tree_parent); + i->is_marked = 0; + } + // clear is_marked and is_outer + i = ARC_HEAD(a0); + while (i != r) + { + i->is_marked = 0; + i->is_outer = 0; + i = ARC_HEAD(i->match); + i->is_outer = 0; + i = ARC_HEAD(i->tree_parent); + } + i = ARC_TAIL(a0); + while (i != r) + { + i->is_marked = 0; + i->is_outer = 0; + i = ARC_HEAD(i->match); + i->is_outer = 0; + i = ARC_HEAD(i->tree_parent); + } + r->is_marked = 0; + r->is_outer = 0; + + return r; +} + + +void PerfectMatching::Shrink(Edge* a0) +{ + //assert(a0->head[0]->is_outer && a0->head[1]->is_outer); + //assert(a0->head[0]->flag == 0 && a0->head[1]->flag == 0); + + double start_time = get_time(); + + int branch, dir; + Node* r; + Node* i; + Node* j; + Edge* a; + Edge** a_inner_ptr; + Arc* a_prev; + Node* b = blossoms->New(); + Edge* a_augment = NULL; + Edge* b_match; + + b->first[0] = b->first[1] = NULL; + + // set is_outer=0 for all nodes in the blossom + r = FindBlossomRoot(a0); + Tree* t = r->tree; + REAL eps = t->eps; + + b->first_tree_child = NULL; + i = ARC_HEAD(a0); + branch = 0; + while ( 1 ) + { + if (i == r && branch) break; + i->is_marked = 1; + if (i == r) + { + branch = 1; + i = ARC_TAIL(a0); + continue; + } + + // remove i from the list of children + REMOVE_FROM_TREE(i); + + // move children of i to the list of children of b + if (i->first_tree_child) + { + j = i->first_tree_child; + if (!b->first_tree_child) b->first_tree_child = j; + else + { + Node* j_last = j->tree_sibling_prev; + j->tree_sibling_prev = b->first_tree_child->tree_sibling_prev; + b->first_tree_child->tree_sibling_prev->tree_sibling_next = j; + b->first_tree_child->tree_sibling_prev = j_last; + } + } + + // go to parent + i = ARC_HEAD(i->match); + i->is_marked = 1; + if (i->is_blossom) + { + a = ARC_TO_EDGE_PTR(i->match); + t->pq_blossoms.Remove(a, pq_buf); + REAL tmp = a->slack; a->slack = i->y; i->y = tmp; + } + i = ARC_HEAD(i->tree_parent); + } + + // move children of r to the list of children of b + if (i->first_tree_child) + { + j = i->first_tree_child; + if (!b->first_tree_child) b->first_tree_child = j; + else + { + Node* j_last = j->tree_sibling_prev; + j->tree_sibling_prev = b->first_tree_child->tree_sibling_prev; + b->first_tree_child->tree_sibling_prev->tree_sibling_next = j; + b->first_tree_child->tree_sibling_prev = j_last; + } + } + + // init b + b->is_removed = 0; + b->is_outer = 1; + b->flag = 0; + b->is_blossom = 1; + b->is_tree_root = r->is_tree_root; + b->is_processed = 1; + b->tree = t; + b->y = -eps; + b->is_marked = 0; + + // replace r with b in the tree + b->tree_sibling_prev = r->tree_sibling_prev; + b->tree_sibling_next = r->tree_sibling_next; + Node* b_parent = NULL; + if (!b->is_tree_root) + { + b_parent = ARC_HEAD(r->match); GET_TREE_PARENT(b_parent, b_parent); + } + if (b->tree_sibling_prev->tree_sibling_next) b->tree_sibling_prev->tree_sibling_next = b; + else b_parent->first_tree_child = b; + if (b->tree_sibling_next) b->tree_sibling_next->tree_sibling_prev = b; + else if (b_parent) b_parent->first_tree_child->tree_sibling_prev = b; + + if (b->is_tree_root) + { + b->tree->root = b; + b_match = NULL; + } + else + { + b->match = r->match; + b_match = ARC_TO_EDGE_PTR(b->match); + } + REAL b_match_slack = 0; // initialize to prevent compiler warning + if (b_match && ARC_HEAD(b->match)->is_blossom) + { + b_match_slack = b_match->slack; + b_match->slack = ARC_HEAD(b->match)->y; + } + + // second pass over nodes in the blossom + branch = 0; + a_prev = EDGE_DIR_TO_ARC(a0, 0); + i = ARC_HEAD(a_prev); + while ( 1 ) + { + // update Arc::next and Arc::head pointers + if (i->flag == 0) i->y += eps; + else i->y -= eps; + i->is_processed = 0; + + if (i->flag == 1) + { + Edge* a_prev; + for (dir=0; dir<2; dir++) + if (i->first[dir]) + { + for (a_inner_ptr=&i->first[dir], a=*a_inner_ptr, a_prev=a->prev[dir], a_prev->next[dir]=NULL; a; a=*a_inner_ptr) + { + Node* j0 = a->head[dir]; + for (j=j0; !j->is_outer && !j->is_marked; j = j->blossom_parent) {} + if (j != j0) { /*assert(j->flag == 0);*/ int dir_rev = 1 - dir; MOVE_EDGE(j0, j, a, dir_rev); } + if (j->is_marked) // "inner" arc + { + a_inner_ptr = &a->next[dir]; + a->prev[dir] = a_prev; + a_prev = a; + + if (j->flag == 1) a->slack += eps; + } + else // "boundary" arc + { + *a_inner_ptr = a->next[dir]; + ADD_EDGE(b, a, dir); + + if (j->flag == 0 && j->tree != t) + { + j->tree->pq_current->pq01[1-j->tree->dir_current].Remove(a, pq_buf); + if (a->slack + eps <= j->tree->eps) a_augment = a; + } + a->slack += 2*eps; + if (j->flag == 2) t->pq0.Add(a); + else if (j->flag == 0) + { + if (!j->tree->pq_current) AddTreeEdge(t, j->tree); + j->tree->pq_current->pq00.Add(a); + } + else if (j->tree != t) + { + if (!j->tree->pq_current) AddTreeEdge(t, j->tree); + j->tree->pq_current->pq01[j->tree->dir_current].Add(a); + } + } + } + if (i->first[dir]) + { + a_prev->next[dir] = i->first[dir]; + i->first[dir]->prev[dir] = a_prev; + } + } + } + + Arc* a_next = (i->flag == 0) ? i->match : i->tree_parent; + i->blossom_parent = b; + i->match = NULL; + i->blossom_grandparent = b; + i->blossom_selfloops = NULL; + if (branch == 0) + { + i->blossom_sibling = a_next; + if (i == r) + { + branch = 1; + a_prev = ARC_REV(a0); + i = ARC_HEAD(a_prev); + if (i == r) break; + } + else + { + a_prev = i->blossom_sibling; + i = ARC_HEAD(a_prev); + } + } + else + { + i->blossom_sibling = ARC_REV(a_prev); + a_prev = a_next; + i = ARC_HEAD(a_prev); + if (i == r) break; + } + } + i->blossom_sibling = ARC_REV(a_prev); + r->is_tree_root = 0; + + for (i=ARC_HEAD(r->blossom_sibling); ; i = ARC_HEAD(i->blossom_sibling)) + { + i->is_marked = 0; + i->blossom_eps = eps; + if (i == r) break; + } + + if (b_match) + { + if (ARC_HEAD(b->match)->is_blossom) + { + b_match->slack = b_match_slack; + } + dir = ARC_TO_EDGE_DIR(b->match); + //assert(b_match->head[1-dir] == r); + MOVE_EDGE(r, b, b_match, dir); + } + + stat.shrink_count ++; + blossom_num ++; + stat.shrink_time += get_time() - start_time; + + if (a_augment) Augment(a_augment); +} + diff --git a/blossom5-v2.05.src/PQ.h b/blossom5-v2.05.src/PQ.h new file mode 100644 index 0000000..8fb37b7 --- /dev/null +++ b/blossom5-v2.05.src/PQ.h @@ -0,0 +1,396 @@ +/* + PQ.h - implements pairing heaps priority queue (with multipass 'delete min') + + Copyright 2008 Vladimir Kolmogorov (vnk@ist.ac.at) + + 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 HFKSJHFKJHARBABDAKFAF +#define HFKSJHFKJHARBABDAKFAF + +// exactly one flag must be defined +//#define PQ_MULTIPASS +#define PQ_INTERLEAVED_MULTIPASS + +#include + +template class PriorityQueue +{ +public: + struct Item + { + REAL slack; + + Item* parentPQ; + union + { + struct + { + Item* leftPQ; + Item* rightPQ; + }; + REAL y_saved; // used in repairs + }; + }; + static void* AllocateBuf(); + static void DeallocateBuf(void* buf); + + static void ResetItem(Item* i); + static bool isReset(Item* i); + + ////////////////////////////////////////////////////////// + + void Reset(); + void Add(Item* i); +#define Remove(i, buf) _Remove(i) + void _Remove(Item* i); + void Decrease(Item* i_old, Item* i_new, void* buf); + Item* GetMin(); + + ////////////////////////////////////////////////////////// + + void Update(REAL delta); + void Merge(PriorityQueue& dest); + + // traversing items in the order they are stored (irrespective of slack). + // The caller must go through all items, no other member functions can be called during the scan. + Item* GetAndResetFirst(); + Item* GetAndResetNext(); + + Item* GetFirst(); + Item* GetNext(Item* i); + + ////////////////////////////////////////////////////////// + +private: + struct Buf + { + }; + Item* rootPQ; + void RemoveRoot(); +}; + +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// + +template inline void* PriorityQueue::AllocateBuf() +{ + return NULL; +} + +template inline void PriorityQueue::DeallocateBuf(void* _buf) +{ +} + +template inline void PriorityQueue::ResetItem(Item* i) +{ + i->parentPQ = NULL; +} + +template inline bool PriorityQueue::isReset(Item* i) +{ + return (i->parentPQ == NULL); +} + +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// + +template inline void PriorityQueue::Reset() +{ + rootPQ = NULL; +} + +/* +template inline void PriorityQueue::RemoveRoot() +{ + Item* r = rootPQ; + PriorityQueue pq; + pq.rootPQ = rootPQ; + rootPQ = NULL; + Item* i; + for (i=pq.GetAndResetFirst(); i; i=pq.GetAndResetNext()) + { + if (i != r) Add(i); + } + r->parentPQ = NULL; +} +*/ + +// sets i = merge(i, j). Ignores parentPQ and rightPQ for i and j. +#define MERGE_PQ(i, j)\ + {\ + if (i->slack <= j->slack)\ + {\ + j->rightPQ = i->leftPQ;\ + if (j->rightPQ) j->rightPQ->parentPQ = j;\ + j->parentPQ = i;\ + i->leftPQ = j;\ + }\ + else\ + {\ + i->rightPQ = j->leftPQ;\ + if (i->rightPQ) i->rightPQ->parentPQ = i;\ + i->parentPQ = j;\ + j->leftPQ = i;\ + i = j;\ + }\ + } + +template inline void PriorityQueue::RemoveRoot() +{ + Item* i = rootPQ->leftPQ; + rootPQ->parentPQ = NULL; + if (i) + { +#ifdef PQ_MULTIPASS + while ( i->rightPQ ) + { + Item** prev_ptr = &rootPQ; + while ( 1 ) + { + if (i->rightPQ) + { + Item* j = i->rightPQ; + Item* next = j->rightPQ; + MERGE_PQ(i, j); + *prev_ptr = i; + if (!next) { i->rightPQ = NULL; break; } + prev_ptr = &i->rightPQ; + i = next; + } + else + { + *prev_ptr = i; + i->rightPQ = NULL; + break; + } + } + i = rootPQ; + } +#endif + +#ifdef PQ_INTERLEAVED_MULTIPASS + while ( i->rightPQ ) + { + Item* prev = NULL; + while ( i ) + { + Item* next; + if (i->rightPQ) + { + Item* j = i->rightPQ; + next = j->rightPQ; + MERGE_PQ(i, j); + } + else next = NULL; + i->rightPQ = prev; + prev = i; + i = next; + } + i = prev; + } +#endif + i->parentPQ = i; + } + rootPQ = i; +} + +template inline void PriorityQueue::Add(Item* i) +{ + if (!rootPQ) + { + rootPQ = i; + i->parentPQ = i; + i->leftPQ = i->rightPQ = NULL; + } + else if (i->slack <= rootPQ->slack) + { + rootPQ->parentPQ = i; + i->leftPQ = rootPQ; + i->rightPQ = NULL; + rootPQ = i; + i->parentPQ = i; + } + else + { + i->leftPQ = NULL; + i->rightPQ = rootPQ->leftPQ; + if (i->rightPQ) i->rightPQ->parentPQ = i; + rootPQ->leftPQ = i; + i->parentPQ = rootPQ; + } +} + + +template inline void PriorityQueue::_Remove(Item* i) +{ + Item* p = i->parentPQ; + if (p == i) RemoveRoot(); + else + { + if (i->rightPQ) i->rightPQ->parentPQ = p; + if (p->leftPQ == i) p->leftPQ = i->rightPQ; + else p->rightPQ = i->rightPQ; + if (i->leftPQ) + { + i->parentPQ = i; + i->rightPQ = NULL; + PriorityQueue pq; + pq.rootPQ = i; + pq.RemoveRoot(); + pq.Merge(*this); + } + else i->parentPQ = NULL; + } +} + +template inline void PriorityQueue::Decrease(Item* i_old, Item* i_new, void* _buf) +{ + if (i_old->parentPQ == i_old) + { + if (i_old != i_new) + { + rootPQ = i_new; + i_new->parentPQ = i_new; + i_new->leftPQ = i_old->leftPQ; + i_new->rightPQ = NULL; + if (i_new->leftPQ) i_new->leftPQ->parentPQ = i_new; + i_old->parentPQ = NULL; + } + } + else + { + Remove(i_old, _buf); + Add(i_new); + } +} + +template inline typename PriorityQueue::Item* PriorityQueue::GetMin() +{ + return rootPQ; +} + +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// + + + +template inline void PriorityQueue::Merge(PriorityQueue& dest) +{ + if (!rootPQ) return; + if (!dest.rootPQ) dest.rootPQ = rootPQ; + else + { + if (rootPQ->slack < dest.rootPQ->slack) + { + Item* j = rootPQ; rootPQ = dest.rootPQ; dest.rootPQ = j; + } + rootPQ->rightPQ = dest.rootPQ->leftPQ; + if (rootPQ->rightPQ) rootPQ->rightPQ->parentPQ = rootPQ; + rootPQ->parentPQ = dest.rootPQ; + dest.rootPQ->leftPQ = rootPQ; + } + rootPQ = NULL; +} + + + +template inline void PriorityQueue::Update(REAL delta) +{ + if (!rootPQ) return; + + Item* i = rootPQ; + while (i->leftPQ) i = i->leftPQ; + + while ( 1 ) + { + i->slack += delta; + + if (i->rightPQ) + { + i = i->rightPQ; + while (i->leftPQ) i = i->leftPQ; + } + else + { + while ( 1 ) + { + Item* j = i; + i = i->parentPQ; + if (i == j) return; + if (i->leftPQ == j) break; + } + } + } +} + + +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// + +template inline typename PriorityQueue::Item* PriorityQueue::GetAndResetFirst() +{ + if (!rootPQ) return NULL; + return GetAndResetNext(); +} + +template inline typename PriorityQueue::Item* PriorityQueue::GetAndResetNext() +{ + if (!rootPQ) return NULL; + Item* result = rootPQ; + result->parentPQ = NULL; + Item* i = rootPQ->leftPQ; + if (!i) rootPQ = result->rightPQ; + else + { + rootPQ = i; + while (i->rightPQ) i = i->rightPQ; + i->rightPQ = result->rightPQ; + } + return result; +} + +template inline typename PriorityQueue::Item* PriorityQueue::GetFirst() +{ + if (!rootPQ) return NULL; + Item* i = rootPQ; + while (i->leftPQ) i = i->leftPQ; + return i; +} + +template inline typename PriorityQueue::Item* PriorityQueue::GetNext(Item* i) +{ + if (i->rightPQ) + { + i = i->rightPQ; + while (i->leftPQ) i = i->leftPQ; + return i; + } + while ( 1 ) + { + Item* j = i; + i = i->parentPQ; + if (i == j) return NULL; + if (i->leftPQ == j) return i; + } +} + +#endif diff --git a/blossom5-v2.05.src/PerfectMatching.h b/blossom5-v2.05.src/PerfectMatching.h new file mode 100644 index 0000000..7b1ae57 --- /dev/null +++ b/blossom5-v2.05.src/PerfectMatching.h @@ -0,0 +1,246 @@ +/* + PerfectMatching.h - interface to min cost perfect matching code + + + + Copyright 2008-2009 UCL Business PLC, Author Vladimir Kolmogorov (vnk@ist.ac.at) + + This software can be used for evaluation and non-commercial research purposes only. Commercial use is prohibited. + Public redistribution of the code or its derivatives is prohibited. + If you use this software for research purposes, you should cite the following paper in any resulting publication: + + Vladimir Kolmogorov. "Blossom V: A new implementation of a minimum cost perfect matching algorithm." + In Mathematical Programming Computation (MPC), July 2009, 1(1):43-67. + + For commercial use of the software not covered by this agreement, you may obtain a licence from + the copyright holders UCL Business via their licensing site: www.e-lucid.com/i/software/optimisation_software/BlossomV.html. + + 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 HALSKDJDFHALSJASFDFASJGLA +#define HALSKDJDFHALSJASFDFASJGLA + +#include +#include +#include "block.h" + + +// if defined, edge costs are of type 'double', otherwise 'int' +//#define PERFECT_MATCHING_DOUBLE + +// Note: with floating point numbers polynomial complexity is not guaranteed; +// the code may even get stuck due to rounding errors. If the code terminates, +// the solution may not be optimal. It may be worth calling CheckPerfectMatchingOptimality() +// to see whether complementary slackness conditions are satisfied. +// +// Using single precision floating point numbers (float) is really not recommended. + + +class PerfectMatching +{ +public: + +#ifdef PERFECT_MATCHING_DOUBLE + typedef double REAL; + #define PM_INFTY ((REAL)1e100) +#else + typedef int REAL; + #define PM_INFTY (INT_MAX/2) +#endif + + typedef int NodeId; + typedef int EdgeId; + + PerfectMatching(int nodeNum, int edgeNumMax); + ~PerfectMatching(); + + // first call returns 0, second 1, and so on. + EdgeId AddEdge(NodeId i, NodeId j, REAL cost); + + // Computes a perfect matching of minimum cost. + // NOTE: a perfect matching of finite cost must exist (otherwise the behaviour is not specified). + // If finish is false, then the final matching is not computed (so GetSolution() cannot be called afterwards). + void Solve(bool finish=true); + + /////////////////////////////////////////////////////////////// + // Read primal solution (can be called after Solve()). + int GetSolution(EdgeId e); // returns 1 if e is in the matching, 0 otherwise + NodeId GetMatch(NodeId i); // alternative way to get the result + + /////////////////////////////////////////////////////////////// + // Read dual solution (can be called after Solve()). + // 'blossom_parents' and 'twice_y' must be arrays of size node_num+GetBlossomNum(). + // The function sets blossom_parent[i] to the parent of i (or to -1 for exterior nodes). + void GetDualSolution(int* blossom_parents, REAL* twice_y); + int GetBlossomNum(); + + /////////////////////////////////////////////////////////////// + // Dynamic graph updates. After calling Solve() you may call // + // StartUpdate(), ..., FinishUpdate() and then Solve() again // + /////////////////////////////////////////////////////////////// + void StartUpdate(); + void FinishUpdate(); + + // 3 functions below can be called only between StartUpdate() and FinishUpdate(). + REAL GetTwiceSum(NodeId i); // if 2*cost(i,j)>=GetTwiceSum(i)+GetTwiceSum(j) then adding new edge (i,j) is not necessary - optimal solution will not change + EdgeId AddNewEdge(NodeId i, NodeId j, REAL cost, bool do_not_add_if_positive_slack=true); // if do_not_add_if_positive_slack is true and the slack of the edge turns out to be non-negative, then the edge will not be added and -1 will be returned + void UpdateCost(EdgeId e, REAL delta_cost); + + + // NOTE: with default options the dual vector is guaranteed to be half-integral + // (if all input weights are integral). However, with some options there is + // no such guarantee, in particular, if dual_greedy_update_option=2 or dual_LP_threshold>0. + // These options can be used only if the code is compiled with REAL=double. + struct Options + { + Options() : fractional_jumpstart(true), + dual_greedy_update_option(0), + dual_LP_threshold(0.00), + update_duals_before(false), + update_duals_after(false), + single_tree_threshold(1.00), + verbose(true) + {} + + bool fractional_jumpstart; // false: greedy, true: compute fractional matching + + int dual_greedy_update_option; // 0: compute connected components (as in Blossom IV) + // 1: compute strongly connected components (discussed by Cook-Rohe, but not implemented) + // 2: single eps for all trees (fixed eps approach) + + double dual_LP_threshold; // if tree_num => dual_updates_threshold*node_num: greedy updates + // if tree_num < dual_updates_threshold*node_num: global updates (solve LP) + + bool update_duals_before; // before tree growth + bool update_duals_after; // after tree growth + + double single_tree_threshold; // if (tree_num => single_tree_threshold*node_num && update_duals_after): try to grow a single tree as long as possible + + bool verbose; + } options; + + + // save problem to a file. format=0 corresponds to DIMACS format, + // format=1 corresponds to the format used by blossom4. + // CANNOT BE CALLED AFTER Solve()!!! + void Save(char* filename, int format=0); + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// +private: + struct Node; + struct Arc; // no such struct, only Arc* is used (the pointer can be odd or even) + struct Edge; // pointer Edge* is always even + struct Tree; + struct TreeEdge; + struct PQPointers; + struct EdgeIterator; + struct TreeEdgeIterator; + struct LCATreeX; + + Node* nodes; + Edge* edges; + char* edges_orig; + DBlock* blossoms; + Tree* trees; + DBlock* tree_edges; + struct ExpandTmpItem + { + Node* i; + Node* blossom_parent; + Node* blossom_grandparent; + }; + Block* expand_tmp_list; // used inside Expand() + + int node_num; + int edge_num, edge_num_max; + int tree_num, tree_num_max; + + Node* removed_first; + int blossom_num; + int removed_num; + + void* pq_buf; + + bool first_solve; + + // stat + struct Stat + { + int shrink_count; + int expand_count; + int grow_count; + double shrink_time; + double expand_time; + double dual_time; + } stat; + + //////////////////////////////////////////////////////////////////// + + void InitGreedy(bool allocate_trees=true); + + void InitGlobal(); // compute fractional matching + Node* FindBlossomRootInit(Edge* a0); + void ShrinkInit(Edge* a0, Node* tree_root); + void ExpandInit(Node* b); + void AugmentBranchInit(Node* i0, Node* tree_root); + + void Finish(); // sets matching for inner nodes + + void ProcessNegativeEdge(Edge* a); + + void GetRealEndpoints(Edge* a, Node*& tail, Node*& head); + Node* FindBlossomRoot(Edge* a0); + void Shrink(Edge* a0); + void Expand(Node* b); + void Augment(Edge* a0); + void AugmentBranch(Node* i0); + void GrowNode(Node* i); + void GrowTree(Node* root, bool new_subtree); + bool ProcessEdge00(Edge* a, bool update_boundary_edge=true); // returns true if boundary edge, false otherwise + void ProcessSelfloop(Node* b, Edge* a); + + void AddTreeEdge(Tree* t0, Tree* t1); + + void ComputeEpsSingle(); // called from UpdateDuals() + void ComputeEpsCC(); // called from UpdateDuals() + void ComputeEpsSCC(); // called from UpdateDuals() + void ComputeEpsGlobal(); // called from UpdateDuals() + bool UpdateDuals(); + + void FreeRemoved(); + void CommitEps(); + + void ReallocateEdges(); + + void PrintAll(); +}; + + +// in the functions below, 'edges' is an array of size 2*edge_num (edge e = (edges[2*e],edges[2*e+1])), +// 'weights' is an array of size edge_num + +// checks complementary slackness conditions. +// returns 0 if success. +// returns 1 if complementary slackness conditions are violated (then the amount of violation is printed - could potentially happen for double's) +// returns 2 if the blossom tree structure is incorrect (or inconsistent with primal solution) +int CheckPerfectMatchingOptimality(int node_num, int edge_num, int* edges, int* weights, PerfectMatching* pm, PerfectMatching::REAL threshold=(PerfectMatching::REAL)(1e-10)); + + +double ComputePerfectMatchingCost(int node_num, int edge_num, int* edges, int* weights, PerfectMatching* pm); + +#endif diff --git a/blossom5-v2.05.src/README.TXT b/blossom5-v2.05.src/README.TXT new file mode 100644 index 0000000..2e65742 --- /dev/null +++ b/blossom5-v2.05.src/README.TXT @@ -0,0 +1,80 @@ +BLOSSOM V - implementation of Edmonds' algorithm for computing a minimum cost perfect matching in a graph +Version 2.05 +http://pub.ist.ac.at/~vnk/software.html + +Details of the implementation are described in + + Vladimir Kolmogorov. "Blossom V: A new implementation of a minimum cost perfect matching algorithm." + In Mathematical Programming Computation (MPC), July 2009, 1(1):43-67. + +Please send comments to vnk@ist.ac.at. +If you use this software for research purposes, you should cite the aforementioned paper in any resulting publication. + +################################################################## + +License & disclaimer: + + Copyright 2008-2009 UCL Business PLC, Author Vladimir Kolmogorov (vnk@ist.ac.at) + + This software can be used for evaluation and non-commercial research purposes only. Commercial use is prohibited. + Public redistribution of the code or its derivatives is prohibited. + If you use this software for research purposes, you should cite the following paper in any resulting publication: + + Vladimir Kolmogorov. "Blossom V: A new implementation of a minimum cost perfect matching algorithm." + In Mathematical Programming Computation (MPC), July 2009, 1(1):43-67. + + For commercial use of the software not covered by this agreement, you may obtain a licence from + the copyright holders UCL Business via their licensing site: www.e-lucid.com/i/software/optimisation_software/BlossomV.html. + + 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. + + +################################################################## + +COMPILATION: + +In unix, type "make". The code should also compile in Windows with Microsoft Visual C++ compiler. +Tested on 32-bit machines. + +################################################################## + +USAGE: + +See PerfectMatching.h for interface functions. Alternatively, +compile the code and run ./blossom5 [options] (see USAGE.TXT for details). + +The code also allows solving complete geometric instances; see GEOM/GeomPerfectMatching.h +for interface functions. + +################################################################## + +PARAMETERS: + +For many types of problems, default parameters should be ok. +But for certain problems (such as structured geometric instances) +you may consider setting +dual_LP_threshold=0.005 (for example). This corresponds to calling +./blossom5 -m0.005 +Type PerfectMatching::REAL in PerfectMatcing.h should then be set to double. + +################################################################## + +EXTERNAL PACKAGE: + +When solving complete geometric instances you need to provide the initial subset of edges. +It may be desirable to use Delaunay triangulation for this purpose. +Then you need to download the "Triangle" package of J. R. Shewchuk from + http://www.cs.cmu.edu/~quake/triangle.html +and extract it to the directory "triange". +Alternatively, you can use nearest neighbours initialization (this does not require external packages). + diff --git a/blossom5-v2.05.src/USAGE.TXT b/blossom5-v2.05.src/USAGE.TXT new file mode 100644 index 0000000..afe93cd --- /dev/null +++ b/blossom5-v2.05.src/USAGE.TXT @@ -0,0 +1,31 @@ +USAGE: ./blossom5 [options] + +At least one of the two flags below must be provided. +If both are provided, then edges in are added to the initial subset of edges for the geometric problem. + + -e read problem (list of edges) from file in DIMACS or blossom4 format + (see GRAPH1.TXT and GRAPH2.TXT for examples) + -g read geometric problem (list of 2D points) from file in TSPLIB format + +perfect matching parameters: + + -j do not use fractional jumpstart + -m# update duals by solving LP if the number of trees is smaller than # times node_num. + (# should belong to [0,1], default is 0.) + -d1 use SCC dual updates instead of CC + -d2 use dual updates with fixed delta instead of CC + -b update duals before processing a tree + -a update duals after processing a tree + +geometric matching parameters: + + -D do NOT add Delaunay triangulation to the initial subset of edges + -K# for each point, add # nearest neighbors to the initial subset of edges + -I do NOT add edges greedily to the initial subset of edges to make sure that a perfect matching exists + -T# use at most # iterations (or run until convergence, if #=0 - default) + +other options: + + -w save result to file + -c check complementary slackness conditions for non-geometric problem + -V no verbose diff --git a/blossom5-v2.05.src/block.h b/blossom5-v2.05.src/block.h new file mode 100644 index 0000000..3bed9d4 --- /dev/null +++ b/blossom5-v2.05.src/block.h @@ -0,0 +1,268 @@ +/* block.h */ +/* + Template classes Block and DBlock + Implement adding and deleting items of the same type in blocks. + + If there there are many items then using Block or DBlock + is more efficient than using 'new' and 'delete' both in terms + of memory and time since + (1) On some systems there is some minimum amount of memory + that 'new' can allocate (e.g., 64), so if items are + small that a lot of memory is wasted. + (2) 'new' and 'delete' are designed for items of varying size. + If all items has the same size, then an algorithm for + adding and deleting can be made more efficient. + (3) All Block and DBlock functions are inline, so there are + no extra function calls. + + Differences between Block and DBlock: + (1) DBlock allows both adding and deleting items, + whereas Block allows only adding items. + (2) Block has an additional operation of scanning + items added so far (in the order in which they were added). + (3) Block allows to allocate several consecutive + items at a time, whereas DBlock can add only a single item. + + Note that no constructors or destructors are called for items. + + Example usage for items of type 'MyType': + + /////////////////////////////////////////////////// + #include "block.h" + #define BLOCK_SIZE 1024 + typedef struct { int a, b; } MyType; + MyType *ptr, *array[10000]; + + ... + + Block *block = new Block(BLOCK_SIZE); + + // adding items + for (int i=0; i New(); + ptr -> a = ptr -> b = rand(); + } + + // reading items + for (ptr=block->ScanFirst(); ptr; ptr=block->ScanNext()) + { + printf("%d %d\n", ptr->a, ptr->b); + } + + delete block; + + ... + + DBlock *dblock = new DBlock(BLOCK_SIZE); + + // adding items + for (int i=0; i New(); + } + + // deleting items + for (int i=0; i Delete(array[i]); + } + + // adding items + for (int i=0; i New(); + } + + delete dblock; + + /////////////////////////////////////////////////// + + Note that DBlock deletes items by marking them as + empty (i.e., by adding them to the list of free items), + so that this memory could be used for subsequently + added items. Thus, at each moment the memory allocated + is determined by the maximum number of items allocated + simultaneously at earlier moments. All memory is + deallocated only when the destructor is called. +*/ + +#ifndef __BLOCK_H__ +#define __BLOCK_H__ + +#include + +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ + +template class Block +{ +public: + /* Constructor. Arguments are the block size and + (optionally) the pointer to the function which + will be called if allocation failed; the message + passed to this function is "Not enough memory!" */ + Block(int size, void (*err_function)(const char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; } + + /* Destructor. Deallocates all items added so far */ + ~Block() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } } + + /* Allocates 'num' consecutive items; returns pointer + to the first item. 'num' cannot be greater than the + block size since items must fit in one block */ + Type *New(int num = 1) + { + Type *t; + + if (!last || last->current + num > last->last) + { + if (last && last->next) last = last -> next; + else + { + block *next = (block *) new char [sizeof(block) + (block_size-1)*sizeof(Type)]; + if (!next) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + if (last) last -> next = next; + else first = next; + last = next; + last -> current = & ( last -> data[0] ); + last -> last = last -> current + block_size; + last -> next = NULL; + } + } + + t = last -> current; + last -> current += num; + return t; + } + + /* Returns the first item (or NULL, if no items were added) */ + Type *ScanFirst() + { + for (scan_current_block=first; scan_current_block; scan_current_block = scan_current_block->next) + { + scan_current_data = & ( scan_current_block -> data[0] ); + if (scan_current_data < scan_current_block -> current) return scan_current_data ++; + } + return NULL; + } + + /* Returns the next item (or NULL, if all items have been read) + Can be called only if previous ScanFirst() or ScanNext() + call returned not NULL. */ + Type *ScanNext() + { + while (scan_current_data >= scan_current_block -> current) + { + scan_current_block = scan_current_block -> next; + if (!scan_current_block) return NULL; + scan_current_data = & ( scan_current_block -> data[0] ); + } + return scan_current_data ++; + } + + /* Marks all elements as empty */ + void Reset() + { + block *b; + if (!first) return; + for (b=first; ; b=b->next) + { + b -> current = & ( b -> data[0] ); + if (b == last) break; + } + last = first; + } + +/***********************************************************************/ + +private: + + typedef struct block_st + { + Type *current, *last; + struct block_st *next; + Type data[1]; + } block; + + int block_size; + block *first; + block *last; + + block *scan_current_block; + Type *scan_current_data; + + void (*error_function)(const char *); +}; + +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ + +template class DBlock +{ +public: + /* Constructor. Arguments are the block size and + (optionally) the pointer to the function which + will be called if allocation failed; the message + passed to this function is "Not enough memory!" */ + DBlock(int size, void (*err_function)(const char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; } + + /* Destructor. Deallocates all items added so far */ + ~DBlock() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } } + + /* Allocates one item */ + Type *New() + { + block_item *item; + + if (!first_free) + { + block *next = first; + first = (block *) new char [sizeof(block) + (block_size-1)*sizeof(block_item)]; + if (!first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + first_free = & (first -> data[0] ); + for (item=first_free; item next_free = item + 1; + item -> next_free = NULL; + first -> next = next; + } + + item = first_free; + first_free = item -> next_free; + return (Type *) item; + } + + /* Deletes an item allocated previously */ + void Delete(Type *t) + { + ((block_item *) t) -> next_free = first_free; + first_free = (block_item *) t; + } + +/***********************************************************************/ + +private: + + typedef union block_item_st + { + Type t; + block_item_st *next_free; + } block_item; + + typedef struct block_st + { + struct block_st *next; + block_item data[1]; + } block; + + int block_size; + block *first; + block_item *first_free; + + void (*error_function)(const char *); +}; + + +#endif + diff --git a/blossom5-v2.05.src/example.cpp b/blossom5-v2.05.src/example.cpp new file mode 100644 index 0000000..16c962e --- /dev/null +++ b/blossom5-v2.05.src/example.cpp @@ -0,0 +1,285 @@ +#include +#include "PerfectMatching.h" +#include "GEOM/GeomPerfectMatching.h" + +void LoadFile(int& node_num, int& edge_num, int*& edges, int*& weights, char* filename) +{ + int e = 0; + char LINE[1000]; + FILE* fp = fopen(filename, "r"); + if (!fp) { printf("Can't open %s\n", filename); exit(1); } + + int format = -1; // 0: DIMACS format. node id's start from 1 + // 1: simpler format (without "p" and "e"). node id's start from 0 + + edge_num = -1; + while (fgets(LINE, sizeof(LINE)-1, fp)) + { + if (LINE[0] == 'c') continue; + if (format < 0) + { + if (LINE[0] == 'p') + { + format = 0; + if (sscanf(LINE, "p edge %d %d\n", &node_num, &edge_num) != 2) { printf("%s: wrong format #1\n", filename); exit(1); } + } + else + { + format = 1; + if (sscanf(LINE, "%d %d\n", &node_num, &edge_num) != 2) { printf("%s: wrong format #1\n", filename); exit(1); } + } + + ////////////////////////////////////////////////////////////////////////////////// + if (node_num <= 0 || edge_num < 0) { printf("# of nodes and edges should be positive\n"); exit(1); } + if (node_num & 1) { printf("# of nodes is odd: perfect matching cannot exist\n"); exit(1); } + edges = new int[2*edge_num]; + weights = new int[edge_num]; + ////////////////////////////////////////////////////////////////////////////////// + } + else + { + int i, j; + char* ptr = LINE; + if (format == 0) { if (LINE[0] != 'e') continue; ptr = &LINE[1]; } + else ptr = &LINE[0]; + + int len; + if (sscanf(ptr, "%d %d %d\n", &i, &j, &len) != 3) continue; + if (format == 0) { i --; j --; } + edges[2*e] = i; + edges[2*e+1] = j; + weights[e] = len; + e ++; + } + } + + if (e != edge_num) { printf("%s: wrong format #3\n", filename); exit(1); } + fclose(fp); +} + +void SaveMatching(int node_num, PerfectMatching* pm, char* filename) +{ + FILE* fp = fopen(filename, "w"); + if (!fp) { printf("Can't open %s\n", filename); exit(1); } + fprintf(fp, "%d %d\n", node_num, node_num/2); + int i, j; + for (i=0; iGetMatch(i); + if (i < j) fprintf(fp, "%d %d\n", i, j); + } + fclose(fp); +} + +void LoadGeomFile(int& node_num, int*& x_array, int*& y_array, char* filename) +{ + int i = 0, i_tmp, x, y, DIM = 0; + char LINE[1000]; + + x_array = y_array = NULL; + + FILE* fp = fopen(filename, "r"); + if (!fp) { printf("Can't open %s\n", filename); exit(1); } + + while (fgets(LINE, sizeof(LINE)-1, fp)) + { + if (sscanf(LINE, "DIMENSION : %d", &node_num) == 1) + { + if (node_num < 1) { printf("too few nodes\n"); exit(1); } + if (node_num & 1) { printf("# of points is odd: perfect matching cannot exist\n"); exit(1); } + if (x_array) { printf("wrong format\n"); exit(1); } + x_array = new int[node_num]; + y_array = new int[node_num]; + continue; + } + if (sscanf(LINE, "EDGE_WEIGHT_TYPE : EUC_%dD", &DIM) == 1) + { + if (DIM != 2) { printf("only EUC_2D is supported"); exit(1); } + continue; + } + if (sscanf(LINE, "%d %d %d", &i_tmp, &x, &y) == 3) + { + i_tmp --; + if (i_tmp != i ++ || i > node_num) { printf("wrong number of points\n"); exit(1); } + x_array[i_tmp] = x; + y_array[i_tmp] = y; + continue; + } + printf("%s", LINE); + } + fclose(fp); + if (i != node_num || !x_array || DIM != 2) { printf("wrong format\n"); exit(1); } +} + + +void SaveMatching(int node_num, GeomPerfectMatching* gpm, char* filename) +{ + FILE* fp = fopen(filename, "w"); + if (!fp) { printf("Can't open %s\n", filename); exit(1); } + fprintf(fp, "%d %d\n", node_num, node_num/2); + int i, j; + for (i=0; iGetMatch(i); + if (i < j) + { + GeomPerfectMatching::REAL len = gpm->Dist(i, j); + if ( ((GeomPerfectMatching::REAL)1 / 2) == 0 ) fprintf(fp, "%d %d %d\n", i, j, len); + else fprintf(fp, "%d %d %f\n", i, j, (double)len); + } + } + fclose(fp); +} + + +void ShowUsage() +{ + printf("Usage: see USAGE.TXT\n"); + exit(1); +} + + +int main(int argc, char* argv[]) +{ + struct PerfectMatching::Options options; + struct GeomPerfectMatching::GPMOptions gpm_options; + char* filename = NULL; + char* geom_filename = NULL; + char* save_filename = NULL; + bool check_perfect_matching = false; + int i, e, node_num, edge_num; + int* edges; + int* weights; + + for (i=1; i1) ShowUsage(); + break; + case 'd': + options.dual_greedy_update_option = atoi(&argv[i][2]); + if (options.dual_greedy_update_option<1 || options.dual_greedy_update_option>2) ShowUsage(); + break; + case 'b': + if (argv[i][2]) ShowUsage(); + options.update_duals_before = true; + break; + case 'a': + if (argv[i][2]) ShowUsage(); + options.update_duals_after = true; + break; + case 'D': + if (argv[i][2]) ShowUsage(); + gpm_options.init_Delaunay = false; + break; + case 'K': + gpm_options.init_KNN = atoi(&argv[i][2]); + if (gpm_options.init_KNN<0) ShowUsage(); + break; + case 'I': + if (argv[i][2]) ShowUsage(); + gpm_options.init_greedy = false; + break; + case 'T': + gpm_options.iter_max = atoi(&argv[i][2]); + if (gpm_options.iter_max<0) ShowUsage(); + break; + case 'w': + if (save_filename || argv[i][2] || ++i == argc) ShowUsage(); + save_filename = argv[i]; + break; + case 'c': + if (argv[i][2]) ShowUsage(); + check_perfect_matching = true; + break; + case 'V': + if (argv[i][2]) ShowUsage(); + options.verbose = false; + break; + default: + printf("Unknown option: %s\n", argv[i]); + ShowUsage(); + break; + } + + } + + if (!filename && !geom_filename) ShowUsage(); + + if (filename) LoadFile(node_num, edge_num, edges, weights, filename); + + if (!geom_filename) + { + PerfectMatching *pm = new PerfectMatching(node_num, edge_num); + for (e=0; eAddEdge(edges[2*e], edges[2*e+1], weights[e]); + pm->options = options; + pm->Solve(); + if (check_perfect_matching) + { + int res = CheckPerfectMatchingOptimality(node_num, edge_num, edges, weights, pm); + printf("check optimality: res=%d (%s)\n", res, (res==0) ? "ok" : ((res==1) ? "error" : "fatal error")); + } + double cost = ComputePerfectMatchingCost(node_num, edge_num, edges, weights, pm); + printf("cost = %.1f\n", cost); + if (save_filename) SaveMatching(node_num, pm, save_filename); + delete pm; + } + else + { + int geom_node_num; + int* x_array; + int* y_array; + LoadGeomFile(geom_node_num, x_array, y_array, geom_filename); + GeomPerfectMatching *gpm = new GeomPerfectMatching(geom_node_num, 2); + for (i=0; iAddPoint(coord); + } + delete [] x_array; + delete [] y_array; + gpm->options = options; + gpm->gpm_options = gpm_options; + if (filename) + { + if (node_num != geom_node_num) { printf("%s and %s don't match!\n", geom_filename, filename); exit(1); } + for (e=0; eDist(edges[2*e], edges[2*e+1])) + { printf("edge lengths in %s and %s don't match!\n", geom_filename, filename); exit(1); } + gpm->AddInitialEdge(edges[2*e], edges[2*e+1]); + } + } + gpm->Solve(); + if (save_filename) SaveMatching(geom_node_num, gpm, save_filename); + delete gpm; + } + + if (filename) + { + delete [] edges; + delete [] weights; + } + + return 0; +} + + diff --git a/blossom5-v2.05.src/misc.cpp b/blossom5-v2.05.src/misc.cpp new file mode 100644 index 0000000..5bdfd3e --- /dev/null +++ b/blossom5-v2.05.src/misc.cpp @@ -0,0 +1,172 @@ +#include +#include "PerfectMatching.h" +#include "LCA.h" + +struct Node +{ + PerfectMatching::REAL sum; // = twice_y[i] + twice_y[i->parent] + twice_y[i->parent->parent] + ... + Node* match; + Node* parent; + Node* child; + Node* sibling; + int lca_preorder; +}; + + +int CheckPerfectMatchingOptimality(int node_num, int edge_num, int* edges, int* weights, PerfectMatching* pm, PerfectMatching::REAL threshold) +{ + int _i, _j, _e; + Node* i; + Node* j; + int blossom_num = pm->GetBlossomNum(); + int* blossom_parents = new int[node_num+blossom_num]; + PerfectMatching::REAL* twice_y = new PerfectMatching::REAL[node_num+blossom_num]; + + PerfectMatching::REAL y_blossom_min = 0; + PerfectMatching::REAL slack_min = 0; + PerfectMatching::REAL active_slack_max = 0; + + // step 1 - read dual solution and construct tree + pm->GetDualSolution(blossom_parents, twice_y); + Node* nodes = new Node[node_num+blossom_num+1]; + memset(nodes, 0, (node_num+blossom_num+1)*sizeof(Node)); + Node* ROOT = nodes+node_num+blossom_num; + for (_i=0, i=nodes; _isum = twice_y[_i]; + if (_i >= node_num && y_blossom_min > i->sum) y_blossom_min = i->sum; + if (blossom_parents[_i] >= 0) + { + if (blossom_parents[_i]=node_num+blossom_num) + { + delete [] nodes; + delete [] blossom_parents; + delete [] twice_y; + return 2; + } + i->parent = nodes + blossom_parents[_i]; + i->sibling = i->parent->child; + i->parent->child = i; + } + } + delete [] blossom_parents; + delete [] twice_y; + + for (i=nodes; iparent) + { + i->parent = ROOT; + i->sibling = ROOT->child; + ROOT->child = i; + } + } + + LCATree* lca_tree = new LCATree(node_num+blossom_num+1); + Node** rev_mapping = new Node*[node_num+blossom_num]; + + i = ROOT; + while ( 1 ) + { + if (i->child) + { + if (i < nodes+node_num) { delete [] nodes; delete lca_tree; delete [] rev_mapping; return 2; } + i->child->sum += i->sum; + i = i->child; + } + else + { + if (i >= nodes+node_num) { delete [] nodes; delete lca_tree; delete [] rev_mapping; return 2; } + while ( 1 ) + { + i->lca_preorder = lca_tree->Add(i, i->parent); + rev_mapping[i->lca_preorder] = i; + if (i->sibling) break; + i = i->parent; + if (i == ROOT) + { + i->lca_preorder = lca_tree->AddRoot(i); + break; + } + } + if (i == ROOT) break; + i = i->sibling; + i->sum += i->parent->sum; + } + } + + int matched_num = 0; + for (_e=0; _e=node_num || _j>=node_num || _i==_j) { delete [] nodes; delete lca_tree; delete [] rev_mapping; return 2; } + + int lca_i = nodes[_i].lca_preorder; + int lca_j = nodes[_j].lca_preorder; + lca_tree->GetPenultimateNodes(lca_i, lca_j); + i = rev_mapping[lca_i]; + j = rev_mapping[lca_j]; + PerfectMatching::REAL twice_slack = 2*weights[_e] - (nodes[_i].sum - i->parent->sum) - (nodes[_j].sum - j->parent->sum); + if (slack_min > twice_slack) slack_min = twice_slack; + if (pm->GetSolution(_e)) + { + if (pm->GetMatch(_i)!=_j || pm->GetMatch(_j)!=_i || i->match || j->match) { delete [] nodes; delete lca_tree; delete [] rev_mapping; return 2; } + i->match = j; + j->match = i; + if (active_slack_max < twice_slack) active_slack_max = twice_slack; + matched_num += 2; + } + } + + delete [] nodes; + delete lca_tree; + delete [] rev_mapping; + + if (matched_num != node_num) return 2; + + if (y_blossom_min < -threshold || slack_min < -threshold || active_slack_max > threshold) + { + printf("ERROR in CheckPerfectMatchingOptimality():\n"); + if ( ((PerfectMatching::REAL)1 / 2) == 0 ) + printf("\ty_blossom_min=%d\n\tslack_min=%d\n\tactive_slack_max=%d\n", (int)y_blossom_min, (int)slack_min, (int)active_slack_max); + else + printf("\ty_blossom_min=%.15f\n\tslack_min=%.15f\n\tactive_slack_max=%.15f\n", (double)y_blossom_min, (double)slack_min, (double)active_slack_max); + return 1; + } + + return 0; +} + +double ComputePerfectMatchingCost(int node_num, int edge_num, int* edges, int* weights, PerfectMatching* pm) +{ + int i; + int j; + int e; + double cost = 0; + + int* nodes = new int[node_num]; + memset(nodes, 0, node_num*sizeof(int)); + for (e=0; eGetSolution(e)) + { + i = edges[2*e]; + j = edges[2*e+1]; + nodes[i] ++; + nodes[j] ++; + cost += weights[e]; + } + } + for (i=0; i + + inline double get_time() + { + LARGE_INTEGER t, frequency; + QueryPerformanceCounter(&t); + QueryPerformanceFrequency(&frequency); + return (double)t.QuadPart/(double)frequency.QuadPart; + } + +#endif + +////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////// + +#ifdef PM_TIMER_CLOCK_GETTIME + + #include + + inline double get_time() + { + struct timespec t; + clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &t); + return (double)t.tv_nsec*1.00E-9 + (double)t.tv_sec; + } + +#endif + +////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////// + +#ifdef PM_TIMER_GETRUSAGE + + #include + + inline double get_time() + { + struct rusage t; + getrusage (RUSAGE_SELF, &t); + return (double)t.ru_utime.tv_usec*1.00E-6 + (double)t.ru_utime.tv_sec; + } + +#endif + +////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////// + +#ifdef PM_TIMER_EXTERNAL + + extern double get_time(); + +#endif + +////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////// + +#ifdef PM_TIMER_NONE + + inline double get_time() { return 0; } + +#endif + +#endif + diff --git a/hungarian.cpp b/hungarian.cpp index 00cc91d..b8cb7ab 100644 --- a/hungarian.cpp +++ b/hungarian.cpp @@ -75,7 +75,7 @@ public: } } bool isTight() const { - return std::abs(s->y + t->y - weight) < 1e-10; + return std::abs(s->y + t->y - weight) < 1e-13; } }; @@ -166,10 +166,10 @@ public: if (e.isTight() && (e.head().index == i && !e.inPath)) { if (!e.tail().inMatching()) { break; + } else { + e.inPath= true; + path.push({e.tail().neighbors.begin(), e.tail().neighbors.end(), e.tail().index}); } - - e.inPath= true; - path.push({e.tail().neighbors.begin(), e.tail().neighbors.end(), e.tail().index}); } else { it++; path.top() = {it, itEnd, i}; @@ -177,12 +177,6 @@ public: } } - - if (path.empty()) { - exit(1); - } - - while (!path.empty()) { auto [it, itEnd, i] = path.top(); Edge& e = it->get(); diff --git a/rbmp.cpp b/rbmp.cpp index ffc7c44..ddd76cb 100644 --- a/rbmp.cpp +++ b/rbmp.cpp @@ -8,6 +8,8 @@ #include "randutils/randutils.hpp" #include "pcg-cpp/include/pcg_random.hpp" +#include "blossom5-v2.05.src/PerfectMatching.h" + using Rng = randutils::random_generator; class Edge; @@ -31,9 +33,8 @@ private: public: const Edge& edge; - double X; - HalfEdge(const Edge& e) : edge(e), X(0) {} + HalfEdge(const Edge& e) : edge(e) {} void setVertices(Vertex& vt, Vertex& vh) { tail = &vt; head = &vh; @@ -58,12 +59,6 @@ public: red.addEdge(halfedges[0]); blue.addEdge(halfedges[1]); } - double belief() const { - return halfedges[0].X + halfedges[1].X; - } - bool active() const { - return belief() >= 0; - } }; class Graph { @@ -86,51 +81,18 @@ public: edges[i].weight = r.variate(1); } } - - void propagateBeliefs() { - for (unsigned i : {0, 1}) { -#pragma omp parallel for - for (Edge& e : edges) { - HalfEdge& h = e.halfedges[i]; - h.X = std::numeric_limits::infinity(); - for (const HalfEdge& hn : h.getHead().neighbors) { - if (h.getTail().index != hn.getHead().index) { - h.X = std::min(hn.edge.weight - hn.X, h.X); - } - } - } - } - } - - double minBelief() const { - double minBelief = std::numeric_limits::infinity(); - - for (const Edge& e : edges) { - minBelief = std::min(std::abs(e.belief()), minBelief); - } - - return minBelief; - } }; int main(int argc, char* argv[]) { unsigned n = 100; - unsigned maxSteps = 1e8; - double beliefThreshold = 1; int opt; - while ((opt = getopt(argc, argv, "n:m:t:")) != -1) { + while ((opt = getopt(argc, argv, "n:")) != -1) { switch (opt) { case 'n': n = atoi(optarg); break; - case 'm': - maxSteps = (unsigned)atof(optarg); - break; - case 't': - beliefThreshold = atof(optarg); - break; default: exit(1); } @@ -139,26 +101,60 @@ int main(int argc, char* argv[]) { Rng r; Graph G(n, r); - unsigned steps = 0; - double minBelief = 0; + PerfectMatching pm(G.vertices.size(), G.edges.size()); - while (minBelief < beliefThreshold && steps < maxSteps) { - G.propagateBeliefs(); - if (steps % 100 == 0) { - std::cerr << (minBelief = G.minBelief()) << std::endl; - } - steps++; + for (const Edge& e : G.edges) { + pm.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, e.weight); } + pm.options.verbose = false; + pm.Solve(); + + std::reference_wrapper eFlip = r.pick(G.edges); + + while (pm.GetMatch(eFlip.get().halfedges[0].getTail().index) != eFlip.get().halfedges[0].getHead().index) { + eFlip = r.pick(G.edges); + } + + eFlip.get().weight = 1e10; + + PerfectMatching pm2(G.vertices.size(), G.edges.size()); + for (const Edge& e : G.edges) { - if (e.active()) { + pm2.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, e.weight); + } + + pm2.options.verbose = false; + pm2.Solve(); + + std::cout << n << std::endl; + + for (unsigned i = 0; i < G.vertices.size() / 2; i++) { + unsigned j1 = pm.GetMatch(i); + + std::cout + << G.vertices[i].coordinate[0] << " " + << G.vertices[i].coordinate[1] << " " + << G.vertices[j1].coordinate[0] << " " + << G.vertices[j1].coordinate[1] << std::endl; + } + + unsigned m = 0; + for (unsigned i = 0; i < G.vertices.size() / 2; i++) { + unsigned j1 = pm.GetMatch(i); + unsigned j2 = pm2.GetMatch(i); + + if (j1 != j2) { + m++; std::cout - << e.halfedges[0].getTail().coordinate[0] << " " - << e.halfedges[0].getTail().coordinate[1] << " " - << e.halfedges[0].getHead().coordinate[0] << " " - << e.halfedges[0].getHead().coordinate[1] << std::endl; + << G.vertices[i].coordinate[0] << " " + << G.vertices[i].coordinate[1] << " " + << G.vertices[j2].coordinate[0] << " " + << G.vertices[j2].coordinate[1] << std::endl; } } + std::cerr << "Size of excitation: " << m << std::endl; + return 0; } -- cgit v1.2.3-70-g09d2