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