C:/Programme/home/Administrator/CD/stableSetCode/ModkSeparator.cpp

Go to the documentation of this file.
00001 
00002 //
00003 //  Module           : ModkSeparator.cpp
00004 //  Description      :
00005 //  Author           : Lars Fricke
00006 //  Email            : lars.fricke@iwr.uni-heidelberg.de
00007 //  Copyright        : University of Heidelberg
00008 //  Created on       : Tue Nov 8 2005
00009 //  Last modified by : Rebennack
00010 //  Last modified on : Wed Apr 05 15:51:14 2006
00011 //  Update count     : 2
00012 //
00014 //
00015 //  Date        Name            Changes/Extensions
00016 //  Apr 04      Rebennack       Disable output; Clean up memory.
00017 //  Apr 05      Rebennack       Fitting to programming guides.
00018 //
00020 
00021 #include "ModkSeparator.defs"
00022 #include "ModkSeparator.hh"
00023 #include <time.h>
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include <stdlib.h>
00027 
00028 
00029 // -----------------------------------------------------------------------------
00030 // C o n s t r u c t o r
00031 // 
00032 // Needs the constraints with slack zero.
00033 // -----------------------------------------------------------------------------
00034 ModkSeparator::ModkSeparator(int modfac, int nConstraints, int nVariables,
00035                              int* nNonZeros, int** NonZeroIndices,
00036                              double** NonZeroCoeff, double* rhs,
00037                              double* lpSolution):
00038     ordGen(false),
00039     constraintsGen(false),
00040     congruenceSystemGen(0),
00041     colsGen(false),
00042     generatedCutsGen(false)
00043 {
00044         K = modfac;
00045         // Order of elements
00046 
00047         if( K == 2 ){
00048           Ord = new int[2];
00049           ordGen = true;
00050           Ord[0]=0;
00051           Ord[1]=1;
00052         }
00053         else if( K == 3 ){
00054           Ord = new int[3];
00055           ordGen = true;
00056           Ord[0]=0;
00057           Ord[1]=1;
00058           Ord[2]=2;
00059         }
00060         else if( K == 5 ){
00061           Ord = new int[5];
00062           ordGen = true;
00063           Ord[0]=0;
00064           Ord[1]=1;
00065           Ord[2]=3;
00066           Ord[3]=2;
00067           Ord[4]=4;
00068         }
00069         else if( K == 7 ){
00070           Ord = new int[7];
00071           ordGen = true;
00072           Ord[0] = 0;
00073           Ord[1] = 1;
00074           Ord[2] = 4;
00075           Ord[3] = 5;
00076           Ord[4] = 2;
00077           Ord[5] = 3;
00078           Ord[6] = 6;
00079         }
00080         else{
00081           printf("%d is not a supported modulo factor\n", K);
00082           exit(1);
00083         }
00084 
00085 
00086         if(IMPORTANTVARIABLES){
00087           int curVar = 0;
00088 
00089           // Analyse variable distribution
00090           // Reset Count
00091           int* count = new int[nVariables];
00092           for(int i=0; i<nVariables; i++)
00093             count[i] = 0; 
00094           // Count variable occurences
00095           for(int i=0; i<nConstraints; i++)
00096             for(int j=0; j<nNonZeros[i]; j++)
00097                 count[NonZeroIndices[i][j]]++;
00098           // Greedy - take variable with most occurences            
00099           for(int i=0; i<nVariables; i++){
00100             if(count[i] > count[curVar]){
00101               curVar = i;
00102             } 
00103           }
00104           // Start with the variable with most occurences
00105 
00106           constraintsGen = true;
00107           NonZeroI = new int*[nConstraints]; 
00108           NonZeroC = new double*[nConstraints];   
00109           nNonZ = new int[nConstraints];
00110           RHS = new double[nConstraints];
00111           
00112           int maxCon = nVariables + MAXGAP;
00113           int* taken = new int[nConstraints];
00114           for(int i=0; i<nConstraints; i++)
00115             taken[i] = 0;
00116           bool doBreak = false;
00117           bool doContinue = false;
00118           nCon = 0;
00119           int nConOld = nCon;
00120         
00121           while(1){
00122             nConOld = nCon;
00123             //printf("Searching for occurences of variable %d\n", curVar);
00124             for(int i=0; i<nConstraints; i++){
00125               if(nCon == maxCon){
00126                 //printf("Maximum number of constrains reached: %d\n", nCon);
00127                 doBreak = true;
00128                 break;
00129               }
00130               if(taken[i] == 0){
00131                 for(int j=0; j<nNonZeros[i]; j++){
00132                   if(NonZeroIndices[i][j] == curVar){
00133                     //printf("Added constraint %d\n", i);
00134                     // Add constraint to pool
00135                     taken[i] = 1;
00136                     nNonZ[nCon] = nNonZeros[i];
00137                     NonZeroI[nCon] = NonZeroIndices[i];
00138                     NonZeroC[nCon] = NonZeroCoeff[i];
00139                     RHS[nCon] = rhs[i];
00140                     nCon++;
00141                   } // end if
00142                 } // end for
00143               } // end if
00144             } // end for               
00145             if(doBreak) // break initiated 
00146               break; // break while        
00147 
00148             if(nConOld == nCon){ // No new constraints found
00149               doContinue = false;
00150               for(int i=0; i<nVariables; i++){
00151                 if(count[i] != -1){
00152                   //printf("Changed variable to %d\n", i);
00153                   curVar = i;   // Take the next available variable
00154                   count[curVar] = -1;
00155                   doContinue = true;
00156                   break;     // and start new search
00157                 }
00158               }
00159               if(doContinue)
00160                 continue;
00161               else{
00162                 //printf("Nothing more to find\n");
00163                 break; // No variable found so lets stop the search
00164               }
00165             }
00166             // Reset Count
00167             count[curVar] = -1;
00168             for(int i=0; i<nVariables; i++){
00169               if(count[i] != -1)
00170                 count[i] = 0;
00171             }
00172             // Count variable occurences
00173             for(int i=0; i<nCon; i++){
00174               for(int j=0; j<nNonZ[i]; j++){
00175                 if(count[NonZeroI[i][j]] >= 0)
00176                   count[NonZeroI[i][j]]++;
00177               }
00178             }
00179             // Greedy - take variable with most occurences            
00180             
00181             for(int i=0; i<nVariables; i++){
00182               if(count[i] > count[curVar]){
00183                 curVar = i;
00184               } 
00185             } 
00186                  
00187               
00188           } // end while
00189         
00190           // Clean up.
00191           delete [] count;
00192           delete [] taken;      
00193 
00194         } // end if IMPORTANTVARIABLES
00195         else{
00196           // Save pointers to important stuff
00197           //printf("Using all constraints\n");
00198           NonZeroI = NonZeroIndices;
00199           NonZeroC = NonZeroCoeff;
00200           RHS = rhs;
00201           nNonZ = nNonZeros;
00202           nCon = nConstraints;
00203         }
00204         
00205         nVar = nVariables;
00206         lpSol = lpSolution;
00207 
00208         nRows = nVar + 1; // + 1 for original rhs
00209         nCols = nCon + 1; // + 1 for congruence rhs
00210         
00211         // Allocate memory for the congruence system
00212         
00213         CongruenceSystem = new int*[nRows];
00214         congruenceSystemGen = nRows;
00215 
00216         for(int i=0; i<=nVar; i++)
00217                 CongruenceSystem[i] = new int[nCols];
00218 
00219         // Create the mod-K congruence system
00220         
00221         int entry=0;
00222 
00223         for(int i=0; i<nCon; i++){
00224                 for(int j=0; j<nVar; j++){
00225                         entry = 0;
00226                         for(int k=0; k<(int) nNonZ[i]; k++){
00227                                 if(NonZeroI[i][k] == j){ // is non-zero
00228                 entry = (K + (((int) NonZeroC[i][k]) % K)) % K; // modulo-K
00229                                 }
00230                         }
00231                         CongruenceSystem[j][i] = entry;
00232                 }       
00233         }
00234 
00235         for(int i=0; i<nRows-1; i++)
00236                 CongruenceSystem[i][nCols-1] = 0;
00237         CongruenceSystem[nRows-1][nCols-1] = K-1;
00238 
00239         for(int i=0; i<nCon; i++){
00240                 CongruenceSystem[nVar][i] = (K + (((int) RHS[i]) % K)) % K; 
00241         }
00242 }// constructor
00243         
00244 
00245 // -----------------------------------------------------------------------------
00246 // D e s t r u c t o r
00247 // -----------------------------------------------------------------------------
00248 ModkSeparator::~ModkSeparator()
00249 {
00250     // Clean up.
00251     if (ordGen) {
00252         delete [] Ord;
00253     }
00254     if (constraintsGen) {
00255         // Do not delete the entries.
00256         // They are only pointers to a object and not copied!
00257         delete [] NonZeroI;
00258         delete [] NonZeroC;
00259         delete [] nNonZ;
00260         delete [] RHS;
00261     }
00262     if (congruenceSystemGen) {
00263         for (int i=0; i<congruenceSystemGen; i++) {
00264             delete [] CongruenceSystem[i];
00265 //          CongruenceSystem[i] = NULL; 
00266         }
00267         delete [] CongruenceSystem;
00268     }
00269     if (colsGen) {
00270         delete [] Cols;
00271     }
00272     if (generatedCutsGen) {
00273         for (int i=0; i<nGen; i++) {
00274             delete [] GeneratedCuts[i];
00275         }
00276         delete [] GeneratedCuts;
00277     }
00278 }
00279 
00280 
00281 // -----------------------------------------------------------------------------
00282 // s e p a r a t e
00283 // -----------------------------------------------------------------------------
00284 int ModkSeparator::separate(int nMaxCuts)
00285 {    
00286         time_t begin = time(NULL);
00287 
00288         bool doBreak = false;
00289         int* swap;
00290         int swap2=0;
00291         int idSize= ((nRows) < (nCols)) ? (nRows) : (nCols);
00292         int factor = 0;
00293 
00294         Cols = new int[nCols];
00295         colsGen = true;
00296         for(int i=0; i<nCols; i++)
00297                 Cols[i] = i;    
00298 
00302         if(RANDOMIZE){
00303           srand(time(NULL));
00304           for (int i=0; i<nCols; i++) {
00305             int r = rand() % nCols;  // generate a random position
00306             int temp = Cols[i]; Cols[i] = Cols[r]; Cols[r] = temp;
00307           }
00308         }
00309 
00310         swap = CongruenceSystem[0];
00311         CongruenceSystem[0] = CongruenceSystem[nRows-1];
00312         CongruenceSystem[nRows-1] = swap;
00313 
00314         
00315         // Gaussian Elimination on the congruence system (forwards)
00316         
00317         int scalFac = 0;
00318         
00319         for(int i=0; (i<nRows) && (i<nCols-1); i++){ // for every row i
00320                 // Check if there is anything greater than 1
00321                 if( CongruenceSystem[i][Cols[i]] > 1 ){
00322                         scalFac = Ord[CongruenceSystem[i][Cols[i]]];
00323                         // If so, multiply row with Ord[N] so that a
00324                         // 1 is on the diagonal
00325                         for(int j=i; j<nCols; j++){    // for every column geq i
00326 CongruenceSystem[i][Cols[j]] = (scalFac*CongruenceSystem[i][Cols[j]]) % K; 
00327                         }
00328                 }
00329                                 
00330                 // Pivot search
00331                 // there is a diagonal 0
00332                 if( CongruenceSystem[i][Cols[i]] == 0 ){         
00333                     // Row pivot search
00334                     for(int k=i+1; k<nRows; k++){ 
00335                         // so search in the other rows             
00336                         if( CongruenceSystem[k][Cols[i]] != 0 ){
00337                             // there is a non-zero  
00338                             swap = CongruenceSystem[i];     // swap the rows
00339                             CongruenceSystem[i] = CongruenceSystem[k];
00340                             CongruenceSystem[k] = swap;
00341                             break;      
00342                             // and break out of the for-loop                                            
00343                                 }
00344                         } // end of row search
00345                 }       
00346                 
00347                 // Still a diagonoal 0 ?
00348                 if( CongruenceSystem[i][Cols[i]] == 0 ){ 
00349                         // Column pivot search
00350                         
00351                         for(int k=i+1; k<nCols-1; k++){
00352                                 // there is a non-zero
00353                                 if( CongruenceSystem[i][Cols[k]] != 0 ){ 
00354                                         swap2 = Cols[k];
00355                                         // swap the index vector entries
00356                                         Cols[k] = Cols[i];
00357                                         Cols[i] = swap2;
00358                                         break;
00359                                         // and break out of for-loop                                            
00360                                 }       
00361                         } 
00362                 } //end of column search
00363 
00364                 // Full pivot search                    
00365                 if( CongruenceSystem[i][Cols[i]] == 0 ){
00366                         // Still a diagonal 0 ?? Jesus tapdancing christ
00367                         doBreak = false;
00368                         for(int k=i+1; k<nRows; k++){            
00369                         // We now have to search the lower right quadrant
00370                         // for any 1
00371                                 for(int l=i+1; l<nCols-1; l++){
00372                                         if( CongruenceSystem[k][Cols[l]] != 0){
00373                                     // Finally !
00374                                    swap = CongruenceSystem[i];
00375                                    // swap the rows
00376                                    CongruenceSystem[i] = CongruenceSystem[k];
00377                                    CongruenceSystem[k] = swap;
00378                                    swap2 = Cols[i];                                     
00379                                    // and then swap the index vector entries
00380                                    Cols[i] = Cols[l];
00381                                    Cols[l] = swap2;     
00382                                    doBreak = true;
00383                                    // break out of the for-loop    
00384                                    break;             
00385                                         }
00386                                         // break out of the for-loop
00387                                         if(doBreak) break;     
00388                                 }
00389                         }
00390                 }  // end of full pivot search
00391                 // end of pivot search
00392 
00393                 // Check if there is anything greater than 1
00394                 if( CongruenceSystem[i][Cols[i]] > 1 ){
00395                         // If so, multiply the row with Ord[N]
00396                         scalFac = Ord[CongruenceSystem[i][Cols[i]]];
00397                         for(int j=i; j<nCols; j++){    
00398                                 // for every column geq i
00399 CongruenceSystem[i][Cols[j]] = (scalFac*CongruenceSystem[i][Cols[j]]) % K;
00400                         }
00401                 }               
00402 
00403                 // NOW THERE IS 1 IN THE DIAGONAL OR NOTHING 
00404 
00405                 if( CongruenceSystem[i][Cols[i]] == 1 ){ 
00406                 // is there a diagonal 1 ?      
00407                 // Elimination
00408                         for(int k=i+1; k<nRows; k++){                                   
00409                             if( CongruenceSystem[k][Cols[i]] != 0 ){
00410         // there is a non-zero, so we substract row i from row k (x-times)
00411                                 factor = CongruenceSystem[k][Cols[i]];
00412                                 for(int j=i; j<nCols; j++){    
00413                                 // for every column geq i
00414 CongruenceSystem[k][Cols[j]] = (K + ((CongruenceSystem[k][Cols[j]] - //
00415 factor*CongruenceSystem[i][Cols[j]])%K))%K; // modulo-K
00416                                         }
00417                                 }
00418                         }
00419                 }
00420                 else{  // no pivot found, elimination is over
00421                         idSize = i;     
00422                         break; // break out of for-loop
00423                 }
00424         } // enf of gaussian elimination
00425 
00426 
00427         // Backward Elimination
00428 
00429         for(int i=idSize-1; i>=0; i--){
00430                 for(int k=i-1; k>=0; k--){
00431                         if( CongruenceSystem[k][Cols[i]] != 0 ){
00432                                 factor = CongruenceSystem[k][Cols[i]];
00433                                 for(int j=i; j<nCols; j++){   
00434                                  // for every column geq i
00435 CongruenceSystem[k][Cols[j]] = (K + ((CongruenceSystem[k][Cols[j]] - //
00436 factor*CongruenceSystem[i][Cols[j]])%K))%K; // modulo-K
00437                                 }
00438                         }
00439                 }
00440         }
00441 
00442         // Check if the Congruence System has a solution
00443         
00444         for(int i=idSize; i<nRows; i++)
00445                 if(CongruenceSystem[i][nCols-1] != 0){
00446                         // BAD BAD BAD
00447                         // The congruence system has no solution
00448                         return 0;
00449                 }
00450 
00451         // Allocate memory for internal cut pool
00452 
00453         GeneratedCuts = new int*[nMaxCuts];
00454         generatedCutsGen = true;
00455 
00456         nGen = 0;
00457 
00458         //      
00459         // GENERATE CUTS
00460         //
00461 
00462         int checksum = 0;
00463         bool feasible = false;
00464 
00465 
00466         // The basic solution
00467 
00468         GeneratedCuts[nGen] = new int[nCols-1];
00469         generatedCutsGen = true;
00470         for(int l=0; l<nCols-1; l++){
00471                 if(l<idSize)
00472                         GeneratedCuts[nGen][Cols[l]] =
00473                                  CongruenceSystem[l][Cols[nCols-1]];
00474                 else
00475                         GeneratedCuts[nGen][Cols[l]] = 0;
00476         }
00477         nGen++;
00478 
00479         // Additional solutions
00480 
00481         bool cutGood = true;
00482         
00483 
00484         int curSum;
00485 
00486         time_t START = time(NULL);
00487         time_t END;
00488         for(int i=1; i<pow(2,nCols-idSize); i++){
00489           END = time(NULL);
00490           if(nGen == nMaxCuts || difftime(END,START) > TIMEOUT)
00491             break;
00492           
00493           GeneratedCuts[nGen] = new int[nCols-1];
00494           for(int l=0; l<nCols-1; l++)
00495             GeneratedCuts[nGen][l] = 0;
00496                                    
00497           for(int j=0; j<idSize; j++){
00498             
00499             curSum = 0;
00500             for(int k=idSize; k<nCols-1; k++){
00501               if( (i & (1 << (k-idSize))) ){
00502                 curSum = curSum + CongruenceSystem[j][Cols[k]];
00503                 GeneratedCuts[nGen][Cols[k]] = 1;
00504               } 
00505            }
00506 
00507 GeneratedCuts[nGen][Cols[j]] = (K + ((CongruenceSystem[j][Cols[nCols-1]] - curSum)%K )) %K;
00508           }
00509           cutGood = true;
00510 
00516           if(MAXDENSITY){
00517             int nNZ = 0;
00518 
00519             double* Cut = new double[nVar];
00520             
00521             for(int j=0; j<nVar; j++)
00522               Cut[j] = 0; // Just in case theres junk in the memory
00523 
00524             double Multiplier = 1.0/((double) K);
00525 
00526             for(int l=0; l<nCon; l++){
00527               if(GeneratedCuts[nGen][l] != 0){ // Column i takes part
00528                 for(int j=0; j<nNonZ[l]; j++){
00529                   Cut[NonZeroI[l][j]] = Cut[NonZeroI[l][j]] + //
00530                         Multiplier*GeneratedCuts[nGen][l]*NonZeroC[l][j]; 
00531                 }
00532               }
00533             }
00534 
00535             for(int l=0; l<nVar; l++){
00536               if(Cut[l] != 0)
00537                 nNZ++;
00538             }
00539             
00540             if(nNZ > MAXDENSITY){
00541               cutGood = false;
00542             }
00543             
00544             // Clean up
00545             delete [] Cut;
00546           }
00547 
00548           if(cutGood)
00549             nGen++;
00550         }
00551 
00552         // END OF CUT GENERATION
00553 
00554         time_t end = time(NULL);
00555 
00556         // END OF SEPARATION
00557 
00558         nCutsLeft = nGen;
00559         return nGen;
00560 
00561 }// separate(..)
00562 
00563 
00564 // -----------------------------------------------------------------------------
00565 // g e t C u t
00566 // -----------------------------------------------------------------------------
00567 int* ModkSeparator::getCut() {
00568 
00569         int nNonZerosCut=0; 
00570         int k=-1;
00571         double* Cut = new double[nVar];
00572         double CutRHS = 0;
00573         int* sparseCut;
00574 
00575         double entry=0;
00576         
00577         if(nCutsLeft){
00578                 nCutsLeft--;
00579                 int cI = nGen-(nCutsLeft+1);
00580         
00581                 // Construct cut from the congruence system solutions
00582         
00583                 for(int j=0; j<nVar; j++)
00584                         Cut[j] = 0; // Just in case theres junk in the memory
00585 
00586                 double Multiplier = 1.0/((double) K);
00587                 
00588 
00589                 for(int i=0; i<nCon; i++){
00590 
00591                         if(GeneratedCuts[cI][i] != 0){ // Column i takes part
00592                                 for(int j=0; j<nNonZ[i]; j++){
00593 Cut[NonZeroI[i][j]] = Cut[NonZeroI[i][j]] +//
00594  Multiplier*GeneratedCuts[cI][i]*NonZeroC[i][j]; 
00595                                 }
00596                                 CutRHS = CutRHS + Multiplier*GeneratedCuts[cI][i]*RHS[i]; 
00597                         }
00598                 }
00599         
00600               
00601                 double lhs=0;
00602 
00603                 // Convert cut to sparse format
00604                 
00605                 for(int i=0; i<nVar; i++)
00606                         if(Cut[i] != 0)
00607                                 nNonZerosCut++;
00608                 sparseCut = new int[(nNonZerosCut*2)+2];
00609                 sparseCut[0] = nNonZerosCut;
00610                 for(int i=0; i<nVar; i++){
00611                         if(Cut[i] != 0){                                
00612                                 lhs = lhs + ((int) Cut[i])*lpSol[i];
00613                                 k++;
00614                                 sparseCut[k+1] = i;
00615                                 sparseCut[k+1+nNonZerosCut] = (int) Cut[i];
00616                         }
00617                 }
00618                 sparseCut[(2*nNonZerosCut)+1] = (int) floor(CutRHS);
00619         
00620                 // Clean up.
00621                 delete [] Cut;
00622 
00623                 return sparseCut;
00624         }
00625         else {
00626             // Clean up.
00627             delete [] Cut;      
00628 
00629                 return NULL;
00630         }
00631 }// getCut()
00632 
00633 
00634 
00635 // -----------------------------------------------------------------------------
00636 // --------------- M e t h o d s  ( p r i v a t e ) ----------------------------
00637 // -----------------------------------------------------------------------------
00638 
00639 
00640 // -----------------------------------------------------------------------------
00641 // h w e i g h t
00642 // -----------------------------------------------------------------------------
00643 unsigned int ModkSeparator::hweight(unsigned int w)
00644 {
00645     unsigned int res = (w & 0x55555555) + ((w >> 1) & 0x55555555);
00646     res = (res & 0x33333333) + ((res >> 2) & 0x33333333);
00647     res = (res & 0x0F0F0F0F) + ((res >> 4) & 0x0F0F0F0F);
00648     res = (res & 0x00FF00FF) + ((res >> 8) & 0x00FF00FF);
00649     return (res & 0x0000FFFF) + ((res >> 16) & 0x0000FFFF);
00650 }
00651 
00652 

Generated on Fri Apr 28 15:49:59 2006 for Branch and Cut algorithm for the Maximum Stable Set Problem by  doxygen 1.4.6-NO