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

Go to the documentation of this file.
00001 
00002 //
00003 //  Module           : LocalCutSeparator.cpp
00004 //  Description      : Separate a local cut.
00005 //  Author           : Steffen Rebennack
00006 //  Email            : srebenac@ix.urz.uni-heidelberg.de;
00007 //                     steffen.rebennack@web.de
00008 //  Copyright        : (C) 2006 by the University of Heidelberg
00009 //  Created on       : Wed Apr 05 15:13:59 2006
00010 //  Last modified by : -
00011 //  Last modified on : -
00012 //  Update count     : 0
00013 //
00015 //
00016 //  Date        Name            Changes/Extensions
00017 //  ----        ----            ------------------
00018 //
00020 
00021 #include "LocalCutSeparator.hh"
00022 #include "abacus/array.h"
00023 #include "abacus/row.h"
00024 #include "abacus/cplexif.h"
00025 #include "GeneralConstraint.hh"
00026 #include "Graph.hh"
00027 #include "MaxCliquesSeparator.hh"
00028 #include "StableSetLPSolution.hh"
00029 #include "StableSetMaster.hh"
00030 
00031 #include <vector>
00032 
00033 using namespace std;
00034 
00035 
00036 //
00037 // The struct NodesAndWeights2 stores an index of a node and a weight.
00038 //
00039 struct NodesAndWeights2 {
00040     int    inode;         // index of node
00041     double wweight;       // weighted weight of node
00042 };// struct NodesAndWeights2
00043 
00044 
00045 //
00046 // This operator is used by the sort function in lowerBoundHeuristic to get
00047 // a sorted vector of variables with increasing weight.
00048 //
00049 bool operator<(const NodesAndWeights2 &a, const NodesAndWeights2 &b) {
00050         return a.wweight > b.wweight;
00051 }
00052 
00053 
00054 
00055 // -----------------------------------------------------------------------------
00056 // --------------- M e t h o d s  ( p u b l i c ) ------------------------------
00057 // -----------------------------------------------------------------------------
00058 
00059 
00060 // -----------------------------------------------------------------------------
00061 // C o n s t r u c t o r
00062 //
00063 // The main part of this constructor is the initialization of the abstract
00064 // template class ABA_SEPARATOR. The second argument is true because we do
00065 // not want double constraints in the buffer and the nuber 500 states the
00066 // maximla number of cuttings planes which are stored.
00067 // -----------------------------------------------------------------------------
00068 LocalCutSeparator::LocalCutSeparator(ABA_MASTER *master,
00069                         StableSetLPSolution *solution, Graph *theGraph):
00070     ABA_SEPARATOR<ABA_CONSTRAINT, ABA_VARIABLE> (solution, true, 500),
00071     master_(master),
00072     lpSolution(solution),
00073     itsGraph(theGraph)
00074 {
00075 }
00076 
00077 
00078 // -----------------------------------------------------------------------------
00079 // D e s t r u c t o r
00080 // -----------------------------------------------------------------------------
00081 LocalCutSeparator::~LocalCutSeparator()
00082 {
00083     for (int i=0; i<adjacentList.size(); i++) {
00084         adjacentList[i].clear();
00085     }
00086 }
00087 
00088 
00089 // -----------------------------------------------------------------------------
00090 // s e p a r a t e
00091 //
00092 // This separation routine computes a subgraph and with help of a LP it creates
00093 // a cuttinf plane. This is called local cut. Only one cut will be generated!
00094 // -----------------------------------------------------------------------------
00095 void LocalCutSeparator::separate() {
00096 
00097     //
00098     // Constrcut a small graph. The graph has to be connected and the number of
00099     // Stable Set solutions should not exceed 10.000. The reason is that the
00100     // number of Stable Set solutions is the number of constraints in a LP.
00101     // The function getPartitialGraph() returns true if a connected graph could
00102     // be found. Otherwise we will stop this separation.
00103     //
00104     if (getPartitialGraph()) {
00105 
00106         //
00107         // Enumerate all Stable Sets of the small graph.
00108         //
00109 
00110         int   i, j, k;          // Loop variables;
00111         int   size;
00112         bool  allFound = false; // If this is true, all Stable Sets have been
00113                                 // computed.
00114         bool  noStableSet;      // Indicates whether a 0-1 sequence is a Stable
00115                                 // Set or not.
00116         // Store the last computed 0-1 sequence, even if it is no Stable Set.
00117         short *oldSolution = new short[numberOfNodes];
00118         // Store a new Stable Set.
00119         short *newSolution = new short[numberOfNodes];
00120         // Store all computed Stable Sets.
00121         ABA_BUFFER<short*> allSolutions(master_, numberOfNodes * numberOfNodes);
00122 
00123         // Compute the first Stable Set wich is the 0 sequence and only the
00124         // last vertex is set to 1.
00125         // -> Only the last vertex is in the Stable Set.
00126         for (i=0; i<numberOfNodes-1; i++) {
00127             oldSolution[i] = 0;
00128             newSolution[i] = 0;
00129         }
00130         oldSolution[numberOfNodes -1] = 1;
00131         newSolution[numberOfNodes -1] = 1;
00132         // Store this Stable Set.
00133         allSolutions.push(newSolution);
00134 
00135 
00136 
00137         //
00138         // Compute all Stable Set.
00139         // The algorithm computes in each while-loop one 0-1 sequence. This
00140         // sequence is unique determined by the last 0-1 sequence.
00141         // To be sure that this sequence is a Stable Set we have to check it.
00142         // If the check is positive we store the new solution.
00143         // The while-loop is finished, ater we have computed all Stable Sets.
00144         //
00145         while (!allFound) {
00146 
00147             // Allocate memory for the new solution.
00148             newSolution = new short[numberOfNodes];
00149             // i is the last vertex.
00150             i = numberOfNodes - 1;
00151             // There are two cases:
00152             // 1.) The last vertex has value 0.
00153             // 2.) The last vertex has value 1.
00154             if (oldSolution[i] == 0) {
00155                 // This is the first case.
00156                 // The next sequence differs from the last one only in the
00157                 // last element: It is set to 1.
00158                 newSolution[i] = 1;
00159                 oldSolution[i] = 1;
00160                 i--;
00161                 // Copy all the other values.
00162                 while (i >= 0) {
00163                     newSolution[i] = oldSolution[i];
00164                     i--;
00165                 }// while
00166             }// if
00167             else {
00168 
00169                 // The second case.
00170                 // We search for the next vertex which has value 0.
00171                 // The search is done backwards.
00172                 while (oldSolution[i] == 1) {
00173                     // Reset all the vertecies to value 0.
00174                     newSolution[i] = 0;
00175                     oldSolution[i] = 0;
00176                     i--;
00177                     // If i is equal to -1 we know that the last 0-1 sequence
00178                     // was the hole graph. In this case we can stop because
00179                     // we have computed all Stable Sets if we add the empty set.
00180                     if (i == -1) {
00181                         allFound = true;
00182                         // Add the empty set to the set of Stable Set solutions.
00183                         if (allSolutions.full()) {
00184                             allSolutions.realloc(allSolutions.size() + 1);
00185                         }
00186                         allSolutions.push(newSolution);
00187                         break;
00188                     }
00189                 } // while
00190                 if (allFound) {
00191                     continue;
00192                 }
00193                 // The new solution selects vertex i into the Stable Set.
00194                 // All vertices which are greater than i are not in the set.
00195                 newSolution[i] = 1;
00196                 oldSolution[i] = 1;
00197                 i--;
00198                 while (i >= 0) {
00199                     newSolution[i] = oldSolution[i];
00200                     i--;
00201                 }
00202             }// else
00203 
00204             // At this point we have computed a new 0-1 sequence.
00205             // The next step is to check if this sequence is a Stable Set.
00206 
00207             noStableSet = false;
00208             // Check each pair of vertecies of the Stable Set if they are
00209             // adjazent.
00210             for (i=0; i<numberOfNodes; i++) {
00211                 for (j=i+1; j<numberOfNodes; j++) {
00212                     k = 0;
00213                     // Look if the two nodes are in the set.
00214                     if ((newSolution[i] == 1) && (newSolution[j] == 1)) {
00215                         // Check if node i and node j are adjazent.
00216                         // We exploit that the adjazens list is upward sorted.
00217                         size = adjacentList[i].size();
00218                         while((k < size) && (j > adjacentList[i][k])) {
00219                             k++;
00220                         }// while
00221                         // This index k is decisive.
00222                         if ((k < size) && (j == adjacentList[i][k])) {
00223                             // -> This is no Stable Set because node i and j
00224                             //    are adjazent.
00225                             noStableSet = true;
00226                             // Update the 0-1 sequence.
00227                             // It can easily be seen that there is no 0-1
00228                             // sequence which is a Stable Set and have the first
00229                             // j vertecies in common with this sequence.
00230                             // Therefore we can jumpt to the set which has all
00231                             // ones after the j vertex.
00232                             for (k=j+1; k<numberOfNodes; k++) {
00233                                 oldSolution[k] = 1;
00234                             }// for
00235                             break;
00236                         }// if
00237                     }// if
00238                 }// for
00239                 if (noStableSet) {
00240                     break;
00241                 }// if
00242             }// for
00243             if (noStableSet) {
00244                 delete [] newSolution;
00245             }// if
00246             else {
00247                 // The sequence is a Stable Set and therefore we store it.
00248                 // Check if we have enough memory. If not we have to allocate
00249                 // some additional.
00250                 if (allSolutions.full()) {
00251                     allSolutions.realloc(allSolutions.size()*2);
00252                 }// if
00253                 allSolutions.push(newSolution);
00254             }// else
00255         }// while
00256 
00257         // Clean up.
00258         delete [] oldSolution;
00259 
00260 
00261 
00262         //
00263         // Create and solve LP.
00264         //
00265 
00266 
00267         // The number of constraints of the LP is the number of Stable Set
00268         // solutions.
00269         int numberOfColumns = allSolutions.number();
00270         // The number of non zero elements in a constraint.
00271         int numberOfNonZeros = 0;
00272         // Right hand side of a constraint.
00273         double rhs;
00274         // Solution value of the last LP for each node.
00275         // Be careful with the indizes. They have not to be equal to the one
00276         // of the small graph!
00277         double *solutionValue = lpSolution->zVal();
00278         // Inner point of the Stable Set polytop.
00279         double innerPoint[numberOfNodes];
00280         // Coefficint to calculate the inner point.
00281         double coeff = 1.0 / allSolutions.number();
00282         // A pointer to an ABACUS constraint.
00283         ABA_ROW *constraint;
00284         // Buffer the non zero elements of a constraint.
00285         ABA_BUFFER<int> nonZeroElementsBuf(master_, numberOfColumns);
00286         // Store the coefficients of each buffered vertex.
00287         ABA_BUFFER<double> coefficientsBuf(master_, numberOfColumns);
00288         // Store the vertecies of a constraint...
00289         ABA_ARRAY<int> *nonZeroElements;
00290         // ...and their coefficients.
00291         ABA_ARRAY<double> *coefficients;
00292         // The coefficients of the target function. They are initialized
00293         // with value 1.
00294         ABA_ARRAY<double> objectiveFunctionCoefficient(master_,
00295                                                        numberOfColumns, 1);
00296         // Lower bound of each variable is 0.
00297         ABA_ARRAY<double> lowerBound(master_, numberOfColumns, 0);
00298         // There is no upper bound for the variables.
00299         ABA_ARRAY<double> upperBound(master_, numberOfColumns, 1000000);
00300         // Store all constraints.
00301         // There are number of nodes constraints:
00302         // For each component of the solution vector one.
00303         ABA_ARRAY<ABA_ROW*> allConstraints(master_, numberOfNodes);
00304         
00305 
00306         // Construct interior point as a convex combination of all Stable Sets.
00307         for (i=0; i<numberOfNodes; i++) {
00308             innerPoint[i] = 0;
00309         }
00310         for (i=0; i<allSolutions.number()-1; i++) {
00311             for (j=0; j<numberOfNodes; j++) {
00312                 innerPoint[j] += coeff * allSolutions[i][j];
00313             }
00314         }
00315 
00316         // Construct all constraints.
00317         for (i=0; i<numberOfNodes; i++) {
00318             for (j=0; j<numberOfColumns; j++) {
00319                 if (allSolutions[j][i] - innerPoint[i] != 0) {
00320                     nonZeroElementsBuf.push(j);
00321                     coefficientsBuf.push(allSolutions[j][i] - innerPoint[i]);
00322                     numberOfNonZeros++;
00323                 }// if
00324             }// for
00325 
00326             nonZeroElements = new ABA_ARRAY<int>(master_, nonZeroElementsBuf);
00327             coefficients    = new ABA_ARRAY<double>(master_, coefficientsBuf);
00328             nonZeroElementsBuf.clear();
00329             coefficientsBuf.clear();
00330 
00331             // Compute the right hand side of this constraint.
00332             rhs = solutionValue[nodes[i]] - innerPoint[i];
00333 
00334             // The new constraint.
00335             constraint = new ABA_ROW(master_,
00336                                      numberOfNonZeros,
00337                                      (*nonZeroElements),
00338                                      (*coefficients),
00339                                      ABA_CSENSE::Equal,
00340                                      rhs);
00341             // Store this constraint.
00342             allConstraints[i] = constraint;
00343         }// for
00344 
00345 
00346         //
00347         // Optimize
00348         //
00349         ABA_CPLEXIF lp(master_ );
00350         lp.initialize(ABA_OPTSENSE::Min,
00351                       numberOfNodes,
00352                       numberOfNodes,
00353                       numberOfColumns,
00354                       numberOfColumns,
00355                       objectiveFunctionCoefficient,
00356                       lowerBound,
00357                       upperBound,
00358                       allConstraints);
00359 
00360         // Optimize and check if optimal solution could be computed.
00361         if (!(ABA_LP::Optimal == lp.optimize(ABA_LP::Primal))) {
00362             master_->err() << "ERROR in LocalCutSeparator::separate()." << endl;
00363             master_->err() << "LP could not be dolved to optimality." << endl;
00364             exit(master_->Fatal);
00365         }
00366 
00367         double eps = master_->machineEps(); // Machine precision.
00368 
00369         //
00370         // Compute cut.
00371         //
00372         if (lp.value() > (1 + eps)) {
00373 
00374             double a;                   // Value of a dual variable
00375             double b;                   // Value of a dual variable.
00376             double gcd;                 // Greates common divisor of two double.
00377             int    dualMult;            // A multiplicator.
00378             int    dualVariable[numberOfNodes]; // Dual variables.
00379             int    rhsConstraint;       // Right hand side of the cut
00380                                         // constraint.
00381             vector<int> vertices;       // Vertices of the cut constraint.
00382             vector<int> dualCoeff;      // Coefficients of the variables of the
00383                                         // cut constraints.
00384 
00385 
00386             i = 0;
00387             a = lp.yVal(i);
00388             dualVariable[i] = -1;
00389             while (((-eps < a) && (a < eps))
00390                    || (a <= -100000000)
00391                    || (a >=  100000000)
00392                   ) {
00393                 dualVariable[i] = 0;
00394                 i++;
00395                 if (i == numberOfNodes) {
00396                     master_->err() << "ERROR in LocalCutSeparator::separate().";
00397                     master_->err() << "\nAll dual variables are zero." << endl;
00398                     exit(master_->Fatal);
00399                 }
00400                 a = lp.yVal(i);
00401                 dualVariable[i] = -1;
00402             }// while
00403 
00404             i++;
00405 
00406             b = lp.yVal(i);
00407             dualVariable[i] = -1;
00408             while (((-eps < b) && (b < eps))
00409                    || (b <= -100000000)
00410                    || (b >=  100000000)
00411                   ) {
00412                 dualVariable[i] = 0;
00413                 i++;
00414                 if (i == numberOfNodes) {
00415                     break;
00416                 }
00417                 b = lp.yVal(i);
00418                 dualVariable[i] = -1;
00419             }// while
00420 
00421             gcd = calculateGCD(a, b);
00422 
00423             i++;
00424             while (i < numberOfNodes) {
00425                 a = lp.yVal(i);
00426                 if (((-eps < a) && (a < eps))
00427                     || (a <= -100000000)
00428                     || (a >=  100000000)
00429                    ) {
00430                     dualVariable[i] = 0;
00431                 }
00432                 else {
00433                     gcd = calculateGCD(a, gcd);
00434                     dualVariable[i] = -1;
00435                 }
00436                 i++;
00437             }// while
00438 
00439             for (i=0; i<numberOfNodes; i++) {
00440                 if (dualVariable[i] == -1) {
00441                     a = lp.yVal(i) / gcd;
00442                     dualMult = (int) ceil(a);
00443                     if ((dualMult - a) > 0.001) {
00444                         dualMult -= 1;
00445                     }
00446                     dualVariable[i] = dualMult;
00447                 }
00448             }// for
00449 
00450             // Calculate rhs.
00451             rhs = 1.0 / gcd;
00452             for (i=0; i<numberOfNodes; i++) {
00453                 rhs += innerPoint[i] * dualVariable[i];
00454             }
00455 
00456             rhsConstraint = (int)floor(rhs + 0.001);
00457 
00458             // Vector vertices is sorted beacuse vector nodes is sorted.
00459             for (i=0; i<numberOfNodes; i++) {
00460                 if (dualVariable[i] != 0) {
00461                     vertices.push_back(nodes[i]);
00462                     dualCoeff.push_back(dualVariable[i]);
00463                 }
00464             }// for
00465 
00466             GeneralConstraint *generalConstr;
00467             generalConstr = new GeneralConstraint(master_,
00468                                                   &vertices,
00469                                                   &dualCoeff,
00470                                                   rhsConstraint,
00471                                                   itsGraph);
00472             cutFound(generalConstr);
00473 
00474         }
00475 
00476         // Clean up.
00477         for (i=0; i<allSolutions.number(); i++) {
00478             delete [] allSolutions[i];
00479         }
00480 
00481 
00482     }// if (getPartitialGraph())
00483     else {
00484         master_->out() << "Local cut separation: Subgraph could ";
00485         master_->out() << "not be computed." << endl;
00486     }
00487 
00488 }// void separate()
00489 
00490 
00491 
00492 // -----------------------------------------------------------------------------
00493 // --------------- M e t h o d s  ( p r i v a t e ) ----------------------------
00494 // -----------------------------------------------------------------------------
00495 
00496 
00497 // -----------------------------------------------------------------------------
00498 // g e t P a r t i t i a l G r a p h
00499 //
00500 // Construct a small graph where the local cut can be computed.
00501 // -----------------------------------------------------------------------------
00502 bool LocalCutSeparator::getPartitialGraph() {
00503 
00504     int i, j, k;
00505     double eps = master_->machineEps(); // Machine precision.
00506     double *solution = lpSolution->zVal();
00507     int    graphSize = itsGraph->numberOfNodes();
00508     vector<NodesAndWeights2> fractionalValue; // Save all fractional variables.
00509 
00510     // Get all fractional variables.
00511     for (i=0; i<graphSize; i++ ) {
00512         if ((solution[i] < -eps)
00513             || (solution[i] > 1 + eps)
00514             || ((solution[i] < 1 -eps) && (solution[i] > eps))
00515            ) {
00516             // This variable is fractional.
00517             NodesAndWeights2 a;
00518             a.inode   = i;
00519             a.wweight = solution[i];
00520             fractionalValue.push_back(a);
00521         }
00522     }
00523 
00524     // Sort the fractional variables:
00525     // the value with the highest weight is the last in this vector.
00526     sort(fractionalValue.begin(), fractionalValue.end());
00527     nodes.clear();
00528 
00529     int sizeOfGraph  = fractionalValue.size() < MAX_NUMBER_OF_NODES
00530                        ? fractionalValue.size()
00531                        : MAX_NUMBER_OF_NODES;
00532     bool found = false;
00533     int cnt;
00534     int maxNum;
00535     int indexOne;
00536     int indexTwo;
00537     int node = fractionalValue.back().inode;
00538     nodes.push_back(node);
00539     fractionalValue.pop_back();
00540 
00541     //
00542     // Increase until the partitial graph is big enough.
00543     //
00544     while (nodes.size() < sizeOfGraph) {
00545         if (found) {
00546             break;
00547         }
00548 
00549         cnt    = 0;
00550         found  = false;
00551         maxNum = stableSetMaster()->getRandomNumber(10) < 4 ? 1 : 2;
00552         indexOne = stableSetMaster()->getRandomNumber(fractionalValue.size());
00553 
00554         for (i=0; i<fractionalValue.size(); i++) {
00555             k = 0;
00556             indexTwo = (i + indexOne) % fractionalValue.size();
00557             if (itsGraph->isMemberOfAdjacentList(node, k, // 
00558                                         fractionalValue[indexTwo].inode)) {
00559                 nodes.push_back(fractionalValue[indexTwo].inode);
00560                 fractionalValue.erase(fractionalValue.begin() + indexTwo);
00561                 cnt++;
00562                 if (cnt == maxNum) {
00563                     break;
00564                 }
00565             }
00566         }
00567         if (cnt == 0) {
00568             if (nodes.size() == 1) {
00569                 if (fractionalValue.size() < MIN_NUMBER_OF_NODES) {
00570                     return false;
00571                 }
00572                 nodes.pop_back();
00573                 node = fractionalValue.back().inode;
00574                 nodes.push_back(node);
00575                 fractionalValue.pop_back();
00576                 // Update max size.
00577                 sizeOfGraph = fractionalValue.size() < MAX_NUMBER_OF_NODES
00578                               ? fractionalValue.size()
00579                               : MAX_NUMBER_OF_NODES;
00580             }
00581             else {
00582                 // Get index of current node.
00583                 k = nodes.size() - 1;
00584                 while (k > 0) {
00585                     if (node == nodes[k]) {
00586                         break;
00587                     }
00588                     k--;
00589                 }
00590 
00591                 if ((k == 0) && (node != nodes[0])) {
00592                     master_->err() << "LocalCutSeparator::";
00593                     master_->err() << "getPartitialGraph()" << endl;
00594                     master_->err() << "node could not be found" << endl;
00595                     exit(master_->Fatal);
00596                 }
00597                 if (k != 0) {
00598                     node = nodes[k-1];
00599                 }
00600                 else { // k == 0
00601 
00602                     if (nodes.size() >= MIN_NUMBER_OF_NODES) {
00603                         found = true;
00604                     }
00605                     else {
00606                         if (fractionalValue.size() < MIN_NUMBER_OF_NODES) {
00607                             sizeOfGraph = fractionalValue.size()
00608                                           < MAX_NUMBER_OF_NODES
00609                                           ? fractionalValue.size()
00610                                           : MAX_NUMBER_OF_NODES;
00611                             return false;
00612                         }
00613                         else {
00614                             nodes.clear();
00615                             node = fractionalValue.back().inode;
00616                             nodes.push_back(node);
00617                             fractionalValue.pop_back();
00618                         }
00619                     }
00620                 }
00621             }
00622         }
00623         else {
00624             node = nodes.back();
00625         }
00626     }
00627 
00628     // A partitional graph could be found.
00629     // Create the adjacent list and store the size of the graph.
00630     numberOfNodes = nodes.size();
00631     sort(nodes.begin(), nodes.end());
00632     adjacentList.resize(numberOfNodes);
00633     for (i=0; i<numberOfNodes; i++) {
00634         k = 0;
00635         for (j=i+1; j<numberOfNodes; j++) {
00636             if (itsGraph->isMemberOfAdjacentList(nodes[i], k, nodes[j])) {
00637                 adjacentList[i].push_back(j);
00638                 adjacentList[j].push_back(i);
00639             }
00640         }
00641     }
00642 
00643     return true;
00644 
00645 }// bool getPartitialGraph()
00646 
00647 
00648 // -----------------------------------------------------------------------------
00649 // c a l c u l a t e G C D
00650 // -----------------------------------------------------------------------------
00651 double LocalCutSeparator::calculateGCD(double a, double b) {
00652 
00653     double eps = master_->machineEps(); // Machine precision.
00654     while ((a > (b + eps)) || (a < (b - eps))) {
00655         if (a > b) {
00656             a -= b;
00657         }
00658         else {
00659             b -= a;
00660         }
00661     }
00662 
00663     return a;
00664 }
00665 
00666 
00667 // -----------------------------------------------------------------------------
00668 // s t a b l e S e t M a s t e r
00669 //
00670 // This function returns a pointer to the corresponding object of the class
00671 // StableSetMaster.
00672 // -----------------------------------------------------------------------------
00673 StableSetMaster *LocalCutSeparator::stableSetMaster() {
00674     return (StableSetMaster *) master_;
00675 }
00676 
00677 
00678 
00679 
00680 
00681 
00682 

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