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/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 ++ 6 files changed, 1141 insertions(+) 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 (limited to 'blossom5-v2.05.src/MinCost') 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 -- cgit v1.2.3-70-g09d2