MRPT  1.9.9
CGraphPartitioner.cpp
Go to the documentation of this file.
1 /* +------------------------------------------------------------------------+
2  | Mobile Robot Programming Toolkit (MRPT) |
3  | https://www.mrpt.org/ |
4  | |
5  | Copyright (c) 2005-2019, Individual contributors, see AUTHORS file |
6  | See: https://www.mrpt.org/Authors - All rights reserved. |
7  | Released under BSD License. See: https://www.mrpt.org/License |
8  +------------------------------------------------------------------------+ */
9 
10 #include "graphs-precomp.h" // Precompiled headers
11 
13 #include <mrpt/math/CMatrixD.h>
14 #include <mrpt/math/CMatrixF.h>
15 #include <mrpt/math/ops_matrices.h> // laplacian()
16 #include <iostream> // cout
17 
18 using namespace mrpt;
19 using namespace mrpt::graphs;
20 using namespace std;
21 
22 /*---------------------------------------------------------------
23 SpectralPartition
24 ---------------------------------------------------------------*/
25 template <class GRAPH_MATRIX, typename num_t>
27  GRAPH_MATRIX& in_A, std::vector<uint32_t>& out_part1,
28  std::vector<uint32_t>& out_part2, num_t& out_cut_value, bool forceSimetry)
29 {
30  size_t nodeCount; // Nodes count
31  GRAPH_MATRIX Adj, eigenVectors;
32  std::vector<typename GRAPH_MATRIX::Scalar> eigenValues;
33 
34  // Check matrix is square:
35  if (in_A.cols() != int(nodeCount = in_A.rows()))
36  THROW_EXCEPTION("Weights matrix is not square!!");
37 
38  // Shi & Malik's method
39  //-----------------------------------------------------------------
40 
41  // forceSimetry?
42  if (forceSimetry)
43  {
44  Adj.setSize(nodeCount, nodeCount);
45  for (size_t i = 0; i < nodeCount; i++)
46  for (size_t j = i; j < nodeCount; j++)
47  Adj(i, j) = Adj(j, i) = 0.5f * (in_A(i, j) + in_A(j, i));
48  }
49  else
50  Adj = in_A;
51 
52  // Compute eigen-vectors of laplacian:
53  GRAPH_MATRIX LAPLACIAN;
54  mrpt::math::laplacian(Adj, LAPLACIAN);
55 
56  LAPLACIAN.eig(eigenVectors, eigenValues);
57 
58  // Execute the bisection
59  // ------------------------------------
60  // Second smallest eigen-vector
61  double mean = 0;
62  size_t colNo = 1; // second smallest
63  size_t nRows = eigenVectors.rows();
64 
65  // for (i=0;i<eigenVectors.cols();i++) mean+=eigenVectors(colNo,i);
66  for (size_t i = 0; i < nRows; i++) mean += eigenVectors(i, colNo);
67  mean /= nRows;
68 
69  out_part1.clear();
70  out_part2.clear();
71 
72  for (size_t i = 0; i < nRows; i++)
73  {
74  if (eigenVectors(i, colNo) >= mean)
75  out_part1.push_back(i);
76  else
77  out_part2.push_back(i);
78  }
79 
80  // Special and strange case: Constant eigenvector: Split nodes in two
81  // equally sized parts arbitrarily:
82  // ------------------------------------------------------------------
83  if (!out_part1.size() || !out_part2.size())
84  {
85  out_part1.clear();
86  out_part2.clear();
87  // Assign 50%-50%:
88  for (int i = 0; i < Adj.cols(); i++)
89  if (i <= Adj.cols() / 2)
90  out_part1.push_back(i);
91  else
92  out_part2.push_back(i);
93  }
94 
95  // Compute the N-cut value
96  out_cut_value = nCut(Adj, out_part1, out_part2);
97 }
98 
99 /*---------------------------------------------------------------
100 RecursiveSpectralPartition
101 ---------------------------------------------------------------*/
102 template <class GRAPH_MATRIX, typename num_t>
104  GRAPH_MATRIX& in_A, std::vector<std::vector<uint32_t>>& out_parts,
105  num_t threshold_Ncut, bool forceSimetry, bool useSpectralBisection,
106  bool recursive, unsigned minSizeClusters, const bool verbose)
107 {
108  MRPT_START
109 
110  size_t nodeCount;
111  std::vector<uint32_t> p1, p2;
112  num_t cut_value;
113  size_t i, j;
114  GRAPH_MATRIX Adj;
115 
116  out_parts.clear();
117 
118  // Check matrix is square:
119  if (in_A.cols() != int(nodeCount = in_A.rows()))
120  THROW_EXCEPTION("Weights matrix is not square!!");
121 
122  if (nodeCount == 1)
123  {
124  // Don't split, there is just a node!
125  p1.push_back(0);
126  out_parts.push_back(p1);
127  return;
128  }
129 
130  // forceSimetry?
131  if (forceSimetry)
132  {
133  Adj.setSize(nodeCount, nodeCount);
134  for (i = 0; i < nodeCount; i++)
135  for (j = i; j < nodeCount; j++)
136  Adj(i, j) = Adj(j, i) = 0.5f * (in_A(i, j) + in_A(j, i));
137  }
138  else
139  Adj = in_A;
140 
141  // Make bisection
142  if (useSpectralBisection)
143  SpectralBisection(Adj, p1, p2, cut_value, false);
144  else
145  exactBisection(Adj, p1, p2, cut_value, false);
146 
147  if (verbose)
148  std::cout << format(
149  "Cut:%u=%u+%u,nCut=%.02f->", (unsigned int)nodeCount,
150  (unsigned int)p1.size(), (unsigned int)p2.size(), cut_value);
151 
152  // Is it a useful partition?
153  if (cut_value > threshold_Ncut || p1.size() < minSizeClusters ||
154  p2.size() < minSizeClusters)
155  {
156  if (verbose) std::cout << "->NO!" << std::endl;
157 
158  // No:
159  p1.clear();
160  for (i = 0; i < nodeCount; i++) p1.push_back(i);
161  out_parts.push_back(p1);
162  }
163  else
164  {
165  if (verbose) std::cout << "->YES!" << std::endl;
166 
167  // Yes:
168  std::vector<std::vector<uint32_t>> p1_parts, p2_parts;
169 
170  if (recursive)
171  {
172  // Split "p1":
173  // --------------------------------------------
174  // sub-matrix:
175  GRAPH_MATRIX A_1(p1.size(), p1.size());
176  for (i = 0; i < p1.size(); i++)
177  for (j = 0; j < p1.size(); j++) A_1(i, j) = in_A(p1[i], p1[j]);
178 
179  RecursiveSpectralPartition(
180  A_1, p1_parts, threshold_Ncut, forceSimetry,
181  useSpectralBisection, recursive, minSizeClusters);
182 
183  // Split "p2":
184  // --------------------------------------------
185  // sub-matrix:
186  GRAPH_MATRIX A_2(p2.size(), p2.size());
187  for (i = 0; i < p2.size(); i++)
188  for (j = 0; j < p2.size(); j++) A_2(i, j) = in_A(p2[i], p2[j]);
189 
190  RecursiveSpectralPartition(
191  A_2, p2_parts, threshold_Ncut, forceSimetry,
192  useSpectralBisection, recursive, minSizeClusters);
193 
194  // Build "out_parts" from "p1_parts" + "p2_parts"
195  // taken care of indexes mapping!
196  // -------------------------------------------------------------------------------------
197  // remap p1 nodes:
198  for (i = 0; i < p1_parts.size(); i++)
199  {
200  for (j = 0; j < p1_parts[i].size(); j++)
201  p1_parts[i][j] = p1[p1_parts[i][j]];
202  out_parts.push_back(p1_parts[i]);
203  }
204 
205  // remap p2 nodes:
206  for (i = 0; i < p2_parts.size(); i++)
207  {
208  for (j = 0; j < p2_parts[i].size(); j++)
209  p2_parts[i][j] = p2[p2_parts[i][j]];
210 
211  out_parts.push_back(p2_parts[i]);
212  }
213  }
214  else
215  {
216  // Force bisection only:
217  out_parts.clear();
218  out_parts.push_back(p1);
219  out_parts.push_back(p2);
220  }
221  }
222 
223  MRPT_END
224 }
225 
226 /*---------------------------------------------------------------
227 nCut
228 ---------------------------------------------------------------*/
229 template <class GRAPH_MATRIX, typename num_t>
231  const GRAPH_MATRIX& in_A, const std::vector<uint32_t>& in_part1,
232  const std::vector<uint32_t>& in_part2)
233 {
234  size_t i, j;
235  const size_t size1 = in_part1.size();
236  const size_t size2 = in_part2.size();
237 
238  // Compute the N-cut value
239  // -----------------------------------------------
240  num_t cut_AB = 0;
241  for (i = 0; i < size1; i++)
242  for (j = 0; j < size2; j++) cut_AB += in_A(in_part1[i], in_part2[j]);
243 
244  num_t assoc_AA = 0;
245  for (i = 0; i < size1; i++)
246  for (j = i; j < size1; j++)
247  if (i != j) assoc_AA += in_A(in_part1[i], in_part1[j]);
248 
249  num_t assoc_BB = 0;
250  for (i = 0; i < size2; i++)
251  for (j = i; j < size2; j++)
252  if (i != j) assoc_BB += in_A(in_part2[i], in_part2[j]);
253 
254  num_t assoc_AV = assoc_AA + cut_AB;
255  num_t assoc_BV = assoc_BB + cut_AB;
256 
257  if (!cut_AB)
258  return 0;
259  else
260  return cut_AB / assoc_AV + cut_AB / assoc_BV;
261 }
262 
263 /*---------------------------------------------------------------
264 exactBisection
265 ---------------------------------------------------------------*/
266 template <class GRAPH_MATRIX, typename num_t>
268  GRAPH_MATRIX& in_A, std::vector<uint32_t>& out_part1,
269  std::vector<uint32_t>& out_part2, num_t& out_cut_value, bool forceSimetry)
270 {
271  size_t nodeCount; // Nodes count
272  size_t i, j;
273  GRAPH_MATRIX Adj;
274  std::vector<bool> partition, bestPartition;
275  std::vector<uint32_t> part1, part2;
276  num_t partCutValue, bestCutValue = std::numeric_limits<num_t>::max();
277  bool end = false;
278 
279  // Check matrix is square:
280  if (in_A.cols() != int(nodeCount = in_A.rows()))
281  THROW_EXCEPTION("Weights matrix is not square!!");
282 
283  ASSERT_(nodeCount >= 2);
284 
285  // forceSimetry?
286  if (forceSimetry)
287  {
288  Adj.setSize(nodeCount, nodeCount);
289  for (i = 0; i < nodeCount; i++)
290  for (j = i; j < nodeCount; j++)
291  Adj(i, j) = Adj(j, i) = 0.5f * (in_A(i, j) + in_A(j, i));
292  }
293  else
294  Adj = in_A;
295 
296  // Brute force: compute all posible partitions:
297  //-----------------------------------------------------------------
298 
299  // First combination: 1000...0000
300  // Last combination: 1111...1110
301  partition.clear();
302  partition.resize(nodeCount, false);
303  partition[0] = true;
304 
305  while (!end)
306  {
307  // Build partitions from binary vector:
308  part1.clear();
309  part2.clear();
310 
311  for (i = 0; i < nodeCount; i++)
312  {
313  if (partition[i])
314  part2.push_back(i);
315  else
316  part1.push_back(i);
317  }
318 
319  // Compute the n-cut:
320  partCutValue = nCut(Adj, part1, part2);
321 
322  if (partCutValue < bestCutValue)
323  {
324  bestCutValue = partCutValue;
325  bestPartition = partition;
326  }
327 
328  // Next combo in the binary vector:
329  i = 0;
330  bool carry = false;
331  do
332  {
333  carry = partition[i];
334  partition[i] = !partition[i];
335  i++;
336  } while (carry && i < nodeCount);
337 
338  // End criterion:
339  end = true;
340  for (i = 0; end && i < nodeCount; i++) end = end && partition[i];
341 
342  }; // End of while
343 
344  // Return the best partition:
345  out_cut_value = bestCutValue;
346 
347  out_part1.clear();
348  out_part2.clear();
349 
350  for (i = 0; i < nodeCount; i++)
351  {
352  if (bestPartition[i])
353  out_part2.push_back(i);
354  else
355  out_part1.push_back(i);
356  }
357 }
358 
359 // Explicit instantiations:
360 namespace mrpt::graphs
361 {
362 // Explicit instantiations:
367 } // namespace mrpt::graphs
#define MRPT_START
Definition: exceptions.h:241
static num_t nCut(const GRAPH_MATRIX &in_A, const std::vector< uint32_t > &in_part1, const std::vector< uint32_t > &in_part2)
Returns the normaliced cut of a graph, given its adjacency matrix A and a bisection: ...
#define THROW_EXCEPTION(msg)
Definition: exceptions.h:67
std::string std::string format(std::string_view fmt, ARGS &&... args)
Definition: format.h:26
Abstract graph and tree data structures, plus generic graph algorithms.
void laplacian(const MATIN &g, MATOUT &ret)
Computes the Laplacian of a square graph weight matrix.
Definition: ops_matrices.h:207
This file implements miscelaneous matrix and matrix/vector operations, and internal functions in mrpt...
STL namespace.
static void exactBisection(GRAPH_MATRIX &in_A, std::vector< uint32_t > &out_part1, std::vector< uint32_t > &out_part2, num_t &out_cut_value, bool forceSimetry=true)
Performs an EXACT minimum n-Cut graph bisection, (Use CGraphPartitioner::SpectralBisection for a fast...
#define ASSERT_(f)
Defines an assertion mechanism.
Definition: exceptions.h:120
const_iterator end() const
Definition: ts_hash_map.h:246
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
static void RecursiveSpectralPartition(GRAPH_MATRIX &in_A, std::vector< std::vector< uint32_t >> &out_parts, num_t threshold_Ncut=1, bool forceSimetry=true, bool useSpectralBisection=true, bool recursive=true, unsigned minSizeClusters=1, const bool verbose=false)
Performs the spectral recursive partition into K-parts for a given graph.
#define MRPT_END
Definition: exceptions.h:245
double mean(const CONTAINER &v)
Computes the mean value of a vector.
static void SpectralBisection(GRAPH_MATRIX &in_A, std::vector< uint32_t > &out_part1, std::vector< uint32_t > &out_part2, num_t &out_cut_value, bool forceSimetry=true)
Performs the spectral bisection of a graph.
Finds the min-normalized-cut of a weighted undirected graph.



Page generated by Doxygen 1.8.14 for MRPT 1.9.9 Git: 7e629e01a Sat Dec 14 00:05:55 2019 +0100 at sáb dic 14 00:15:10 CET 2019