TTK
Loading...
Searching...
No Matches
GabowTarjan.cpp
Go to the documentation of this file.
1#include <GabowTarjan.h>
2
3#include <iostream>
4#include <queue>
5#include <vector>
6
7bool ttk::GabowTarjan::DFS(int v) {
8 if(v < 0)
9 return true;
10
11 // for every adjacent vertex u of v
12 for(unsigned int i = 0; i < Connections[v].size(); ++i) {
13 int u = Connections[v][i];
14 if(Layers[Pair[u] + 1] == Layers[v + 1] + 1) {
15 if(DFS(Pair[u])) {
16 Pair[u] = v;
17 Pair[v] = u;
18 return true;
19 }
20 }
21 }
22
23 Layers[v + 1] = -1;
24
25 return false;
26}
27
28bool ttk::GabowTarjan::BFS() {
29 std::queue<int> vertexQueue;
30
31 // For every vertex v given by PersistencePairs1
32 for(unsigned int v = 0; v < MaxSize; v++) {
33 // If its not paired to vertex in PersistencePairs2
34 if(Pair[v] < 0) {
35 // Set its layer to 0 and put it in the queue
36 Layers[v + 1] = 0;
37 vertexQueue.push(v);
38 } else {
39 // Otherwise mark it as matched (set Layer to NIL)
40 Layers[v + 1] = -1;
41 }
42 }
43
44 // Set layer for NIL
45 Layers[0] = -1;
46
47 // Search the vertices in the queue
48 while(!vertexQueue.empty()) {
49 int v = vertexQueue.front();
50 vertexQueue.pop();
51 if(Layers[v + 1] > Layers[0]) {
52 for(unsigned int i = 0; i < Connections[v].size(); ++i) {
53 int u = Connections[v][i];
54 // Check if the vertex has an edge to the match vertex
55 if(Layers[Pair[u] + 1] < 0) {
56 // Set the layer of the vertex (it can be NIL) which is matched to the
57 // matched vertex u
58 Layers[Pair[u] + 1] = Layers[v + 1] + 1;
59 // If the pairing vertex is not NIL add it into the queue
60 if(Pair[u] != -1)
61 vertexQueue.push(Pair[u]);
62 }
63 }
64 }
65 }
66
67 return Layers[0] != -1;
68}
69
70void ttk::GabowTarjan::HopcroftKarp(unsigned int &matching) {
71 while(BFS())
72 for(unsigned int vertex = 0; vertex < MaxSize; ++vertex) {
73 if(Pair[vertex] == -1) {
74 if(DFS(vertex))
75 ++matching;
76 }
77 }
78}
79
81 // Clear the pairing
82 Pair.clear();
83 Pair.assign(2 * MaxSize, -1);
84
85 // Clearing Layers
86 Layers.clear();
87 Layers.resize(MaxSize + 1);
88
89 // No vertices are matched
90 unsigned int matching = 0;
91 unsigned int matching0 = 0;
92
93 // Clear the connection matrix and set it to the right size
94 Connections.clear();
95 Connections.resize(MaxSize);
96
97 // The maximal weight of the edges which are used for the matching
98 double currentWeight = 0;
99
100 // First non added edge is an iterator pointing to the first edge
101 // in Edges which was added to the Connections
102 unsigned int firstNotAddedEdge = 0;
103 unsigned int nbEdges = (unsigned int)Edges.size();
104
105 unsigned int lowerBound = 0;
106 unsigned int upperBound = nbEdges;
107 unsigned int guessEdge = (lowerBound + nbEdges) / 2;
108
109 std::map<int, int> offPairings;
110
111 // Repeat till all the vertices are matched
112 while(matching < MaxSize) {
113 // Save initial matching.
114 matching0 = matching;
115 unsigned int oldGuessEdge = guessEdge;
116
117 currentWeight = Edges[guessEdge].weight;
118 while(Edges[guessEdge].weight == currentWeight && guessEdge < nbEdges)
119 ++guessEdge;
120
121 if(guessEdge >= nbEdges) {
122 this->printMsg("ran out of edges.");
123 }
124
125 // Add the edges with the current weight (distance) to the connection matrix
126 while(firstNotAddedEdge >= guessEdge) {
127 int v1 = Edges[firstNotAddedEdge].v1;
128 int v2 = Edges[firstNotAddedEdge].v2;
129
130 std::vector<int> vec1 = Connections[v1];
131 vec1.erase(std::remove(vec1.begin(), vec1.end(), v2), vec1.end());
132 Connections[v1] = vec1;
133
134 --firstNotAddedEdge;
135 }
136
137 // Edges[firstNotAddedEdge].weight == currentWeight && firstNotAddedEdge <
138 // nbEdges
139 while(firstNotAddedEdge < guessEdge) {
140 int v1 = Edges[firstNotAddedEdge].v1;
141 int v2 = Edges[firstNotAddedEdge].v2;
142
143 Connections[v1].push_back(v2);
144 ++firstNotAddedEdge;
145 }
146
147 // Clear the pairing
148 Pair.clear();
149 Pair.assign(2 * MaxSize, -1);
150 // Clearing Layers
151 Layers.clear();
152 Layers.resize(MaxSize + 1);
153
154 // Temporarily augment matching.
155 this->printMsg("Guessing for " + std::to_string(guessEdge) + "...");
156
157 matching = 0;
158 HopcroftKarp(matching);
159 // printCurrentMatching<dataType>();
160
161 if(matching >= MaxSize) {
162
163 // Reset matching.
164 matching = matching0;
165 upperBound = guessEdge;
166 guessEdge = (lowerBound + guessEdge) / 2;
167
168 // Ended binary search.
169 if(oldGuessEdge == guessEdge) {
170 this->printMsg("Binary search success.");
171 return Edges[(guessEdge > 0 ? guessEdge - 1 : guessEdge)].weight;
172 }
173
174 } else {
175
176 // Check if we did not run out of edges. This should never happen.
177 if(firstNotAddedEdge == nbEdges) {
178 std::stringstream msg;
179 ttk::Debug d;
180 this->printMsg("Not enough edges to find the matching!");
181 // return -1;
182 return currentWeight;
183 }
184
185 lowerBound = guessEdge;
186 guessEdge = (guessEdge + upperBound) / 2;
187 // Increase the value of the current weight
188 // currentWeight = Edges[firstNotAddedEdge].weight;
189
190 // Ended binary search.
191 if(oldGuessEdge == guessEdge) {
192 this->printMsg("Binary search success.");
193 return Edges[(guessEdge > 0 ? guessEdge - 1 : guessEdge)].weight;
194 }
195 }
196 }
197
198 return -1;
199}
200
202 int size = 2 * MaxSize;
203 std::vector<int> missedPlaces;
204
205 {
206 std::stringstream msg;
207 msg << "Assignment matrix: " << std::endl;
208 for(int i = 0; i < size; ++i) {
209 int k = Pair[i];
210 if(k < 0 || k > size)
211 missedPlaces.push_back(i);
212 for(int j = 0; j < size; ++j) {
213 msg << (j == k ? "1 " : "0 ");
214 }
215 msg << std::endl;
216 }
217 msg << "/Assignment matrix." << std::endl << std::endl;
218 this->printMsg(msg.str(), debug::Priority::VERBOSE);
219 }
220
221 {
222 std::stringstream msg;
223 ttk::Debug d;
224 msg << "Missed:" << std::endl;
225 for(unsigned int i = 0; i < missedPlaces.size(); ++i) {
226 msg << missedPlaces.at(i) << " ";
227 }
228 msg << std::endl << std::endl;
229 this->printMsg(msg.str());
230 }
231}
232
233int ttk::GabowTarjan::run(std::vector<MatchingType> &matchings) {
234 // Compute distance.
235 double dist = this->Distance();
236 this->printMsg("Computed distance " + std::to_string(dist));
237
238 // Fill matchings.
239 matchings.clear();
240
241 for(unsigned int i = 0; i < Size1; ++i) {
242 if(Pair[i] == -1)
243 continue;
244
245 int j = Pair[i] - MaxSize;
246 if(j <= -1 || (j < (int)Size2 && Pair[j + MaxSize] != (int)i)) {
247 this->printErr("Hopcroft-Karp built an invalid matching.");
248 // return -1;
249 }
250
251 if(j >= (int)Size2) {
252 matchings.emplace_back(i, j, (*Cptr)[i][Size2]);
253 } else {
254 matchings.emplace_back(i, j, (*Cptr)[i][j]);
255 }
256 }
257
258 for(unsigned int j = Size1; j < MaxSize; ++j) {
259 if(Pair[j] == -1)
260 continue;
261
262 int i = Pair[j] - MaxSize - Size2;
263 if(i > -1 && (i >= (int)Size1 || Pair[i + MaxSize + Size2] != (int)j)) {
264 this->printErr("Hopcroft-Karp built an invalid matching.");
265 // return -1;
266 }
267
268 if(i <= -1) {
269 matchings.emplace_back(i, j - Size1, (*Cptr)[Size1][j - Size1]);
270 } else {
271 // Already added by symmetry.
272 matchings.emplace_back(i, j - Size1, (*Cptr)[i][j - Size1]);
273 }
274 }
275
276 return 0;
277}
Minimalist debugging class.
Definition: Debug.h:88
int run(std::vector< MatchingType > &matchings)
void printCurrentMatching()
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)