Update the changelog
[opencv] / cvaux / src / cvclique.cpp
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 //  By downloading, copying, installing or using the software you agree to this license.
6 //  If you do not agree to this license, do not download, install,
7 //  copy or use the software.
8 //
9 //
10 //                        Intel License Agreement
11 //                For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
15 //
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
18 //
19 //   * Redistribution's of source code must retain the above copyright notice,
20 //     this list of conditions and the following disclaimer.
21 //
22 //   * Redistribution's in binary form must reproduce the above copyright notice,
23 //     this list of conditions and the following disclaimer in the documentation
24 //     and/or other materials provided with the distribution.
25 //
26 //   * The name of Intel Corporation may not be used to endorse or promote products
27 //     derived from this software without specific prior written permission.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 //M*/
41
42 #include "_cvaux.h"
43
44 #if 0
45
46 #include <float.h>
47 #include <limits.h>
48 #include <stdio.h>
49
50
51 #include "_cvutils.h"
52 #include "_cvwrap.h"
53
54 /*typedef struct CvCliqueFinder
55 {   
56     CvGraph* graph;
57     int**    adj_matr;
58     int N; //graph size
59
60     // stacks, counters etc/
61     int k; //stack size
62     int* current_comp;
63     int** All;
64     
65     int* ne;
66     int* ce;
67     int* fixp; //node with minimal disconnections
68     int* nod;
69     int* s; //for selected candidate
70     int status;
71     int best_score;
72
73 } CvCliqueFinder;
74 */
75
76 #define GO 1
77 #define BACK 2
78 #define PEREBOR 3
79 #define NEXT PEREBOR
80 #define END 4
81
82
83 #define  CV_GET_ADJ_VTX( vertex, edge ) \
84 (                                       \
85     assert(edge->vtx[0]==vertex||edge->vtx[1] == vertex ), \
86     (edge->vtx[0] == vertex)?edge->vtx[1]:edge->vtx[0]     \
87 )
88
89
90 #define NUMBER( v ) ((v)->flags >> 1 )
91
92 void _MarkNodes( CvGraph* graph )
93 {
94     //set number of vertices to their flags
95     for( int i = 0; i < graph->total; i++ )
96     {
97         CvGraphVtx* ver = cvGetGraphVtx( graph, i );
98         if( ver )
99         {
100             ver->flags = i<<1;
101         }
102     }  
103 }
104
105 void _FillAdjMatrix( CvGraph* graph, int** connected, int reverse )
106 {
107     //assume all vertices are marked
108     for( int i = 0; i < graph->total; i++ )
109     {
110         for( int j = 0; j < graph->total; j++ )
111         {
112             connected[i][j] = 0|reverse;
113         } 
114         //memset( connected[i], 0, sizeof(int)*graph->total );
115         CvGraphVtx* ver = cvGetGraphVtx( graph, i );
116         if( ver )
117         {
118             connected[i][i] = 1;
119             for( CvGraphEdge* e = ver->first; e ; e = CV_NEXT_GRAPH_EDGE( e, ver ) )
120             {
121                CvGraphVtx* v = CV_GET_ADJ_VTX( ver, e );
122                connected[i][NUMBER(v)] = 1^reverse;
123             }
124         }
125     }
126 }
127
128
129 void cvStartFindCliques( CvGraph* graph, CvCliqueFinder* finder, int reverse, int weighted, int weighted_edges )
130 {
131     int i;
132
133     if (weighted)
134     {
135         finder->weighted = 1;
136         finder->best_weight = 0;
137         finder->vertex_weights = (float*)malloc( sizeof(float)*(graph->total+1));
138         finder->cur_weight = (float*)malloc( sizeof(float)*(graph->total+1));
139         finder->cand_weight = (float*)malloc( sizeof(float)*(graph->total+1));
140
141         finder->cur_weight[0] = 0;
142         finder->cand_weight[0] = 0;
143         for( i = 0 ; i < graph->total; i++ )
144         {
145             CvGraphWeightedVtx* ver = (CvGraphWeightedVtx*)cvGetGraphVtx( graph, i );
146             assert(ver);
147             assert(ver->weight>=0);
148             finder->vertex_weights[i] = ver->weight;
149             finder->cand_weight[0] += ver->weight;             
150         }
151     }
152     else finder->weighted = 0;
153
154     if (weighted_edges)
155     {
156         finder->weighted_edges = 1;
157         //finder->best_weight = 0;
158         finder->edge_weights = (float*)malloc( sizeof(float)*(graph->total)*(graph->total));
159         //finder->cur_weight = (float*)malloc( sizeof(float)*(graph->total+1));
160         //finder->cand_weight = (float*)malloc( sizeof(float)*(graph->total+1));
161
162         //finder->cur_weight[0] = 0;
163         //finder->cand_weight[0] = 0;
164         memset( finder->edge_weights, 0, sizeof(float)*(graph->total)*(graph->total) );
165         for( i = 0 ; i < graph->total; i++ )
166         {
167             CvGraphVtx* ver1 = cvGetGraphVtx( graph, i );
168             if(!ver1) continue;
169             for( int j = i ; j < graph->total; j++ )
170             {
171                 CvGraphVtx* ver2 = cvGetGraphVtx( graph, j );
172                 if(!ver2) continue;
173                 CvGraphEdge* edge = cvFindGraphEdgeByPtr( graph, ver1, ver2 );
174                 if( edge )
175                 {
176                     assert( ((CvGraphWeightedEdge*)edge)->weight >= 0 );
177                     finder->edge_weights[ i * graph->total + j ] = 
178                     finder->edge_weights[ j * graph->total + i ] = ((CvGraphWeightedEdge*)edge)->weight;
179                 }
180             }
181         }        
182     }
183     else finder->weighted_edges = 0;
184                                
185
186     //int* Compsub; //current component (vertex stack)
187     finder->k = 0; //counter of steps
188     int N = finder->N = graph->total;
189     finder->current_comp = new int[N];
190     finder->All = new int*[N];
191     for( i = 0 ; i < finder->N; i++ )
192     {
193         finder->All[i] = new int[N];
194     }
195
196     finder->ne = new int[N+1];
197     finder->ce = new int[N+1];
198     finder->fixp = new int[N+1]; //node with minimal disconnections
199     finder->nod = new int[N+1];
200     finder->s = new int[N+1]; //for selected candidate
201     
202     //form adj matrix
203     finder->adj_matr = new int*[N]; //assume filled with 0
204     for( i = 0 ; i < N; i++ )
205     {
206         finder->adj_matr[i] = new int[N];
207     }
208
209     //set number to vertices
210     _MarkNodes( graph );
211     _FillAdjMatrix( graph, finder->adj_matr, reverse );
212
213     //init all arrays
214     int k = finder->k = 0; //stack size
215     memset( finder->All[k], 0, sizeof(int) * N );
216     for( i = 0; i < N; i++ )  finder->All[k][i] = i;
217     finder->ne[0] = 0;
218     finder->ce[0] = N;
219
220     finder->status = GO;
221     finder->best_score = 0;
222
223 }   
224
225 void cvEndFindCliques( CvCliqueFinder* finder )
226 {
227     int i;
228     
229     //int* Compsub; //current component (vertex stack)
230     delete finder->current_comp;
231     for( i = 0 ; i < finder->N; i++ )
232     {
233         delete finder->All[i];
234     }
235     delete finder->All;
236     
237     delete finder->ne;
238     delete finder->ce;
239     delete finder->fixp; //node with minimal disconnections
240     delete finder->nod;
241     delete finder->s; //for selected candidate
242     
243     //delete adj matrix
244     for( i = 0 ; i < finder->N; i++ )
245     {
246         delete finder->adj_matr[i];
247     }  
248     delete finder->adj_matr;
249     
250     if(finder->weighted)
251     {
252         free(finder->vertex_weights);
253         free(finder->cur_weight);
254         free(finder->cand_weight); 
255     }
256     if(finder->weighted_edges)
257     {
258         free(finder->edge_weights);
259     }
260 }
261
262 int cvFindNextMaximalClique( CvCliqueFinder* finder ) 
263 {
264     int**  connected = finder->adj_matr;
265 //    int N = finder->N; //graph size
266
267     // stacks, counters etc/
268     int k = finder->k; //stack size
269     int* Compsub = finder->current_comp;
270     int** All = finder->All;
271     
272     int* ne = finder->ne;
273     int* ce = finder->ce;
274     int* fixp = finder->fixp; //node with minimal disconnections
275     int* nod = finder->nod;
276     int* s = finder->s; //for selected candidate   
277     
278     //START
279     while( k >= 0)
280     {
281         int* old = All[k];
282         switch(finder->status)
283         {    
284         case GO://Forward step
285         /* we have all sets and will choose fixed point */ 
286             {   
287                 //check potential size of clique
288                 if( (!finder->weighted) && (k + ce[k] - ne[k] < finder->best_score) ) 
289                 {
290                     finder->status  = BACK;
291                     break;
292                 } 
293                 //check potential weight
294                 if( finder->weighted && !finder->weighted_edges &&
295                     finder->cur_weight[k] + finder->cand_weight[k] < finder->best_weight )
296                 {
297                     finder->status  = BACK;
298                     break;
299                 }
300
301                 int minnod = ce[k];
302                 nod[k] = 0;
303
304                 //for every vertex of All determine counter value and choose minimum
305                 for( int i = 0; i < ce[k] && minnod != 0; i++) 
306                 {   
307                     int p = old[i]; //current point
308                     int count = 0;  //counter
309                     int pos = 0;
310
311                     /* Count disconnections with candidates */
312                     for (int j = ne[k]; j < ce[k] && (count < minnod); j++) 
313                     {
314                         if ( !connected[p][old[j]] ) 
315                         {
316                             count++;
317                             /* Save position of potential candidate */
318                             pos = j;
319                         }
320                     }
321                     
322                     /* Test new minimum */
323                     if (count < minnod) 
324                     {
325                         fixp[k] = p;     //set current point as fixed
326                         minnod = count;  //new value for minnod
327                         if (i < ne[k])      //if current fixed point belongs to 'not'
328                         {
329                             s[k] = pos;     //s - selected candidate
330                         } 
331                         else 
332                         {
333                             s[k] = i;        //selected candidate is fixed point itself
334                             /* preincr */
335                             nod[k] = 1;      //nod is aux variable, 1 means fixp == s
336                         }
337                     }
338                 }//for
339
340                 nod[k] = minnod + nod[k];
341                 finder->status = NEXT;//go to backtrackcycle 
342             }
343             break;
344         case NEXT:
345             //here we will look for candidate to translate into not    
346             //s[k] now contains index of choosen candidate
347             {
348                 int* new_ = All[k+1];
349                 if( nod[k] != 0 )
350                 {
351                     //swap selected and first candidate
352                     int i, p = old[s[k]];
353                     old[s[k]] = old[ne[k]];
354                     int sel = old[ne[k]] = p;
355
356                     int newne = 0;
357                     //fill new set 'not'
358                     for ( i = 0; i < ne[k]; i++) 
359                     {
360                         if (connected[sel][old[i]]) 
361                         {
362                             new_[newne] = old[i];
363                             newne++;
364                             
365                         }
366                     }
367                     //fill new set 'candidate'
368                     int newce = newne;
369                     i++;//skip selected candidate
370                     
371                     float candweight = 0;
372                     
373                     for (; i < ce[k]; i++) 
374                     {
375                         if (connected[sel][old[i]]) 
376                         {
377                             new_[newce] = old[i];
378                             
379                             if( finder->weighted ) 
380                                 candweight += finder->vertex_weights[old[i]];                            
381                             
382                             newce++; 
383                         }
384                     }
385
386                     nod[k]--;
387
388                     //add selected to stack 
389                     Compsub[k] = sel;
390
391                     k++;
392                     assert( k <= finder->N );
393                     if( finder->weighted )
394                     {
395                         //update weights of current clique and candidates
396                         finder->cur_weight[k] = finder->cur_weight[k-1] + finder->vertex_weights[sel]; 
397                         finder->cand_weight[k] = candweight;
398                     }       
399                     if( finder->weighted_edges )
400                     {
401                         //update total weight by edge weights
402                         float added = 0;
403                         for( int ind = 0; ind < k-1; ind++ )
404                         {
405                             added += finder->edge_weights[ Compsub[ind] * finder->N + sel ];
406                         }
407                         finder->cur_weight[k] += added;
408                     }                    
409                                         
410                     //check if 'not' and 'cand' are both empty
411                     if( newce == 0 )
412                     {
413                         finder->best_score = MAX(finder->best_score, k );
414
415                         if( finder->weighted )
416                             finder->best_weight = MAX( finder->best_weight, finder->cur_weight[k] );
417                         /*FILE* file  = fopen("cliques.txt", "a" );
418
419                         for (int t=0; t<k; t++)
420                         {
421                           fprintf(file, "%d ", Compsub[t]);
422                         }
423                         fprintf(file, "\n");
424
425                         fclose(file);
426                         */
427                         
428                         //output new clique//************************
429                         finder->status = BACK;
430                         finder->k = k;
431
432                         return CLIQUE_FOUND;
433
434                     }
435                     else //check nonempty set of candidates
436                     if( newne < newce )
437                     {
438                         //go further
439                         ne[k] = newne;
440                         ce[k] = newce;
441                         finder->status  = GO;
442                         break;
443                     }
444                     
445                 }
446                 else 
447                     finder->status  = BACK;
448
449             }
450             break;
451
452         case BACK:
453             {         
454                 //decrease stack
455                 k--;         
456                 old = All[k];
457                 if( k < 0 ) break;
458
459                 //add to not
460                 ne[k]++;
461
462                 if( nod[k] > 0 )
463                 {
464                     //select next candidate
465                     for( s[k] = ne[k]; s[k] < ce[k]; s[k]++ )
466                     {
467                         if( !connected[fixp[k]][old[s[k]]])
468                             break;
469                     }
470                     assert( s[k] < ce[k] );
471                     finder->status = NEXT;
472                 }
473                 else
474                     finder->status = BACK;
475
476             }
477             break;
478         case END: assert(0);       
479
480         }
481     }//end while
482
483     finder->status = END;
484     return CLIQUE_END;
485 }
486
487
488                                      
489
490 void cvBronKerbosch( CvGraph* graph )
491 {
492     int* Compsub; //current component (vertex stack)
493     int k = 0; //counter of steps
494     int N = graph->total;
495     int i;
496     Compsub = new int[N];
497     int** All = new int*[N];
498     for( i = 0 ; i < N; i++ )
499     {
500         All[i] = new int[N];
501     }
502
503     int* ne = new int[N];
504     int* ce = new int[N];
505     int* fixp = new int[N]; //node with minimal disconnections
506     int* nod = new int[N];
507     int* s = new int[N]; //for selected candidate
508     
509     //form adj matrix
510     int** connected = new int*[N]; //assume filled with 0
511     for( i = 0 ; i < N; i++ )
512     {
513         connected[i] = new int[N];
514     }
515
516
517
518     //set number to vertices
519     _MarkNodes( graph );
520     _FillAdjMatrix( graph, connected, 0 );
521
522     //init all arrays
523     k = 0; //stack size
524     memset( All[k], 0, sizeof(int) * N );
525     for( i = 0; i < N; i++ )  All[k][i] = i;
526     ne[0] = 0;
527     ce[0] = N;
528
529     int status = GO;
530     int best_score = 0;
531
532     //START
533     while( k >= 0)
534     {
535         int* old = All[k];
536         switch(status)
537         {    
538         case GO://Forward step
539         /* we have all sets and will choose fixed point */ 
540             {   
541
542                 if( k + ce[k] - ne[k] < best_score ) 
543                 {
544                     status  = BACK;
545                     break;
546                 } 
547
548                 int minnod = ce[k];
549                 nod[k] = 0;
550
551                 //for every vertex of All determine counter value and choose minimum
552                 for( int i = 0; i < ce[k] && minnod != 0; i++) 
553                 {   
554                     int p = old[i]; //current point
555                     int count = 0;  //counter
556                     int pos = 0;
557
558                     /* Count disconnections with candidates */
559                     for (int j = ne[k]; j < ce[k] && (count < minnod); j++) 
560                     {
561                         if ( !connected[p][old[j]] ) 
562                         {
563                             count++;
564                             /* Save position of potential candidate */
565                             pos = j;
566                         }
567                     }
568                     
569                     /* Test new minimum */
570                     if (count < minnod) 
571                     {
572                         fixp[k] = p;     //set current point as fixed
573                         minnod = count;  //new value for minnod
574                         if (i < ne[k])      //if current fixed point belongs to 'not'
575                         {
576                             s[k] = pos;     //s - selected candidate
577                         } 
578                         else 
579                         {
580                             s[k] = i;        //selected candidate is fixed point itself
581                             /* preincr */
582                             nod[k] = 1;      //nod is aux variable, 1 means fixp == s
583                         }
584                     }
585                 }//for
586
587                 nod[k] = minnod + nod[k];
588                 status = NEXT;//go to backtrackcycle 
589             }
590             break;
591         case NEXT:
592             //here we will look for candidate to translate into not    
593             //s[k] now contains index of choosen candidate
594             {
595                 int* new_ = All[k+1];
596                 if( nod[k] != 0 )
597                 {
598                     //swap selected and first candidate
599                     int p = old[s[k]];
600                     old[s[k]] = old[ne[k]];
601                     int sel = old[ne[k]] = p;
602
603                     int newne = 0;
604                     //fill new set 'not'
605                     for ( i = 0; i < ne[k]; i++) 
606                     {
607                         if (connected[sel][old[i]]) 
608                         {
609                             new_[newne] = old[i];
610                             newne++;
611                             
612                         }
613                     }
614                     //fill new set 'candidate'
615                     int newce = newne;
616                     i++;//skip selected candidate
617                     for (; i < ce[k]; i++) 
618                     {
619                         if (connected[sel][old[i]]) 
620                         {
621                             new_[newce] = old[i];
622                             newce++;
623                         }
624                     }
625
626                     nod[k]--;
627
628                     //add selected to stack 
629                     Compsub[k] = sel;
630                     k++;
631
632                     //check if 'not' and 'cand' are both empty
633                     if( newce == 0 )
634                     {
635                         best_score = MAX(best_score, k );
636
637                         FILE* file  = fopen("cliques.txt", "a" );
638
639                         for (int t=0; t<k; t++)
640                         {
641                           fprintf(file, "%d ", Compsub[t]);
642                         }
643                         fprintf(file, "\n");
644
645                         fclose(file);
646
647                         /*for( int t = 0; t < k; t++ )
648                         {
649                             printf("%d ", Compsub[t] );
650                         }
651                         printf("\n"); */
652
653                         //printf("found %d\n", k);
654
655                         //output new clique//************************
656                         status = BACK;
657                     }
658                     else //check nonempty set of candidates
659                     if( newne < newce )
660                     {
661                         //go further
662                         ne[k] = newne;
663                         ce[k] = newce;
664                         status  = GO;
665                         break;
666                     }
667                     
668                 }
669                 else 
670                     status  = BACK;
671
672             }
673             break;
674
675         case BACK:
676             {         
677                 //decrease stack
678                 k--;         
679                 old = All[k];
680                 if( k < 0 ) break;
681
682                 //add to not
683                 ne[k]++;
684
685                 if( nod[k] > 0 )
686                 {
687                     //select next candidate
688                     for( s[k] = ne[k]; s[k] < ce[k]; s[k]++ )
689                     {
690                         if( !connected[fixp[k]][old[s[k]]])
691                             break;
692                     }
693                     assert( s[k] < ce[k] );
694                     status = NEXT;
695                 }
696                 else
697                     status = BACK;
698
699             }
700             break;
701        
702
703         }
704     }//end while
705
706 }//end cvBronKerbosch
707
708 #endif
709