LocalCutSeparator Class Reference

#include <LocalCutSeparator.hh>

List of all members.

Public Member Functions

 LocalCutSeparator (ABA_MASTER *master, StableSetLPSolution *solution, Graph *theGraph)
 ~LocalCutSeparator ()
virtual void separate ()


Detailed Description

Separate a local cut.

Definition at line 38 of file LocalCutSeparator.hh.


Constructor & Destructor Documentation

LocalCutSeparator::LocalCutSeparator ABA_MASTER *  master,
StableSetLPSolution solution,
Graph theGraph
 

Constructor.

Parameters:
*master Pointer to the master of the problem.
*solution Pointer to the current LP-Solution to be separated.
*theGraph Pointer to the graph of the problem instance.

Definition at line 68 of file LocalCutSeparator.cpp.

00069                                                                        :
00070     ABA_SEPARATOR<ABA_CONSTRAINT, ABA_VARIABLE> (solution, true, 500),
00071     master_(master),
00072     lpSolution(solution),
00073     itsGraph(theGraph)
00074 {
00075 }

LocalCutSeparator::~LocalCutSeparator  ) 
 

Destructor.

Definition at line 81 of file LocalCutSeparator.cpp.

00082 {
00083     for (int i=0; i<adjacentList.size(); i++) {
00084         adjacentList[i].clear();
00085     }
00086 }


Member Function Documentation

void LocalCutSeparator::separate  )  [virtual]
 

This method provides the separation procedure itself and overrides the pure virtual function of ABA_SEPARATOR.

Parameters:
void 
Returns:
void

Definition at line 95 of file LocalCutSeparator.cpp.

00095                                  {
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()


The documentation for this class was generated from the following files:
Generated on Fri Apr 28 15:50:00 2006 for Branch and Cut algorithm for the Maximum Stable Set Problem by  doxygen 1.4.6-NO