00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00015
00016
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
00038
00039 struct NodesAndWeights2 {
00040 int inode;
00041 double wweight;
00042 };
00043
00044
00045
00046
00047
00048
00049 bool operator<(const NodesAndWeights2 &a, const NodesAndWeights2 &b) {
00050 return a.wweight > b.wweight;
00051 }
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
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
00080
00081 LocalCutSeparator::~LocalCutSeparator()
00082 {
00083 for (int i=0; i<adjacentList.size(); i++) {
00084 adjacentList[i].clear();
00085 }
00086 }
00087
00088
00089
00090
00091
00092
00093
00094
00095 void LocalCutSeparator::separate() {
00096
00097
00098
00099
00100
00101
00102
00103
00104 if (getPartitialGraph()) {
00105
00106
00107
00108
00109
00110 int i, j, k;
00111 int size;
00112 bool allFound = false;
00113
00114 bool noStableSet;
00115
00116
00117 short *oldSolution = new short[numberOfNodes];
00118
00119 short *newSolution = new short[numberOfNodes];
00120
00121 ABA_BUFFER<short*> allSolutions(master_, numberOfNodes * numberOfNodes);
00122
00123
00124
00125
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
00133 allSolutions.push(newSolution);
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145 while (!allFound) {
00146
00147
00148 newSolution = new short[numberOfNodes];
00149
00150 i = numberOfNodes - 1;
00151
00152
00153
00154 if (oldSolution[i] == 0) {
00155
00156
00157
00158 newSolution[i] = 1;
00159 oldSolution[i] = 1;
00160 i--;
00161
00162 while (i >= 0) {
00163 newSolution[i] = oldSolution[i];
00164 i--;
00165 }
00166 }
00167 else {
00168
00169
00170
00171
00172 while (oldSolution[i] == 1) {
00173
00174 newSolution[i] = 0;
00175 oldSolution[i] = 0;
00176 i--;
00177
00178
00179
00180 if (i == -1) {
00181 allFound = true;
00182
00183 if (allSolutions.full()) {
00184 allSolutions.realloc(allSolutions.size() + 1);
00185 }
00186 allSolutions.push(newSolution);
00187 break;
00188 }
00189 }
00190 if (allFound) {
00191 continue;
00192 }
00193
00194
00195 newSolution[i] = 1;
00196 oldSolution[i] = 1;
00197 i--;
00198 while (i >= 0) {
00199 newSolution[i] = oldSolution[i];
00200 i--;
00201 }
00202 }
00203
00204
00205
00206
00207 noStableSet = false;
00208
00209
00210 for (i=0; i<numberOfNodes; i++) {
00211 for (j=i+1; j<numberOfNodes; j++) {
00212 k = 0;
00213
00214 if ((newSolution[i] == 1) && (newSolution[j] == 1)) {
00215
00216
00217 size = adjacentList[i].size();
00218 while((k < size) && (j > adjacentList[i][k])) {
00219 k++;
00220 }
00221
00222 if ((k < size) && (j == adjacentList[i][k])) {
00223
00224
00225 noStableSet = true;
00226
00227
00228
00229
00230
00231
00232 for (k=j+1; k<numberOfNodes; k++) {
00233 oldSolution[k] = 1;
00234 }
00235 break;
00236 }
00237 }
00238 }
00239 if (noStableSet) {
00240 break;
00241 }
00242 }
00243 if (noStableSet) {
00244 delete [] newSolution;
00245 }
00246 else {
00247
00248
00249
00250 if (allSolutions.full()) {
00251 allSolutions.realloc(allSolutions.size()*2);
00252 }
00253 allSolutions.push(newSolution);
00254 }
00255 }
00256
00257
00258 delete [] oldSolution;
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 int numberOfColumns = allSolutions.number();
00270
00271 int numberOfNonZeros = 0;
00272
00273 double rhs;
00274
00275
00276
00277 double *solutionValue = lpSolution->zVal();
00278
00279 double innerPoint[numberOfNodes];
00280
00281 double coeff = 1.0 / allSolutions.number();
00282
00283 ABA_ROW *constraint;
00284
00285 ABA_BUFFER<int> nonZeroElementsBuf(master_, numberOfColumns);
00286
00287 ABA_BUFFER<double> coefficientsBuf(master_, numberOfColumns);
00288
00289 ABA_ARRAY<int> *nonZeroElements;
00290
00291 ABA_ARRAY<double> *coefficients;
00292
00293
00294 ABA_ARRAY<double> objectiveFunctionCoefficient(master_,
00295 numberOfColumns, 1);
00296
00297 ABA_ARRAY<double> lowerBound(master_, numberOfColumns, 0);
00298
00299 ABA_ARRAY<double> upperBound(master_, numberOfColumns, 1000000);
00300
00301
00302
00303 ABA_ARRAY<ABA_ROW*> allConstraints(master_, numberOfNodes);
00304
00305
00306
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
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 }
00324 }
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
00332 rhs = solutionValue[nodes[i]] - innerPoint[i];
00333
00334
00335 constraint = new ABA_ROW(master_,
00336 numberOfNonZeros,
00337 (*nonZeroElements),
00338 (*coefficients),
00339 ABA_CSENSE::Equal,
00340 rhs);
00341
00342 allConstraints[i] = constraint;
00343 }
00344
00345
00346
00347
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
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();
00368
00369
00370
00371
00372 if (lp.value() > (1 + eps)) {
00373
00374 double a;
00375 double b;
00376 double gcd;
00377 int dualMult;
00378 int dualVariable[numberOfNodes];
00379 int rhsConstraint;
00380
00381 vector<int> vertices;
00382 vector<int> dualCoeff;
00383
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 }
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 }
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 }
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 }
00449
00450
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
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 }
00465
00466 GeneralConstraint *generalConstr;
00467 generalConstr = new GeneralConstraint(master_,
00468 &vertices,
00469 &dualCoeff,
00470 rhsConstraint,
00471 itsGraph);
00472 cutFound(generalConstr);
00473
00474 }
00475
00476
00477 for (i=0; i<allSolutions.number(); i++) {
00478 delete [] allSolutions[i];
00479 }
00480
00481
00482 }
00483 else {
00484 master_->out() << "Local cut separation: Subgraph could ";
00485 master_->out() << "not be computed." << endl;
00486 }
00487
00488 }
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 bool LocalCutSeparator::getPartitialGraph() {
00503
00504 int i, j, k;
00505 double eps = master_->machineEps();
00506 double *solution = lpSolution->zVal();
00507 int graphSize = itsGraph->numberOfNodes();
00508 vector<NodesAndWeights2> fractionalValue;
00509
00510
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
00517 NodesAndWeights2 a;
00518 a.inode = i;
00519 a.wweight = solution[i];
00520 fractionalValue.push_back(a);
00521 }
00522 }
00523
00524
00525
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
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
00577 sizeOfGraph = fractionalValue.size() < MAX_NUMBER_OF_NODES
00578 ? fractionalValue.size()
00579 : MAX_NUMBER_OF_NODES;
00580 }
00581 else {
00582
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 {
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
00629
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 }
00646
00647
00648
00649
00650
00651 double LocalCutSeparator::calculateGCD(double a, double b) {
00652
00653 double eps = master_->machineEps();
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
00669
00670
00671
00672
00673 StableSetMaster *LocalCutSeparator::stableSetMaster() {
00674 return (StableSetMaster *) master_;
00675 }
00676
00677
00678
00679
00680
00681
00682