#include <LocalCutSeparator.hh>
Public Member Functions | |
LocalCutSeparator (ABA_MASTER *master, StableSetLPSolution *solution, Graph *theGraph) | |
~LocalCutSeparator () | |
virtual void | separate () |
Definition at line 38 of file LocalCutSeparator.hh.
|
Constructor.
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 }
|
|
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 }
|
|
This method provides the separation procedure itself and overrides the pure virtual function of ABA_SEPARATOR.
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()
|