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. --- blossom5-v2.05.src/GEOM/GeomPerfectMatching.h | 177 ++++++++++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 blossom5-v2.05.src/GEOM/GeomPerfectMatching.h (limited to 'blossom5-v2.05.src/GEOM/GeomPerfectMatching.h') 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 -- cgit v1.2.3-54-g00ecf