Main MRPT website > C++ reference for MRPT 1.5.7
CGraphPartitioner_impl.h
Go to the documentation of this file.
1 /* +---------------------------------------------------------------------------+
2  | Mobile Robot Programming Toolkit (MRPT) |
3  | http://www.mrpt.org/ |
4  | |
5  | Copyright (c) 2005-2017, Individual contributors, see AUTHORS file |
6  | See: http://www.mrpt.org/Authors - All rights reserved. |
7  | Released under BSD License. See details in http://www.mrpt.org/License |
8  +---------------------------------------------------------------------------+ */
9 
10 
11 #if !defined(CGRAPHPARTITIONER_H)
12 #error "This file can't be included from outside of CGraphPartitioner.h"
13 #endif
14 
15 namespace mrpt
16 {
17 namespace graphs
18 {
19 
20 /*---------------------------------------------------------------
21  SpectralPartition
22  ---------------------------------------------------------------*/
23 template <class GRAPH_MATRIX, typename num_t>
24 void CGraphPartitioner<GRAPH_MATRIX,num_t>::SpectralBisection(
25  GRAPH_MATRIX &in_A,
26  vector_uint &out_part1,
27  vector_uint &out_part2,
28  num_t &out_cut_value,
29  bool forceSimetry )
30 {
31  size_t nodeCount; // Nodes count
32  GRAPH_MATRIX Adj, eigenVectors,eigenValues;
33 
34  // Check matrix is square:
35  if (in_A.getColCount() != (nodeCount = in_A.getRowCount()) )
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 Adj = in_A;
50 
51  // Compute eigen-vectors of laplacian:
52  GRAPH_MATRIX LAPLACIAN;
53  Adj.laplacian(LAPLACIAN);
54 
55 
56  LAPLACIAN.eigenVectors( 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.getRowCount();
64 
65  //for (i=0;i<eigenVectors.getColCount();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 out_part2.push_back( i );
77  }
78 
79  // Special and strange case: Constant eigenvector: Split nodes in two
80  // equally sized parts arbitrarily:
81  // ------------------------------------------------------------------
82  if (!out_part1.size() || !out_part2.size())
83  {
84  out_part1.clear();
85  out_part2.clear();
86  // Assign 50%-50%:
87  for (size_t i=0;i<Adj.getColCount();i++)
88  if (i<=Adj.getColCount()/2)
89  out_part1.push_back(i);
90  else out_part2.push_back(i);
91  }
92 
93  // Compute the N-cut value
94  out_cut_value = nCut( Adj, out_part1, out_part2);
95 }
96 
97 /*---------------------------------------------------------------
98  RecursiveSpectralPartition
99  ---------------------------------------------------------------*/
100 template <class GRAPH_MATRIX, typename num_t>
101 void CGraphPartitioner<GRAPH_MATRIX,num_t>::RecursiveSpectralPartition(
102  GRAPH_MATRIX &in_A,
103  std::vector<vector_uint> &out_parts,
104  num_t threshold_Ncut,
105  bool forceSimetry,
106  bool useSpectralBisection,
107  bool recursive,
108  unsigned minSizeClusters,
109  const bool verbose )
110 {
111  MRPT_START
112 
113  size_t nodeCount;
114  vector_uint p1,p2;
115  num_t cut_value;
116  size_t i,j;
117  GRAPH_MATRIX Adj;
118 
119  out_parts.clear();
120 
121  // Check matrix is square:
122  if (in_A.getColCount() != (nodeCount = in_A.getRowCount()) )
123  THROW_EXCEPTION("Weights matrix is not square!!");
124 
125  if (nodeCount==1)
126  {
127  // Don't split, there is just a node!
128  p1.push_back(0);
129  out_parts.push_back(p1);
130  return;
131  }
132 
133  // forceSimetry?
134  if (forceSimetry)
135  {
136  Adj.setSize(nodeCount,nodeCount);
137  for (i=0;i<nodeCount;i++)
138  for (j=i;j<nodeCount;j++)
139  Adj(i,j)=Adj(j,i)= 0.5f*(in_A(i,j)+in_A(j,i));
140  }
141  else Adj = in_A;
142 
143  // Make bisection
144  if (useSpectralBisection)
145  SpectralBisection( Adj, p1, p2, cut_value, false);
146  else exactBisection(Adj, p1, p2, cut_value, false);
147 
148  if (verbose)
149  std::cout << format("Cut:%u=%u+%u,nCut=%.02f->",(unsigned int)nodeCount,(unsigned int)p1.size(),(unsigned int)p2.size(),cut_value);
150 
151  // Is it a useful partition?
152  if (cut_value>threshold_Ncut || p1.size()<minSizeClusters || p2.size()<minSizeClusters )
153  {
154  if (verbose)
155  std::cout << "->NO!" << std::endl;
156 
157  // No:
158  p1.clear();
159  for (i=0;i<nodeCount;i++) p1.push_back(i);
160  out_parts.push_back(p1);
161  }
162  else
163  {
164  if (verbose)
165  std::cout << "->YES!" << std::endl;
166 
167  // Yes:
168  std::vector<vector_uint> 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++)
178  A_1(i,j)= in_A(p1[i],p1[j]);
179 
180  RecursiveSpectralPartition(A_1,p1_parts, threshold_Ncut, forceSimetry,useSpectralBisection, recursive, minSizeClusters);
181 
182  // Split "p2":
183  // --------------------------------------------
184  // sub-matrix:
185  GRAPH_MATRIX A_2( p2.size(),p2.size() );
186  for (i=0;i<p2.size();i++)
187  for (j=0;j<p2.size();j++)
188  A_2(i,j)= in_A(p2[i],p2[j]);
189 
190  RecursiveSpectralPartition(A_2,p2_parts, threshold_Ncut, forceSimetry,useSpectralBisection, recursive, minSizeClusters);
191 
192  // Build "out_parts" from "p1_parts" + "p2_parts"
193  // taken care of indexes mapping!
194  // -------------------------------------------------------------------------------------
195  // remap p1 nodes:
196  for (i=0;i<p1_parts.size();i++)
197  {
198  for (j=0;j<p1_parts[i].size();j++)
199  p1_parts[i][j] = p1[ p1_parts[i][j] ];
200  out_parts.push_back(p1_parts[i]);
201  }
202 
203  // remap p2 nodes:
204  for (i=0;i<p2_parts.size();i++)
205  {
206  for (j=0;j<p2_parts[i].size();j++)
207  p2_parts[i][j] = p2[ p2_parts[i][j] ];
208 
209  out_parts.push_back(p2_parts[i]);
210  }
211  }
212  else
213  {
214  // Force bisection only:
215  out_parts.clear();
216  out_parts.push_back(p1);
217  out_parts.push_back(p2);
218  }
219 
220  }
221 
222  MRPT_END
223 }
224 
225 /*---------------------------------------------------------------
226  nCut
227  ---------------------------------------------------------------*/
228 template <class GRAPH_MATRIX, typename num_t>
229 num_t CGraphPartitioner<GRAPH_MATRIX,num_t>::nCut(
230  const GRAPH_MATRIX &in_A,
231  const vector_uint &in_part1,
232  const vector_uint &in_part2)
233 {
234  unsigned int i,j;
235  size_t size1=in_part1.size();
236  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++)
243  cut_AB += in_A(in_part1[i],in_part2[j]);
244 
245  num_t assoc_AA = 0;
246  for(i=0;i<size1;i++)
247  for(j=i;j<size1;j++)
248  if (i!=j)
249  assoc_AA += in_A(in_part1[i],in_part1[j]);
250 
251  num_t assoc_BB = 0;
252  for(i=0;i<size2;i++)
253  for(j=i;j<size2;j++)
254  if (i!=j)
255  assoc_BB += in_A(in_part2[i],in_part2[j]);
256 
257  num_t assoc_AV = assoc_AA + cut_AB;
258  num_t assoc_BV = assoc_BB + cut_AB;
259 
260  if (!cut_AB)
261  return 0;
262  else return cut_AB/assoc_AV + cut_AB/assoc_BV;
263 
264 }
265 
266 
267 /*---------------------------------------------------------------
268  exactBisection
269  ---------------------------------------------------------------*/
270 template <class GRAPH_MATRIX, typename num_t>
271 void CGraphPartitioner<GRAPH_MATRIX,num_t>::exactBisection(
272  GRAPH_MATRIX &in_A,
273  vector_uint &out_part1,
274  vector_uint &out_part2,
275  num_t &out_cut_value,
276  bool forceSimetry )
277 {
278  size_t nodeCount; // Nodes count
279  size_t i,j;
280  GRAPH_MATRIX Adj;
281  vector_bool partition, bestPartition;
282  vector_uint part1,part2;
283  num_t partCutValue, bestCutValue = std::numeric_limits<num_t>::max();
284  bool end = false;
285 
286  // Check matrix is square:
287  if (in_A.getColCount() != (nodeCount = in_A.getRowCount()) )
288  THROW_EXCEPTION("Weights matrix is not square!!");
289 
290  ASSERT_(nodeCount>=2);
291 
292  // forceSimetry?
293  if (forceSimetry)
294  {
295  Adj.setSize(nodeCount,nodeCount);
296  for (i=0;i<nodeCount;i++)
297  for (j=i;j<nodeCount;j++)
298  Adj(i,j)=Adj(j,i)= 0.5f*(in_A(i,j)+in_A(j,i));
299  }
300  else Adj = in_A;
301 
302 
303  // Brute force: compute all posible partitions:
304  //-----------------------------------------------------------------
305 
306  // First combination: 1000...0000
307  // Last combination: 1111...1110
308  partition.clear();
309  partition.resize(nodeCount,false);
310  partition[0] = true;
311 
312  while (!end)
313  {
314  // Build partitions from binary vector:
315  part1.clear();
316  part2.clear();
317 
318  for (i=0;i<nodeCount;i++)
319  {
320  if (partition[i])
321  part2.push_back(i);
322  else part1.push_back(i);
323  }
324 
325  // Compute the n-cut:
326  partCutValue = nCut(Adj,part1,part2);
327 
328  if (partCutValue<bestCutValue)
329  {
330  bestCutValue = partCutValue;
331  bestPartition = partition;
332  }
333 
334  // Next combo in the binary vector:
335  i = 0;
336  bool carry = false;
337  do
338  {
339  carry = partition[i];
340  partition[i]=!partition[i];
341  i++;
342  } while ( carry && i<nodeCount );
343 
344  // End criterion:
345  end = true;
346  for (i=0;end && i<nodeCount;i++)
347  end = end && partition[i];
348 
349  }; // End of while
350 
351 
352  // Return the best partition:
353  out_cut_value = bestCutValue;
354 
355  out_part1.clear();
356  out_part2.clear();
357 
358  for (i=0;i<nodeCount;i++)
359  {
360  if (bestPartition[i])
361  out_part2.push_back(i);
362  else out_part1.push_back(i);
363  }
364 
365 }
366 
367 } // end NS
368 } // end NS
mrpt::vector_uint
std::vector< uint32_t > vector_uint
Definition: types_simple.h:28
mrpt::graphs::CGraphPartitioner::RecursiveSpectralPartition
static void RecursiveSpectralPartition(GRAPH_MATRIX &in_A, std::vector< vector_uint > &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.
Definition: CGraphPartitioner_impl.h:101
eigenVectors
EIGEN_STRONG_INLINE bool eigenVectors(MATRIX1 &eVecs, MATRIX2 &eVals) const
[For square matrices only] Compute the eigenvectors and eigenvalues (sorted), both returned as matric...
THROW_EXCEPTION
#define THROW_EXCEPTION(msg)
Definition: maps/CColouredPointsMap.h:154
mrpt::vector_bool
std::vector< bool > vector_bool
A type for passing a vector of bools.
Definition: types_simple.h:29
MRPT_END
#define MRPT_END
Definition: maps/CColouredPointsMap.h:370
end
GLuint GLuint end
Definition: glext.h:3512
mrpt::graphs::CGraphPartitioner::SpectralBisection
static void SpectralBisection(GRAPH_MATRIX &in_A, vector_uint &out_part1, vector_uint &out_part2, num_t &out_cut_value, bool forceSimetry=true)
Performs the spectral bisection of a graph.
Definition: CGraphPartitioner_impl.h:24
mrpt::format
std::string BASE_IMPEXP format(const char *fmt,...) MRPT_printf_format_check(1
A std::string version of C sprintf.
Definition: format.cpp:21
MRPT_START
#define MRPT_START
Definition: maps/CColouredPointsMap.h:366
mrpt
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
Definition: CParticleFilter.h:16
mrpt::graphs::CGraphPartitioner::exactBisection
static void exactBisection(GRAPH_MATRIX &in_A, vector_uint &out_part1, vector_uint &out_part2, num_t &out_cut_value, bool forceSimetry=true)
Performs an EXACT minimum n-Cut graph bisection, (Use CGraphPartitioner::SpectralBisection for a fast...
Definition: CGraphPartitioner_impl.h:271
ASSERT_
#define ASSERT_(f)
Definition: maps/CColouredPointsMap.h:278
eigenValues
EIGEN_STRONG_INLINE void eigenValues(VECTOR &eVals) const
[For square matrices only] Compute the eigenvectors and eigenvalues (sorted), and return only the eig...
Definition: eigen_plugins.h:647
mrpt::graphs::CGraphPartitioner::nCut
static num_t nCut(const GRAPH_MATRIX &in_A, const vector_uint &in_part1, const vector_uint &in_part2)
Returns the normaliced cut of a graph, given its adjacency matrix A and a bisection: ...
Definition: CGraphPartitioner_impl.h:229
mean
EIGEN_STRONG_INLINE double mean() const
Computes the mean of the entire matrix.
Definition: eigen_plugins.h:355


Page generated by Doxygen 1.8.14 for MRPT 1.5.7 Git: 5902e14cc Wed Apr 24 15:04:01 2019 +0200 at lun oct 28 01:39:17 CET 2019