Main MRPT website > C++ reference for MRPT 1.4.0
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-2016, 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
15namespace mrpt
16{
17namespace graphs
18{
19
20/*---------------------------------------------------------------
21 SpectralPartition
22 ---------------------------------------------------------------*/
23template <class GRAPH_MATRIX, typename num_t>
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 ---------------------------------------------------------------*/
100template <class GRAPH_MATRIX, typename num_t>
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{
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
223}
224
225/*---------------------------------------------------------------
226 nCut
227 ---------------------------------------------------------------*/
228template <class GRAPH_MATRIX, typename num_t>
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 ---------------------------------------------------------------*/
270template <class GRAPH_MATRIX, typename num_t>
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
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.
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:
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.
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...
EIGEN_STRONG_INLINE void eigenValues(VECTOR &eVals) const
[For square matrices only] Compute the eigenvectors and eigenvalues (sorted), and return only the eig...
EIGEN_STRONG_INLINE void eigenVectors(MATRIX1 &eVecs, MATRIX2 &eVals) const
[For square matrices only] Compute the eigenvectors and eigenvalues (sorted), both returned as matric...
EIGEN_STRONG_INLINE double mean() const
Computes the mean of the entire matrix.
EIGEN_STRONG_INLINE iterator end()
Definition: eigen_plugins.h:27
std::vector< bool > vector_bool
A type for passing a vector of bools.
Definition: types_simple.h:29
std::vector< uint32_t > vector_uint
Definition: types_simple.h:28
#define MRPT_START
Definition: mrpt_macros.h:349
#define ASSERT_(f)
Definition: mrpt_macros.h:261
#define MRPT_END
Definition: mrpt_macros.h:353
#define THROW_EXCEPTION(msg)
Definition: mrpt_macros.h:110
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
std::string BASE_IMPEXP format(const char *fmt,...) MRPT_printf_format_check(1
A std::string version of C sprintf.



Page generated by Doxygen 1.9.4 for MRPT 1.4.0 SVN: at Sun Aug 14 11:34:44 UTC 2022