modules/legacy/src/clique.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 "precomp.hpp"
  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 chosen 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 chosen 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