00001
00002
00003 #if defined(_MSC_VER)
00004
00005 # pragma warning(disable:4786)
00006 #endif
00007 #include <string>
00008
00009
00010 #include <cassert>
00011 #include <cmath>
00012 #include <cfloat>
00013
00014 #include "OsiSolverInterface.hpp"
00015 #include "CoinWarmStartBasis.hpp"
00016 #include "CoinPackedMatrix.hpp"
00017 #include "SbbCompareActual.hpp"
00018 #include "SbbBranchActual.hpp"
00019 #include "SbbHeuristic.hpp"
00020 #include "SbbModel.hpp"
00021 #include "SbbMessage.hpp"
00022 #include "OsiRowCut.hpp"
00023 #include "OsiColCut.hpp"
00024 #include "OsiRowCutDebugger.hpp"
00025 #include "OsiCuts.hpp"
00026 #include "SbbCountRowCut.hpp"
00027 #include "SbbCutGenerator.hpp"
00028
00029 #include "CglProbing.hpp"
00030
00031 #define COIN_USE_CLP
00032 #ifdef COIN_USE_CLP
00033
00034 #include "ClpPresolve.hpp"
00035
00036 #include "OsiClpSolverInterface.hpp"
00037 #endif
00038
00039 #include "CoinTime.hpp"
00040
00041 using namespace std;
00042
00043
00044
00045 #include <vector>
00046
00053 class tree : public std::vector <SbbNode *> {
00054
00055 SbbCompare comparison_;
00056
00057 public:
00058
00059 tree() {};
00060
00063
00065 void setComparison(SbbCompareBase &compare) {
00066 comparison_.test_ = &compare;
00067 make_heap(begin(), end(), comparison_);
00068 };
00069
00071 SbbNode * top() {return front();};
00072
00074 void push(SbbNode * x) {
00075 push_back(x);
00076 push_heap(begin(), end(), comparison_);
00077 };
00078
00080 void pop() {
00081 pop_heap(begin(), end(), comparison_);
00082 pop_back();
00083 };
00084
00086
00089
00096 void cleanTree(SbbModel * model, double cutoff)
00097 {
00098 int j;
00099 int nNodes = size();
00100 SbbNode ** nodeArray = new SbbNode * [nNodes];
00101 int * depth = new int [nNodes];
00102 int k=0;
00103 int kDelete=nNodes;
00104
00105
00106
00107
00108
00109 for (j=0;j<nNodes;j++) {
00110 SbbNode * node = top();
00111 pop();
00112 if (node->objectiveValue() >= cutoff) {
00113 nodeArray[--kDelete] = node;
00114 depth[kDelete] = node->depth();
00115 } else {
00116 nodeArray[k++]=node;
00117 }
00118 }
00119
00120
00121
00122 for (j=0;j<k;j++) { push(nodeArray[j]); }
00123
00124
00125
00126 CoinSort_2(depth+kDelete,depth+nNodes,nodeArray+kDelete);
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 for (j=nNodes-1;j >= kDelete;j--) {
00151 SbbNode * node = nodeArray[j];
00152 CoinWarmStartBasis *lastws = model->getEmptyBasis() ;
00153
00154 model->addCuts1(node,lastws);
00155
00156 int numberLeft = node->nodeInfo()->numberBranchesLeft();
00157 int i;
00158 for (i=0;i<model->currentNumberCuts();i++) {
00159
00160 CoinWarmStartBasis::Status status =
00161 lastws->getArtifStatus(i+model->numberRowsAtContinuous());
00162 if (status != CoinWarmStartBasis::basic&&
00163 model->addedCuts()[i]) {
00164 if (!model->addedCuts()[i]->decrement(numberLeft))
00165 delete model->addedCuts()[i];
00166 }
00167 }
00168
00169 node->nodeInfo()->throwAway();
00170 delete node ;
00171 delete lastws ;
00172 }
00173 delete [] nodeArray;
00174 delete [] depth;
00175 };
00176
00178
00179 };
00180
00181
00182
00183
00184
00185 namespace {
00186
00187
00188
00189
00190
00191 static int gcd(int a, int b)
00192 {
00193 int remainder = -1;
00194
00195 if(a > b) {
00196
00197 int temp = a;
00198 a = b;
00199 b = temp;
00200 }
00201
00202 if (!a) {
00203 if (b) {
00204 return b;
00205 } else {
00206 printf("**** gcd given two zeros!!\n");
00207 abort();
00208 }
00209 }
00210 while (remainder) {
00211 remainder = b % a;
00212 b = a;
00213 a = remainder;
00214 }
00215 return b;
00216 }
00217
00218
00219
00220 #ifdef CHECK_NODE_FULL
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234 void verifyTreeNodes (const tree &branchingTree, const SbbModel &model)
00235
00236 { printf("*** CHECKING tree after %d nodes\n",model.getNodeCount()) ;
00237
00238 int j ;
00239 int nNodes = branchingTree.size() ;
00240 # define MAXINFO 1000
00241 int *count = new int [MAXINFO] ;
00242 SbbNodeInfo **info = new SbbNodeInfo*[MAXINFO] ;
00243 int nInfo = 0 ;
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 for (j = 0 ; j < nNodes ; j++)
00255 { SbbNode *node = branchingTree[j] ;
00256 SbbNodeInfo *nodeInfo = node->nodeInfo() ;
00257 int change = node->nodeInfo()->numberBranchesLeft() ;
00258 assert(change) ;
00259 while (nodeInfo)
00260 { int k ;
00261 for (k = 0 ; k < nInfo ; k++)
00262 { if (nodeInfo == info[k]) break ; }
00263 if (k == nInfo)
00264 { assert(nInfo < MAXINFO) ;
00265 nInfo++ ;
00266 info[k] = nodeInfo ;
00267 count[k] = 0 ; }
00268 nodeInfo = nodeInfo->parent() ; } }
00269
00270
00271
00272
00273 for (j = 0 ; j < nInfo ; j++)
00274 { SbbNodeInfo *nodeInfo = info[j] ;
00275 nodeInfo = nodeInfo->parent() ;
00276 if (nodeInfo)
00277 { int k ;
00278 for (k = 0 ; k < nInfo ; k++)
00279 { if (nodeInfo == info[k]) break ; }
00280 assert (k < nInfo) ;
00281 count[k]++ ; } }
00282
00283
00284
00285
00286
00287
00288 for (j = 0;j < nInfo;j++) {
00289 SbbNodeInfo * nodeInfo = info[j] ;
00290 if (nodeInfo) {
00291 int k ;
00292 for (k = 0;k < nInfo;k++)
00293 if (nodeInfo == info[k])
00294 break ;
00295 printf("Nodeinfo %x - %d left, %d count\n",
00296 nodeInfo,
00297 nodeInfo->numberBranchesLeft(),
00298 nodeInfo->numberPointingToThis()) ;
00299 assert(nodeInfo->numberPointingToThis() ==
00300 count[k]+nodeInfo->numberBranchesLeft()) ; } }
00301
00302 delete [] count ;
00303 delete [] info ;
00304
00305 return ; }
00306
00307 #endif
00308
00309
00310
00311 #ifdef CHECK_CUT_COUNTS
00312
00313
00314
00315
00316 void verifyCutCounts (const tree &branchingTree, SbbModel &model)
00317
00318 { printf("*** CHECKING cuts after %d nodes\n",model.getNodeCount()) ;
00319
00320 int j ;
00321 int nNodes = branchingTree.size() ;
00322
00323
00324
00325
00326
00327
00328 for (j = 0 ; j < nNodes ; j++)
00329 { SbbNode *node = branchingTree[j] ;
00330 SbbNodeInfo * nodeInfo = node->nodeInfo() ;
00331 int change = node->nodeInfo()->numberBranchesLeft() ;
00332 assert(change) ;
00333 while (nodeInfo)
00334 { int k ;
00335 for (k = 0 ; k < nodeInfo->numberCuts() ; k++)
00336 { SbbCountRowCut *cut = nodeInfo->cuts()[k] ;
00337 if (cut) cut->tempNumber_ = 0; }
00338 nodeInfo = nodeInfo->parent() ; } }
00339
00340
00341
00342
00343
00344 for (j = 0 ; j < nNodes ; j++)
00345 { CoinWarmStartBasis *debugws = model.getEmptyBasis() ;
00346 SbbNode *node = branchingTree[j] ;
00347 SbbNodeInfo *nodeInfo = node->nodeInfo();
00348 int change = node->nodeInfo()->numberBranchesLeft() ;
00349 printf("Node %d %x (info %x) var %d way %d obj %g",j,node,
00350 node->nodeInfo(),node->variable(),node->way(),
00351 node->objectiveValue()) ;
00352
00353 model.addCuts1(node,debugws) ;
00354
00355 int i ;
00356 int numberRowsAtContinuous = model.numberRowsAtContinuous() ;
00357 SbbCountRowCut **addedCuts = model.addedCuts() ;
00358 for (i = 0 ; i < model.currentNumberCuts() ; i++)
00359 { CoinWarmStartBasis::Status status =
00360 debugws->getArtifStatus(i+numberRowsAtContinuous) ;
00361 if (status != CoinWarmStartBasis::basic && addedCuts[i])
00362 { addedCuts[i]->tempNumber_ += change ; } }
00363
00364 while (nodeInfo)
00365 { nodeInfo = nodeInfo->parent() ;
00366 if (nodeInfo) printf(" -> %x",nodeInfo); }
00367 printf("\n") ;
00368 delete debugws ; }
00369
00370
00371
00372
00373
00374
00375 for (j = 0 ; j < nNodes ; j++)
00376 { SbbNode *node = branchingTree[j] ;
00377 SbbNodeInfo *nodeInfo = node->nodeInfo();
00378 while (nodeInfo)
00379 { int k ;
00380 for (k = 0 ; k < nodeInfo->numberCuts() ; k++)
00381 { SbbCountRowCut *cut = nodeInfo->cuts()[k] ;
00382 if (cut && cut->tempNumber_ >= 0)
00383 { if (cut->tempNumber_ != cut->numberPointingToThis())
00384 printf("mismatch %x %d %x %d %d\n",nodeInfo,k,
00385 cut,cut->tempNumber_,cut->numberPointingToThis()) ;
00386 else
00387 printf(" match %x %d %x %d %d\n", nodeInfo,k,
00388 cut,cut->tempNumber_,cut->numberPointingToThis()) ;
00389 cut->tempNumber_ = -1 ; } }
00390 nodeInfo = nodeInfo->parent() ; } }
00391
00392 return ; }
00393
00394 #endif
00395
00396
00397
00398 void analyseObjective (SbbModel &model)
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 { const double *objective = model.getObjCoefficients() ;
00413 const double *lower = model.getColLower() ;
00414 const double *upper = model.getColUpper() ;
00415
00416
00417
00418
00419
00420
00421
00422
00423 double maximumCost = 0.0 ;
00424 bool possibleMultiple = true ;
00425 int iColumn ;
00426 int numberColumns = model.getNumCols() ;
00427 for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
00428 { if (upper[iColumn] > lower[iColumn]+1.0e-8)
00429 { if (model.isInteger(iColumn))
00430 maximumCost = max(maximumCost,fabs(objective[iColumn])) ;
00431 else if (objective[iColumn])
00432 possibleMultiple = false ; } }
00433 model.setIntParam(SbbModel::SbbFathomDiscipline,possibleMultiple) ;
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443 if (possibleMultiple)
00444 { int increment = 0 ;
00445 double multiplier = 2520.0 ;
00446 while (10.0*multiplier*maximumCost < 1.0e8)
00447 multiplier *= 10.0 ;
00448
00449 for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
00450 { if (upper[iColumn] > lower[iColumn]+1.0e-8)
00451 { if (model.isInteger(iColumn)&&objective[iColumn])
00452 { double value = fabs(objective[iColumn])*multiplier ;
00453 int nearest = (int) floor(value+0.5) ;
00454 if (fabs(value-floor(value+0.5)) > 1.0e-8)
00455 { increment = 0 ;
00456 break ; }
00457 else if (!increment)
00458 { increment = nearest ; }
00459 else
00460 { increment = gcd(increment,nearest) ; } } } }
00461
00462
00463
00464 if (increment)
00465 { double value = increment ;
00466 double cutoff = model.getDblParam(SbbModel::SbbCutoffIncrement) ;
00467 value /= multiplier ;
00468 if (value*0.999 > cutoff)
00469 { model.messageHandler()->message(SBB_INTEGERINCREMENT,
00470 model.messages())
00471 << value << CoinMessageEol ;
00472 model.setDblParam(SbbModel::SbbCutoffIncrement,value) ; } } }
00473
00474 return ; }
00475
00476
00477 }
00478
00479
00480
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507 void SbbModel::branchAndBound()
00508
00509 {
00510 # ifdef SBB_DEBUG
00511 std::string problemName ;
00512 solver_->getStrParam(OsiProbName,problemName) ;
00513 printf("Problem name - %s\n",problemName.c_str()) ;
00514 solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
00515 # endif
00516
00517
00518
00519 status_ = 0 ;
00520
00521
00522
00523
00524 findIntegers(false) ;
00525
00526
00527
00528
00529 synchronizeModel() ;
00530
00531
00532
00533 double time1 = CoinCpuTime() ;
00534
00535
00536
00537
00538
00539
00540
00541 bool feasible = resolve() ;
00542
00543
00544
00545
00546 if (!feasible)
00547 { handler_->message(SBB_INFEAS,messages_)<< CoinMessageEol ;
00548 status_ = 1 ;
00549 return ; }
00550
00551 # ifdef SBB_DEBUG
00552
00553
00554
00555
00556
00557 solver_->activateRowCutDebugger(problemName.c_str()) ;
00558 # endif
00559
00560
00561
00562
00563 bestObjective_ = 1.0e50 ;
00564 numberSolutions_ = 0 ;
00565 numberHeuristicSolutions_ = 0 ;
00566
00567 double cutoff=getCutoff() ;
00568 double direction = solver_->getObjSense() ;
00569 if (cutoff < 1.0e20&&direction<0.0)
00570 messageHandler()->message(SBB_CUTOFF_WARNING1,
00571 messages())
00572 << cutoff << -cutoff << CoinMessageEol ;
00573 if (cutoff > bestObjective_)
00574 cutoff = bestObjective_ ;
00575 setCutoff(cutoff) ;
00576
00577
00578
00579 int numberColumns = getNumCols() ;
00580 if (!currentSolution_)
00581 currentSolution_ = new double[numberColumns] ;
00582
00583
00584
00585
00586 continuousSolver_ = solver_->clone() ;
00587 numberRowsAtContinuous_ = getNumRows() ;
00588
00589
00590
00591
00592
00593 analyseObjective(*this) ;
00594
00595
00596
00597
00598
00599 int maximumWhich = 1000 ;
00600 int * whichGenerator = new int[maximumWhich] ;
00601 int currentNumberCuts = 0 ;
00602 maximumNumberCuts_ = 0 ;
00603 currentNumberCuts_ = 0 ;
00604 delete [] addedCuts_ ;
00605 addedCuts_ = NULL ;
00606
00607
00608
00609
00610 tree branchingTree ;
00611 SbbCompareObjective compareObjective ;
00612 SbbCompareDepth compareDepth ;
00613 SbbCompareDefault compareDefault(dblParam_[SbbInfeasibilityWeight]) ;
00614 if (!nodeCompare_)
00615 branchingTree.setComparison(compareDepth) ;
00616 else
00617 branchingTree.setComparison(*nodeCompare_) ;
00618
00619
00620
00621
00622 maximumDepth_ = 10 ;
00623 delete [] walkback_ ;
00624 walkback_ = new SbbNodeInfo * [maximumDepth_] ;
00625
00626
00627
00628 double * lowerBefore = new double [numberColumns] ;
00629 double * upperBefore = new double [numberColumns] ;
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645 OsiCuts cuts ;
00646 int anyAction = -1 ;
00647 int numberOldActiveCuts = 0 ;
00648 int numberNewCuts = 0 ;
00649 { int iObject ;
00650 int preferredWay ;
00651 double otherWay ;
00652 int numberUnsatisfied = 0 ;
00653 double integerTolerance = getIntegerTolerance() ;
00654
00655 memcpy(currentSolution_,solver_->getColSolution(),
00656 numberColumns*sizeof(double)) ;
00657
00658 for (iObject = 0 ; iObject < numberObjects_ ; iObject++)
00659 { double infeasibility =
00660 object_[iObject]->infeasibility(preferredWay,otherWay) ;
00661 if (infeasibility > integerTolerance) numberUnsatisfied++ ; }
00662 if (numberUnsatisfied)
00663 { feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
00664 NULL,numberOldActiveCuts,numberNewCuts,
00665 maximumWhich, whichGenerator) ; } }
00666 currentNumberCuts_ = numberNewCuts ;
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682 bool resolved = false ;
00683 SbbNode *newNode = NULL ;
00684
00685 if (feasible)
00686 { newNode = new SbbNode ;
00687 newNode->setObjectiveValue(direction*solver_->getObjValue()) ;
00688 anyAction = -1 ;
00689 while (anyAction == -1)
00690 { anyAction = newNode->chooseBranch(this,NULL) ;
00691 if (anyAction == -1)
00692 { feasible = resolve() ;
00693 resolved = true ;
00694 # ifdef SBB_DEBUG
00695 printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
00696 solver_->getObjValue(),
00697 solver_->getNumRows()) ;
00698 # endif
00699 if (!feasible) anyAction = -2 ; }
00700 if (anyAction == -2)
00701 { delete newNode ;
00702 newNode = NULL ;
00703 feasible = false ; } } }
00704
00705
00706
00707
00708
00709 assert (!newNode || newNode->objectiveValue() < getCutoff()) ;
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720 CoinWarmStartBasis *lastws = 0 ;
00721 if (feasible && newNode->variable() >= 0)
00722 { if (resolved)
00723 { bool needValidSolution = (newNode->variable() < 0) ;
00724 takeOffCuts(cuts,whichGenerator,numberOldActiveCuts,numberNewCuts,
00725 needValidSolution) ;
00726 # ifdef CHECK_CUT_COUNTS
00727 { printf("Number of rows after chooseBranch fix (root)"
00728 "(active only) %d\n",
00729 numberRowsAtContinuous_+numberNewCuts+numberOldActiveCuts) ;
00730 const CoinWarmStartBasis* debugws =
00731 dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
00732 debugws->print() ;
00733 delete debugws ; }
00734 # endif
00735 }
00736 newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
00737 newNode->nodeInfo()->addCuts(cuts,
00738 newNode->numberBranches(),whichGenerator) ;
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748 if (basis_) delete basis_ ;
00749 basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
00750 if (lastws) delete lastws ;
00751 lastws = dynamic_cast<CoinWarmStartBasis*>(basis_->clone()) ; }
00752
00753
00754
00755 continuousObjective_ = 0.0 ;
00756 continuousInfeasibilities_ = 0 ;
00757 if (newNode)
00758 { continuousObjective_ = newNode->objectiveValue() ;
00759 continuousInfeasibilities_ = newNode->numberUnsatisfied() ; }
00760
00761
00762
00763 { int i ;
00764 for (i = 0;i < numberObjects_;i++)
00765 object_[i]->resetBounds() ; }
00766
00767
00768
00769 numberIterations_ = 0 ;
00770 numberNodes_ = 0 ;
00771 bool stoppedOnGap = false ;
00772
00773
00774
00775
00776
00777
00778
00779
00780 if (newNode) {
00781 if (newNode->variable() >= 0) {
00782 newNode->initializeInfo() ;
00783 branchingTree.push(newNode) ;
00784 # ifdef CHECK_NODE
00785 printf("Node %x on tree\n",newNode) ;
00786 # endif
00787 } else {
00788
00789 double objectiveValue = newNode->objectiveValue();
00790 setBestSolution(SBB_SOLUTION,objectiveValue,
00791 solver_->getColSolution()) ;
00792 delete newNode ;
00793 newNode = NULL ;
00794 }
00795 }
00796
00797 if (printFrequency_ <= 0) {
00798 printFrequency_ = 1000 ;
00799 if (getNumCols() > 2000)
00800 printFrequency_ = 100 ;
00801 }
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812 double bestValue = 0.0 ;
00813 while (!branchingTree.empty())
00814 { if (cutoff > getCutoff()) {
00815
00816 branchingTree.cleanTree(this, getCutoff()) ;
00817 if (!nodeCompare_) {
00818 if (!compareDefault.getWeight()) {
00819
00820 double costPerInteger =
00821 (bestObjective_-continuousObjective_)/
00822 ((double) continuousInfeasibilities_) ;
00823 compareDefault.setWeight(0.98*costPerInteger) ;
00824
00825
00826 }
00827 branchingTree.setComparison(compareDefault) ;
00828
00829 } else {
00830 nodeCompare_->newSolution(this) ;
00831 nodeCompare_->newSolution(this,continuousObjective_,
00832 continuousInfeasibilities_) ;
00833 branchingTree.setComparison(*nodeCompare_) ;
00834 }
00835 if (branchingTree.empty())
00836 break;
00837 }
00838 cutoff = getCutoff() ;
00839
00840
00841
00842
00843
00844
00845 if ((numberNodes_%1000) == 0&&nodeCompare_)
00846 nodeCompare_->every1000Nodes(this, numberNodes_) ;
00847 if ((numberNodes_%printFrequency_) == 0) {
00848 int j ;
00849 int nNodes = branchingTree.size() ;
00850 bestValue = 1.0e100 ;
00851 for (j = 0;j < nNodes;j++) {
00852 SbbNode * node = branchingTree[j] ;
00853 if (node->objectiveValue() < bestValue)
00854 bestValue = node->objectiveValue() ;
00855 }
00856 messageHandler()->message(SBB_STATUS,messages())
00857 << numberNodes_<< nNodes<< bestObjective_<< bestValue
00858 << CoinMessageEol ;
00859
00860 if (bestObjective_-bestValue < dblParam_[SbbAllowableGap]) {
00861 stoppedOnGap = true ;
00862 }
00863 }
00864
00865 # ifdef CHECK_NODE_FULL
00866 verifyTreeNodes(branchingTree,*this) ;
00867 # endif
00868 # ifdef CHECK_CUT_COUNTS
00869 verifyCutCounts(branchingTree,*this) ;
00870 # endif
00871
00872
00873
00874
00875
00876
00877
00878 SbbNode *node = branchingTree.top() ;
00879 branchingTree.pop() ;
00880
00881 # ifdef CHECK_NODE
00882
00883
00884
00885
00886
00887 printf("Node %x popped from tree - %d left, %d count\n",node,
00888 node->nodeInfo()->numberBranchesLeft(),
00889 node->nodeInfo()->numberPointingToThis()) ;
00890 printf("\tdepth = %d, z = %g, unsat = %d, var = %d.\n",
00891 node->depth(),node->objectiveValue(),
00892 node->numberUnsatisfied(),
00893 integerVariable_[node->variable()]) ;
00894 # endif
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907 SbbNodeInfo * nodeInfo = node->nodeInfo() ;
00908 newNode = NULL ;
00909 if (!addCuts(node,lastws))
00910 { int i ;
00911 const double * lower = getColLower() ;
00912 const double * upper = getColUpper() ;
00913 for (i = 0 ; i < numberColumns ; i++)
00914 { lowerBefore[i]= lower[i] ;
00915 upperBefore[i]= upper[i] ; }
00916 bool deleteNode ;
00917 if (node->branch())
00918 { branchingTree.push(node) ;
00919 deleteNode = false ;
00920 # ifdef CHECK_NODE
00921 printf("Node %x pushed back on tree - %d left, %d count\n",node,
00922 nodeInfo->numberBranchesLeft(),
00923 nodeInfo->numberPointingToThis()) ;
00924 # endif
00925 }
00926 else
00927 { deleteNode = true ; }
00928
00929 # ifdef SBB_DEBUG
00930
00931
00932
00933
00934
00935
00936 const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
00937 if (debugger)
00938 { if(debugger->onOptimalPath(*solver_))
00939 printf("On optimal path\n") ;
00940 else
00941 printf("Not on optimal path\n") ; }
00942 # endif
00943
00944
00945
00946
00947
00948 cuts = OsiCuts() ;
00949 currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_ ;
00950 feasible = solveWithCuts(cuts,maximumCutPasses_,node,
00951 numberOldActiveCuts,numberNewCuts,
00952 maximumWhich,whichGenerator) ;
00953
00954
00955
00956 double totalTime = CoinCpuTime()-time1 ;
00957 if (numberNodes_ < intParam_[SbbMaxNumNode] &&
00958 numberSolutions_ < intParam_[SbbMaxNumSol] &&
00959 totalTime < dblParam_[SbbMaximumSeconds] &&
00960 !stoppedOnGap)
00961 {
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976 if (feasible)
00977 { newNode = new SbbNode ;
00978 newNode->setObjectiveValue(direction*solver_->getObjValue()) ;
00979 anyAction =-1 ;
00980 resolved = false ;
00981 while (anyAction == -1)
00982 { anyAction = newNode->chooseBranch(this,node) ;
00983 if (anyAction == -1)
00984 { feasible = resolve() ;
00985 resolved = true ;
00986 if (feasible)
00987 { newNode->setObjectiveValue(direction*
00988 solver_->getObjValue()) ; }
00989 else
00990 { anyAction = -2 ; } } }
00991 if (anyAction >= 0)
00992 { if (resolved)
00993 { bool needValidSolution = (newNode->variable() < 0) ;
00994 takeOffCuts(cuts,whichGenerator,numberOldActiveCuts,
00995 numberNewCuts,needValidSolution) ;
00996 # ifdef CHECK_CUT_COUNTS
00997 { printf("Number of rows after chooseBranch fix (node)"
00998 "(active only) %d\n",
00999 numberRowsAtContinuous_+numberNewCuts+
01000 numberOldActiveCuts) ;
01001 const CoinWarmStartBasis* debugws =
01002 dynamic_cast<const CoinWarmStartBasis*>
01003 (solver_->getWarmStart()) ;
01004 debugws->print() ;
01005 delete debugws ; }
01006 # endif
01007 }
01008 newNode->createInfo(this,node,lastws,lowerBefore,upperBefore,
01009 numberOldActiveCuts,numberNewCuts) ;
01010 if (newNode->numberUnsatisfied())
01011 newNode->nodeInfo()->addCuts(cuts,newNode->numberBranches(),
01012 whichGenerator) ; } }
01013 else
01014 { anyAction = -2 ; }
01015
01016 if ( anyAction >=0 ) {
01017 assert (newNode);
01018 if (newNode->objectiveValue() > getCutoff())
01019 anyAction = -2;
01020 }
01021
01022
01023
01024
01025
01026
01027 if (anyAction == -2)
01028 { delete newNode ;
01029 newNode = NULL ;
01030 for (i = 0 ; i < currentNumberCuts_ ; i++)
01031 { if (addedCuts_[i])
01032 { if (!addedCuts_[i]->decrement(1))
01033 delete addedCuts_[i] ; } } }
01034 else
01035 { nodeInfo->increment() ; }
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050 assert (!newNode || newNode->objectiveValue() <= getCutoff()) ;
01051 if (newNode)
01052 { if (newNode->variable() >= 0)
01053 { handler_->message(SBB_BRANCH,messages_)
01054 << numberNodes_<< newNode->objectiveValue()
01055 << newNode->numberUnsatisfied()<< newNode->depth()
01056 << integerVariable_[newNode->variable()]
01057 << newNode->variable()
01058 << CoinMessageEol ;
01059
01060 int numberLeft = newNode->numberBranches() ;
01061 for (i = 0;i < currentNumberCuts_;i++)
01062 { if (addedCuts_[i])
01063 {
01064 # ifdef CHECK_CUT_COUNTS
01065 printf("Count on cut %x increased by %d\n",addedCuts_[i],
01066 numberLeft-1) ;
01067 # endif
01068 addedCuts_[i]->increment(numberLeft-1) ; } }
01069
01070 double estValue = 1.0e100 ;
01071 int found = -1 ;
01072 solver_->resolve() ;
01073 double * newSolution = new double [numberColumns] ;
01074 double heurValue = getCutoff() ;
01075 int iHeur ;
01076 for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++)
01077 { double saveValue = heurValue ;
01078 int ifSol = heuristic_[iHeur]->solution(heurValue,newSolution) ;
01079 if (ifSol > 0)
01080 { found = iHeur ; }
01081 else
01082 if (ifSol < 0)
01083 { estValue = min(heurValue,estValue) ;
01084 heurValue = saveValue ; } }
01085 if (found >= 0)
01086 { setBestSolution(SBB_ROUNDING,heurValue,newSolution) ; }
01087 delete [] newSolution ;
01088 newNode->setGuessedObjectiveValue(estValue) ;
01089 branchingTree.push(newNode) ;
01090 # ifdef CHECK_NODE
01091 printf("Node %x pushed on tree c\n",newNode) ;
01092 # endif
01093 }
01094 else
01095 { for (i = 0 ; i < currentNumberCuts_ ; i++)
01096 { if (addedCuts_[i])
01097 { if (!addedCuts_[i]->decrement(1))
01098 delete addedCuts_[i] ; } }
01099 double objectiveValue = newNode->objectiveValue();
01100 setBestSolution(SBB_SOLUTION,objectiveValue,
01101 solver_->getColSolution()) ;
01102 assert(nodeInfo->numberPointingToThis() <= 2) ;
01103
01104 nodeInfo->increment();
01105 delete newNode ;
01106 nodeInfo->decrement() ; } }
01107
01108
01109
01110
01111 if (deleteNode)
01112 delete node ; }
01113
01114
01115
01116
01117
01118 else
01119 { branchingTree.cleanTree(this,-DBL_MAX) ;
01120 if (stoppedOnGap)
01121 { messageHandler()->message(SBB_GAP,messages())
01122 << bestObjective_-bestValue
01123 << dblParam_[SbbAllowableGap]
01124 << CoinMessageEol ;
01125 status_ = 0 ; }
01126 else
01127 if (isNodeLimitReached())
01128 { handler_->message(SBB_MAXNODES,messages_) << CoinMessageEol ;
01129 status_ = 1 ; }
01130 else
01131 if (totalTime >= dblParam_[SbbMaximumSeconds])
01132 { handler_->message(SBB_MAXTIME,messages_) << CoinMessageEol ; }
01133 else
01134 { handler_->message(SBB_MAXSOLS,messages_) << CoinMessageEol ;
01135 status_ = 1 ; }
01136 break ; }
01137
01138
01139
01140
01141
01142
01143 int numberToDelete = getNumRows()-numberRowsAtContinuous_ ;
01144 if (numberToDelete)
01145 { int * delRows = new int[numberToDelete] ;
01146 int i ;
01147 for (i = 0 ; i < numberToDelete ; i++)
01148 { delRows[i] = i+numberRowsAtContinuous_ ; }
01149 solver_->deleteRows(numberToDelete,delRows) ;
01150 delete [] delRows ; } }
01151
01152
01153
01154 else
01155 { delete node ; } }
01156
01157
01158
01159
01160 if (!status_)
01161 handler_->message(SBB_END_GOOD,messages_)
01162 << bestObjective_ << numberIterations_ << numberNodes_
01163 << CoinMessageEol ;
01164 else
01165 handler_->message(SBB_END,messages_)
01166 << numberIterations_ << numberNodes_ << CoinMessageEol ;
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182 if (bestSolution_)
01183 { setCutoff(1.0e50) ;
01184 setBestSolution(SBB_SOLUTION,bestObjective_,bestSolution_,true) ;
01185 continuousSolver_->resolve() ;
01186 if (!continuousSolver_->isProvenOptimal())
01187 { continuousSolver_->messageHandler()->setLogLevel(2) ;
01188 continuousSolver_->initialSolve() ; }
01189 delete solver_ ;
01190 solver_ = continuousSolver_ ;
01191 continuousSolver_ = NULL ; }
01192
01193
01194
01195 delete lastws ;
01196 delete [] whichGenerator ;
01197 delete [] lowerBefore ;
01198 delete [] upperBefore ;
01199 delete [] walkback_ ;
01200 walkback_ = NULL ;
01201 delete [] addedCuts_ ;
01202 addedCuts_ = NULL ;
01203 if (continuousSolver_)
01204 { delete continuousSolver_ ;
01205 continuousSolver_ = NULL ; }
01206
01207
01208
01209 globalCuts_= OsiCuts() ;
01210
01211 return ; }
01212
01213
01214
01215
01216 void
01217 SbbModel::initialSolve()
01218 {
01219 assert (solver_);
01220 solver_->initialSolve();
01221 }
01222
01232 CoinWarmStartBasis *SbbModel::getEmptyBasis (int ns, int na) const
01233
01234 { CoinWarmStartBasis *emptyBasis ;
01235
01236
01237
01238 if (emptyWarmStart_ == 0)
01239 { if (solver_ == 0)
01240 { throw CoinError("Cannot construct basis without solver!",
01241 "getEmptyBasis","SbbModel") ; }
01242 emptyBasis =
01243 dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
01244 if (emptyBasis == 0)
01245 { throw CoinError(
01246 "Solver does not appear to use a basis-oriented warm start.",
01247 "getEmptyBasis","SbbModel") ; }
01248 emptyBasis->setSize(0,0) ;
01249 emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ; }
01250
01251
01252
01253 emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
01254 assert(emptyBasis) ;
01255 if (ns != 0 || na != 0) emptyBasis->setSize(ns,na) ;
01256
01257 return (emptyBasis) ; }
01258
01259
01264 SbbModel::SbbModel()
01265
01266 :
01267 solver_(NULL),
01268 ourSolver_(true),
01269 continuousSolver_(NULL),
01270 defaultHandler_(true),
01271 emptyWarmStart_(NULL),
01272 basis_(NULL),
01273 bestObjective_(DBL_MAX),
01274 bestSolution_(NULL),
01275 currentSolution_(NULL),
01276 minimumDrop_(1.0e-4),
01277 numberSolutions_(0),
01278 forcePriority_(-1),
01279 numberHeuristicSolutions_(0),
01280 numberNodes_(0),
01281 numberIterations_(0),
01282 status_(0),
01283 numberIntegers_(0),
01284 numberRowsAtContinuous_(0),
01285 maximumNumberCuts_(0),
01286 currentNumberCuts_(0),
01287 maximumDepth_(0),
01288 walkback_(NULL),
01289 addedCuts_(NULL),
01290 integerVariable_(NULL),
01291 strategy_(0),
01292 numberStrong_(5),
01293 printFrequency_(0),
01294 numberCutGenerators_(0),
01295 generator_(NULL),
01296 numberHeuristics_(0),
01297 heuristic_(NULL),
01298 numberObjects_(0),
01299 object_(NULL),
01300 originalColumns_(NULL),
01301 priority_(NULL),
01302 howOftenGlobalScan_(1),
01303 continuousObjective_(DBL_MAX),
01304 continuousInfeasibilities_(INT_MAX),
01305 maximumCutPassesAtRoot_(20),
01306 maximumCutPasses_(10)
01307 {
01308 intParam_[SbbMaxNumNode] = 9999999;
01309 intParam_[SbbMaxNumSol] = 9999999;
01310 intParam_[SbbFathomDiscipline] = 0;
01311
01312 dblParam_[SbbIntegerTolerance] = 1e-6;
01313 dblParam_[SbbInfeasibilityWeight] = 0.0;
01314 dblParam_[SbbCutoffIncrement] = 1e-5;
01315 dblParam_[SbbAllowableGap] = 1.0e-10;
01316 dblParam_[SbbMaximumSeconds] = 1.0e100;
01317 nodeCompare_=NULL;
01318 branchingMethod_=NULL;
01319 handler_ = new CoinMessageHandler();
01320 handler_->setLogLevel(2);
01321 messages_ = SbbMessage();
01322 }
01323
01329 SbbModel::SbbModel(const OsiSolverInterface &rhs)
01330 :
01331 continuousSolver_(NULL),
01332 defaultHandler_(true),
01333 emptyWarmStart_(NULL),
01334 basis_(NULL) ,
01335 bestObjective_(DBL_MAX),
01336 minimumDrop_(1.0e-4),
01337 numberSolutions_(0),
01338 forcePriority_(-1),
01339 numberHeuristicSolutions_(0),
01340 numberNodes_(0),
01341 numberIterations_(0),
01342 status_(0),
01343 numberRowsAtContinuous_(0),
01344 maximumNumberCuts_(0),
01345 currentNumberCuts_(0),
01346 maximumDepth_(0),
01347 walkback_(NULL),
01348 addedCuts_(NULL),
01349 strategy_(0),
01350 numberStrong_(5),
01351 printFrequency_(0),
01352 numberCutGenerators_(0),
01353 generator_(NULL),
01354 numberHeuristics_(0),
01355 heuristic_(NULL),
01356 numberObjects_(0),
01357 object_(NULL),
01358 originalColumns_(NULL),
01359 priority_(NULL),
01360 howOftenGlobalScan_(-1),
01361 continuousObjective_(DBL_MAX),
01362 continuousInfeasibilities_(INT_MAX),
01363 maximumCutPassesAtRoot_(20),
01364 maximumCutPasses_(10)
01365 {
01366 intParam_[SbbMaxNumNode] = 9999999;
01367 intParam_[SbbMaxNumSol] = 9999999;
01368 intParam_[SbbFathomDiscipline] = 0;
01369
01370 dblParam_[SbbIntegerTolerance] = 1e-6;
01371 dblParam_[SbbInfeasibilityWeight] = 0.0;
01372 dblParam_[SbbCutoffIncrement] = 1e-5;
01373 dblParam_[SbbAllowableGap] = 1.0e-10;
01374 dblParam_[SbbMaximumSeconds] = 1.0e100;
01375
01376 nodeCompare_=NULL;
01377 branchingMethod_=NULL;
01378 handler_ = new CoinMessageHandler();
01379 handler_->setLogLevel(2);
01380 messages_ = SbbMessage();
01381 solver_ = rhs.clone();
01382 ourSolver_ = true ;
01383
01384
01385 bestSolution_ = NULL;
01386 numberIntegers_=0;
01387 int numberColumns = solver_->getNumCols();
01388 int iColumn;
01389 if (numberColumns) {
01390
01391 currentSolution_ = new double[numberColumns];
01392 for (iColumn=0;iColumn<numberColumns;iColumn++) {
01393 if( solver_->isInteger(iColumn))
01394 numberIntegers_++;
01395 }
01396 } else {
01397
01398 currentSolution_=NULL;
01399 }
01400 if (numberIntegers_) {
01401 integerVariable_ = new int [numberIntegers_];
01402 numberIntegers_=0;
01403 for (iColumn=0;iColumn<numberColumns;iColumn++) {
01404 if( solver_->isInteger(iColumn))
01405 integerVariable_[numberIntegers_++]=iColumn;
01406 }
01407 } else {
01408 integerVariable_ = NULL;
01409 }
01410 }
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430 void
01431 SbbModel::assignSolver(OsiSolverInterface *&solver)
01432
01433 {
01434
01435 if (solver_)
01436 solver->messageHandler()->setLogLevel(solver_->messageHandler()->logLevel()) ;
01437
01438 if (ourSolver_) delete solver_ ;
01439 solver_ = solver;
01440 solver = NULL ;
01441 ourSolver_ = true ;
01442
01443
01444
01445 if (basis_)
01446 { delete basis_ ;
01447 basis_ = 0 ; }
01448 if (emptyWarmStart_)
01449 { delete emptyWarmStart_ ;
01450 emptyWarmStart_ = 0 ; }
01451
01452
01453
01454 numberIntegers_=0;
01455 int numberColumns = solver_->getNumCols();
01456 int iColumn;
01457 for (iColumn=0;iColumn<numberColumns;iColumn++) {
01458 if( solver_->isInteger(iColumn))
01459 numberIntegers_++;
01460 }
01461 if (numberIntegers_) {
01462 integerVariable_ = new int [numberIntegers_];
01463 numberIntegers_=0;
01464 for (iColumn=0;iColumn<numberColumns;iColumn++) {
01465 if( solver_->isInteger(iColumn))
01466 integerVariable_[numberIntegers_++]=iColumn;
01467 }
01468 } else {
01469 integerVariable_ = NULL;
01470 }
01471
01472 return ;
01473 }
01474
01475
01476
01477 SbbModel::SbbModel(const SbbModel & rhs)
01478 :
01479 continuousSolver_(NULL),
01480 defaultHandler_(rhs.defaultHandler_),
01481 emptyWarmStart_(NULL),
01482 basis_(NULL),
01483 bestObjective_(rhs.bestObjective_),
01484 minimumDrop_(rhs.minimumDrop_),
01485 numberSolutions_(rhs.numberSolutions_),
01486 forcePriority_(rhs.forcePriority_),
01487 numberHeuristicSolutions_(rhs.numberHeuristicSolutions_),
01488 numberNodes_(rhs.numberNodes_),
01489 numberIterations_(rhs.numberIterations_),
01490 status_(rhs.status_),
01491 strategy_(rhs.strategy_),
01492 numberStrong_(rhs.numberStrong_),
01493 printFrequency_(rhs.printFrequency_),
01494 howOftenGlobalScan_(rhs.howOftenGlobalScan_),
01495 continuousObjective_(rhs.continuousObjective_),
01496 continuousInfeasibilities_(rhs.continuousInfeasibilities_),
01497 maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_),
01498 maximumCutPasses_( rhs.maximumCutPasses_)
01499 {
01500 intParam_[SbbMaxNumNode] = rhs.intParam_[SbbMaxNumNode];
01501 intParam_[SbbMaxNumSol] = rhs.intParam_[SbbMaxNumSol];
01502 intParam_[SbbFathomDiscipline] = rhs.intParam_[SbbFathomDiscipline];
01503 dblParam_[SbbIntegerTolerance] = rhs.dblParam_[SbbIntegerTolerance];
01504 dblParam_[SbbInfeasibilityWeight] = rhs.dblParam_[SbbInfeasibilityWeight];
01505 dblParam_[SbbCutoffIncrement] = rhs.dblParam_[SbbCutoffIncrement];
01506 dblParam_[SbbAllowableGap] = rhs.dblParam_[SbbAllowableGap];
01507 dblParam_[SbbMaximumSeconds] = rhs.dblParam_[SbbMaximumSeconds];
01508 if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
01509 if (rhs.basis_) basis_ =
01510 dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ;
01511 if (defaultHandler_) {
01512 handler_ = new CoinMessageHandler();
01513 handler_->setLogLevel(2);
01514 } else {
01515 handler_ = rhs.handler_;
01516 }
01517 messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
01518 numberCutGenerators_ = rhs.numberCutGenerators_;
01519 if (numberCutGenerators_) {
01520 generator_ = new SbbCutGenerator * [numberCutGenerators_];
01521 int i;
01522 for (i=0;i<numberCutGenerators_;i++) {
01523 generator_[i]=new SbbCutGenerator(*rhs.generator_[i]);
01524 }
01525 } else {
01526 generator_=NULL;
01527 }
01528 globalCuts_ = rhs.globalCuts_;
01529 numberHeuristics_ = rhs.numberHeuristics_;
01530 if (numberHeuristics_) {
01531 heuristic_ = new SbbHeuristic * [numberHeuristics_];
01532 memcpy(heuristic_,rhs.heuristic_,
01533 numberHeuristics_*sizeof(SbbHeuristic *));
01534 } else {
01535 heuristic_=NULL;
01536 }
01537 numberObjects_=rhs.numberObjects_;
01538 if (numberObjects_) {
01539 object_ = new SbbObject * [numberObjects_];
01540 int i;
01541 for (i=0;i<numberObjects_;i++)
01542 object_[i]=(rhs.object_[i])->clone();
01543 } else {
01544 object_=NULL;
01545 }
01546 if (rhs.priority_) {
01547 priority_= new int [numberObjects_];
01548 memcpy(priority_,rhs.priority_,numberObjects_*sizeof(int));
01549 } else {
01550 priority_=NULL;
01551 }
01552 if (rhs.originalColumns_) {
01553 int numberColumns = solver_->getNumCols();
01554 originalColumns_= new int [numberColumns];
01555 memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
01556 } else {
01557 originalColumns_=NULL;
01558 }
01559 nodeCompare_=rhs.nodeCompare_;
01560 branchingMethod_=rhs.branchingMethod_;
01561 messages_ = rhs.messages_;
01562 solver_ = rhs.solver_->clone();
01563 ourSolver_ = true ;
01564 messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
01565 numberIntegers_=rhs.numberIntegers_;
01566 if (numberIntegers_) {
01567 integerVariable_ = new int [numberIntegers_];
01568 memcpy(integerVariable_,rhs.integerVariable_,numberIntegers_*sizeof(int));
01569 } else {
01570 integerVariable_ = NULL;
01571 }
01572 if (rhs.bestSolution_) {
01573 int numberColumns = solver_->getNumCols();
01574 bestSolution_ = new double[numberColumns];
01575 memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
01576 } else {
01577 bestSolution_=NULL;
01578 }
01579 if (rhs.currentSolution_) {
01580 int numberColumns = solver_->getNumCols();
01581 currentSolution_ = new double[numberColumns];
01582 memcpy(currentSolution_,rhs.currentSolution_,numberColumns*sizeof(double));
01583 } else {
01584 currentSolution_=NULL;
01585 }
01586 numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
01587 maximumNumberCuts_=rhs.maximumNumberCuts_;
01588 currentNumberCuts_=rhs.currentNumberCuts_;
01589 maximumDepth_= rhs.maximumDepth_;
01590
01591 if (maximumNumberCuts_) {
01592 addedCuts_ = new SbbCountRowCut * [maximumNumberCuts_];
01593 } else {
01594 addedCuts_ = NULL;
01595 }
01596 if (maximumDepth_)
01597 walkback_ = new SbbNodeInfo * [maximumDepth_];
01598 else
01599 walkback_ = NULL;
01600 synchronizeModel();
01601 }
01602
01603
01604 SbbModel &
01605 SbbModel::operator=(const SbbModel& rhs)
01606 {
01607 if (this!=&rhs) {
01608
01609 if (defaultHandler_)
01610 { delete handler_;
01611 handler_ = NULL; }
01612 defaultHandler_ = rhs.defaultHandler_;
01613 if (defaultHandler_)
01614 { handler_ = new CoinMessageHandler();
01615 handler_->setLogLevel(2); }
01616 else
01617 { handler_ = rhs.handler_; }
01618 messages_ = rhs.messages_;
01619 messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
01620
01621 delete solver_;
01622 if (rhs.solver_)
01623 { solver_ = rhs.solver_->clone() ; }
01624 else
01625 { solver_ = 0 ; }
01626 ourSolver_ = true ;
01627 delete continuousSolver_ ;
01628 if (rhs.continuousSolver_)
01629 { solver_ = rhs.continuousSolver_->clone() ; }
01630 else
01631 { continuousSolver_ = 0 ; }
01632
01633 delete emptyWarmStart_ ;
01634 if (rhs.emptyWarmStart_)
01635 { emptyWarmStart_ = rhs.emptyWarmStart_->clone() ; }
01636 else
01637 { emptyWarmStart_ = 0 ; }
01638 delete basis_ ;
01639 if (rhs.basis_)
01640 { basis_ = dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ; }
01641 else
01642 { basis_ = 0 ; }
01643
01644 bestObjective_ = rhs.bestObjective_;
01645 delete [] bestSolution_;
01646 if (rhs.bestSolution_) {
01647 int numberColumns = rhs.getNumCols();
01648 bestSolution_ = new double[numberColumns];
01649 memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
01650 } else {
01651 bestSolution_=NULL;
01652 }
01653 if (rhs.currentSolution_) {
01654 int numberColumns = rhs.getNumCols();
01655 currentSolution_ = new double[numberColumns];
01656 memcpy(currentSolution_,rhs.currentSolution_,numberColumns*sizeof(double));
01657 } else {
01658 currentSolution_=NULL;
01659 }
01660 minimumDrop_ = rhs.minimumDrop_;
01661 numberSolutions_=rhs.numberSolutions_;
01662 forcePriority_=rhs.forcePriority_;
01663 numberHeuristicSolutions_=rhs.numberHeuristicSolutions_;
01664 numberNodes_ = rhs.numberNodes_;
01665 numberIterations_ = rhs.numberIterations_;
01666 status_ = rhs.status_;
01667 strategy_ = rhs.strategy_;
01668 numberStrong_ = rhs.numberStrong_;
01669 printFrequency_ = rhs.printFrequency_;
01670 howOftenGlobalScan_=rhs.howOftenGlobalScan_;
01671 continuousObjective_=rhs.continuousObjective_;
01672 continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
01673 maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
01674 maximumCutPasses_ = rhs.maximumCutPasses_;
01675 intParam_[SbbMaxNumNode] = rhs.intParam_[SbbMaxNumNode];
01676 intParam_[SbbMaxNumSol] = rhs.intParam_[SbbMaxNumSol];
01677 intParam_[SbbFathomDiscipline] = rhs.intParam_[SbbFathomDiscipline];
01678 dblParam_[SbbIntegerTolerance] = rhs.dblParam_[SbbIntegerTolerance];
01679 dblParam_[SbbInfeasibilityWeight] = rhs.dblParam_[SbbInfeasibilityWeight];
01680 dblParam_[SbbCutoffIncrement] = rhs.dblParam_[SbbCutoffIncrement];
01681 dblParam_[SbbAllowableGap] = rhs.dblParam_[SbbAllowableGap];
01682 dblParam_[SbbMaximumSeconds] = rhs.dblParam_[SbbMaximumSeconds];
01683 globalCuts_ = rhs.globalCuts_;
01684 int i;
01685 for (i=0;i<numberCutGenerators_;i++)
01686 delete generator_[i];
01687 delete [] generator_;
01688 delete [] heuristic_;
01689 numberCutGenerators_ = rhs.numberCutGenerators_;
01690 if (numberCutGenerators_) {
01691 generator_ = new SbbCutGenerator * [numberCutGenerators_];
01692 int i;
01693 for (i=0;i<numberCutGenerators_;i++) {
01694 generator_[i]=new SbbCutGenerator(*rhs.generator_[i]);
01695 }
01696 } else {
01697 generator_=NULL;
01698 }
01699 numberHeuristics_ = rhs.numberHeuristics_;
01700 if (numberHeuristics_) {
01701 heuristic_ = new SbbHeuristic * [numberHeuristics_];
01702 memcpy(heuristic_,rhs.heuristic_,
01703 numberHeuristics_*sizeof(SbbHeuristic *));
01704 } else {
01705 heuristic_=NULL;
01706 }
01707 for (i=0;i<numberObjects_;i++)
01708 delete object_[i];
01709 delete [] object_;
01710 numberObjects_=rhs.numberObjects_;
01711 if (numberObjects_) {
01712 object_ = new SbbObject * [numberObjects_];
01713 int i;
01714 for (i=0;i<numberObjects_;i++)
01715 object_[i]=(rhs.object_[i])->clone();
01716 } else {
01717 object_=NULL;
01718 }
01719 delete [] priority_;
01720 if (rhs.priority_) {
01721 priority_= new int [numberObjects_];
01722 memcpy(priority_,rhs.priority_,numberObjects_*sizeof(int));
01723 } else {
01724 priority_=NULL;
01725 }
01726 delete [] originalColumns_;
01727 if (rhs.originalColumns_) {
01728 int numberColumns = rhs.getNumCols();
01729 originalColumns_= new int [numberColumns];
01730 memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
01731 } else {
01732 originalColumns_=NULL;
01733 }
01734 nodeCompare_=rhs.nodeCompare_;
01735 branchingMethod_=rhs.branchingMethod_;
01736
01737 delete [] integerVariable_;
01738 numberIntegers_=rhs.numberIntegers_;
01739 if (numberIntegers_) {
01740 integerVariable_ = new int [numberIntegers_];
01741 memcpy(integerVariable_,rhs.integerVariable_,
01742 numberIntegers_*sizeof(int));
01743 } else {
01744 integerVariable_ = NULL;
01745 }
01746 numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
01747 maximumNumberCuts_=rhs.maximumNumberCuts_;
01748 currentNumberCuts_=rhs.currentNumberCuts_;
01749 maximumDepth_= rhs.maximumDepth_;
01750 delete [] addedCuts_;
01751 delete [] walkback_;
01752
01753 if (maximumNumberCuts_) {
01754 addedCuts_ = new SbbCountRowCut * [maximumNumberCuts_];
01755 } else {
01756 addedCuts_ = NULL;
01757 }
01758 if (maximumDepth_)
01759 walkback_ = new SbbNodeInfo * [maximumDepth_];
01760 else
01761 walkback_ = NULL;
01762 synchronizeModel();
01763 }
01764 return *this;
01765 }
01766
01767
01768 SbbModel::~SbbModel ()
01769 {
01770 if (defaultHandler_) {
01771 delete handler_;
01772 handler_ = NULL;
01773 }
01774 if (ourSolver_) delete solver_;
01775 gutsOfDestructor();
01776 }
01777
01778 void
01779 SbbModel::gutsOfDestructor()
01780 {
01781 delete emptyWarmStart_ ;
01782 emptyWarmStart_ =NULL;
01783 delete basis_ ;
01784 basis_ =NULL;
01785 delete continuousSolver_;
01786 continuousSolver_=NULL;
01787 delete [] bestSolution_;
01788 bestSolution_=NULL;
01789 delete [] currentSolution_;
01790 currentSolution_=NULL;
01791 delete [] integerVariable_;
01792 integerVariable_=NULL;
01793 int i;
01794 for (i=0;i<numberCutGenerators_;i++)
01795 delete generator_[i];
01796 delete [] generator_;
01797 generator_=NULL;
01798 delete [] heuristic_;
01799 heuristic_=NULL;
01800 delete [] addedCuts_;
01801 addedCuts_=NULL;
01802 delete [] walkback_;
01803 walkback_=NULL;
01804 for (i=0;i<numberObjects_;i++)
01805 delete object_[i];
01806 delete [] object_;
01807 object_=NULL;
01808 delete [] priority_;
01809 priority_=NULL;
01810 delete [] originalColumns_;
01811 originalColumns_=NULL;
01812 }
01813
01814 bool
01815 SbbModel::isAbandoned() const
01816 {
01817 return status_ == 2;
01818 }
01819
01820 bool
01821 SbbModel::isProvenOptimal() const
01822 {
01823 if (!status_ && bestObjective_<1.0e30)
01824 return true;
01825 else
01826 return false;
01827 }
01828
01829 bool
01830 SbbModel::isProvenInfeasible() const
01831 {
01832 if (!status_ && bestObjective_==1.0e30)
01833 return true;
01834 else
01835 return false;
01836 }
01837
01838 bool
01839 SbbModel::isNodeLimitReached() const
01840 {
01841 return numberNodes_ >= intParam_[SbbMaxNumNode];
01842 }
01843
01844 bool
01845 SbbModel::isSolutionLimitReached() const
01846 {
01847 return numberSolutions_ >= intParam_[SbbMaxNumSol];
01848 }
01849
01850 void
01851 SbbModel::newLanguage(CoinMessages::Language language)
01852 {
01853 messages_ = SbbMessage(language);
01854 }
01855 void
01856 SbbModel::setNumberStrong(int number)
01857 {
01858 if (number<0)
01859 numberStrong_=0;
01860 else
01861 numberStrong_=number;
01862 }
01863
01864 void
01865 SbbModel::addCutGenerator(CglCutGenerator * generator,
01866 int howOften, const char * name,
01867 bool normal, bool atSolution,
01868 bool whenInfeasible)
01869 {
01870 SbbCutGenerator ** temp = generator_;
01871 generator_ = new SbbCutGenerator * [numberCutGenerators_+1];
01872 memcpy(generator_,temp,numberCutGenerators_*sizeof(SbbCutGenerator *));
01873 delete[] temp ;
01874 generator_[numberCutGenerators_++]=
01875 new SbbCutGenerator(this,generator, howOften, name,
01876 normal,atSolution,whenInfeasible);
01877
01878 }
01879
01880 void
01881 SbbModel::addHeuristic(SbbHeuristic * generator)
01882 {
01883 SbbHeuristic ** temp = heuristic_;
01884 heuristic_ = new SbbHeuristic * [numberHeuristics_+1];
01885 memcpy(heuristic_,temp,numberHeuristics_*sizeof(SbbHeuristic *));
01886 delete [] temp;
01887 heuristic_[numberHeuristics_++]=generator;
01888 }
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909 void SbbModel::addCuts1 (SbbNode * node, CoinWarmStartBasis *&lastws)
01910 { int i;
01911 int nNode=0;
01912 int numberColumns = getNumCols();
01913 SbbNodeInfo * nodeInfo = node->nodeInfo();
01914
01915
01916
01917
01918
01919
01920 int currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_;
01921 int *which = new int[currentNumberCuts];
01922 for (i = 0 ; i < currentNumberCuts ; i++)
01923 which[i] = i+numberRowsAtContinuous_;
01924 solver_->deleteRows(currentNumberCuts,which);
01925 delete [] which;
01926
01927
01928
01929
01930
01931
01932
01933 currentNumberCuts = 0;
01934 while (nodeInfo) {
01935
01936 walkback_[nNode++]=nodeInfo;
01937 currentNumberCuts += nodeInfo->numberCuts() ;
01938 nodeInfo = nodeInfo->parent() ;
01939 if (nNode==maximumDepth_) {
01940 maximumDepth_ *= 2;
01941 SbbNodeInfo ** temp = new SbbNodeInfo * [maximumDepth_];
01942 for (i=0;i<nNode;i++)
01943 temp[i] = walkback_[i];
01944 delete [] walkback_;
01945 walkback_ = temp;
01946 }
01947 }
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960 currentNumberCuts_=currentNumberCuts;
01961 if (currentNumberCuts >= maximumNumberCuts_) {
01962 maximumNumberCuts_ = currentNumberCuts;
01963 delete [] addedCuts_;
01964 addedCuts_ = new SbbCountRowCut * [maximumNumberCuts_];
01965 }
01966 lastws->setSize(numberColumns,numberRowsAtContinuous_+currentNumberCuts);
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977 currentNumberCuts=0;
01978 while (nNode) {
01979 --nNode;
01980 walkback_[nNode]->applyToModel(this,lastws,addedCuts_,currentNumberCuts);
01981 }
01982 }
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997 int SbbModel::addCuts (SbbNode *node, CoinWarmStartBasis *&lastws)
01998 {
01999
02000
02001
02002
02003 addCuts1(node,lastws);
02004 int i;
02005 int numberColumns = getNumCols();
02006 SbbNodeInfo * nodeInfo = node->nodeInfo();
02007 double cutoff = getCutoff() ;
02008 int currentNumberCuts=currentNumberCuts_;
02009
02010
02011
02012
02013 if (node->objectiveValue() < cutoff)
02014 { int numberToAdd = 0;
02015 OsiRowCut * addCuts;
02016 if (currentNumberCuts == 0)
02017 addCuts = NULL;
02018 else
02019 addCuts = new OsiRowCut[currentNumberCuts];
02020 # ifdef CHECK_CUT_COUNTS
02021 printf("addCuts: expanded basis; rows %d+%d\n",
02022 numberRowsAtContinuous_,currentNumberCuts);
02023 lastws->print();
02024 # endif
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041 if (basis_) delete basis_ ;
02042 basis_= dynamic_cast<CoinWarmStartBasis *>(lastws->clone()) ;
02043 for (i=0;i<currentNumberCuts;i++) {
02044 CoinWarmStartBasis::Status status =
02045 lastws->getArtifStatus(i+numberRowsAtContinuous_);
02046 if (status != CoinWarmStartBasis::basic&&addedCuts_[i]) {
02047 # ifdef CHECK_CUT_COUNTS
02048 printf("Using cut %d %x as row %d\n",i,addedCuts_[i],
02049 numberRowsAtContinuous_+numberToAdd);
02050 # endif
02051 lastws->setArtifStatus(numberToAdd+numberRowsAtContinuous_,status);
02052 addCuts[numberToAdd++] = OsiRowCut(*addedCuts_[i]);
02053 } else {
02054 # ifdef CHECK_CUT_COUNTS
02055 printf("Dropping cut %d %x\n",i,addedCuts_[i]);
02056 # endif
02057 addedCuts_[i]=NULL;
02058 }
02059 }
02060 int numberRowsNow=numberRowsAtContinuous_+numberToAdd;
02061 lastws->resize(numberRowsNow,numberColumns);
02062 #ifdef FULL_DEBUG
02063 printf("addCuts: stripped basis; rows %d + %d\n",
02064 numberRowsAtContinuous_,numberToAdd);
02065 lastws->print();
02066 #endif
02067
02068
02069
02070 solver_->applyRowCuts(numberToAdd,addCuts);
02071 solver_->setWarmStart(lastws);
02072
02073
02074
02075 delete lastws ;
02076 lastws = basis_ ;
02077 basis_ = 0 ;
02078
02079 #if 0
02080 if ((numberNodes_%printFrequency_)==0) {
02081 printf("Objective %g, depth %d, unsatisfied %d\n",
02082 node->objectiveValue(),
02083 node->depth(),node->numberUnsatisfied());
02084 }
02085 #endif
02086
02087
02088
02089 delete [] addCuts;
02090 numberNodes_++;
02091 return 0;
02092 }
02093
02094
02095
02096
02097
02098
02099 else
02100 { int i;
02101 int numberLeft = nodeInfo->numberBranchesLeft();
02102 for (i = 0 ; i < currentNumberCuts ; i++)
02103 { if (addedCuts_[i])
02104 { if (!addedCuts_[i]->decrement(numberLeft))
02105 { delete addedCuts_[i];
02106 addedCuts_[i] = NULL; } } }
02107 return 1 ; }
02108 }
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118 void SbbModel::reducedCostFix ()
02119
02120 { double cutoff = getCutoff() ;
02121 double direction = solver_->getObjSense() ;
02122 double gap = cutoff - solver_->getObjValue()*direction ;
02123 double integerTolerance = getDblParam(SbbIntegerTolerance) ;
02124
02125 const double *lower = solver_->getColLower() ;
02126 const double *upper = solver_->getColUpper() ;
02127 const double *solution = solver_->getColSolution() ;
02128 const double *reducedCost = solver_->getReducedCost() ;
02129
02130 int numberFixed = 0 ;
02131 for (int i = 0 ; i < numberIntegers_ ; i++)
02132 { int iColumn = integerVariable_[i] ;
02133 double djValue = direction*reducedCost[iColumn] ;
02134 if (upper[iColumn]-lower[iColumn] > integerTolerance)
02135 { if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap)
02136 { solver_->setColUpper(iColumn,lower[iColumn]) ;
02137 numberFixed++ ; }
02138 else
02139 if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap)
02140 { solver_->setColLower(iColumn,upper[iColumn]) ;
02141 numberFixed++ ; } } }
02142
02143 return ; }
02144
02145
02146
02147
02162 bool
02163 SbbModel::solveWithCuts (OsiCuts &cuts, int numberTries, SbbNode *node,
02164 int &numberOldActiveCuts, int &numberNewCuts,
02165 int &maximumWhich, int *&whichGenerator)
02166
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186 { bool feasible = true ;
02187 int lastNumberCuts = 0 ;
02188 double lastObjective = -1.0e100 ;
02189 int violated = 0 ;
02190 int numberRowsAtStart = solver_->getNumRows() ;
02191 int numberColumns = solver_->getNumCols() ;
02192
02193 numberOldActiveCuts = numberRowsAtStart-numberRowsAtContinuous_ ;
02194 numberNewCuts = 0 ;
02195
02196 # ifdef SBB_DEBUG
02197
02198
02199
02200
02201
02202 bool onOptimalPath = false ;
02203 const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
02204 if (debugger)
02205 onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
02206 # endif
02207
02208
02209
02210
02211
02212 feasible = resolve() ;
02213
02214 #ifdef SBB_DEBUG
02215 if (feasible)
02216 { printf("Obj value %g (%s) %d rows\n",solver_->getObjValue(),
02217 (solver_->isProvenOptimal())?"proven":"unproven",
02218 solver_->getNumRows()) ; }
02219 else
02220 { printf("Infeasible %d rows\n",solver_->getNumRows()) ; }
02221
02222
02223
02224
02225
02226 if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
02227 assert(feasible) ;
02228 # endif
02229
02230 if (!feasible) return (false) ;
02231
02232
02233
02234 reducedCostFix() ;
02235 const double *lower = solver_->getColLower() ;
02236 const double *upper = solver_->getColUpper() ;
02237 const double *solution = solver_->getColSolution() ;
02238
02239
02240
02241
02242
02243 double minimumDrop = minimumDrop_ ;
02244 if (numberTries<0)
02245 { numberTries = -numberTries ;
02246 minimumDrop = -1.0 ; }
02247
02248
02249
02250
02251 # define SCANCUTS 100
02252 int *countColumnCuts = NULL ;
02253 int *countRowCuts = NULL ;
02254 bool fullScan = false ;
02255 if ((numberNodes_%SCANCUTS) == 0)
02256 { fullScan = true ;
02257 countColumnCuts = new int[numberCutGenerators_+numberHeuristics_] ;
02258 countRowCuts = new int[numberCutGenerators_+numberHeuristics_] ;
02259 memset(countColumnCuts,0,
02260 (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ;
02261 memset(countRowCuts,0,
02262 (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ; }
02263
02264 double direction = solver_->getObjSense() ;
02265 double startObjective = solver_->getObjValue()*direction ;
02266
02267 int numberPasses = 0 ;
02268 double primalTolerance = 1.0e-7 ;
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282 do
02283 { numberPasses++ ;
02284 numberTries-- ;
02285 OsiCuts theseCuts ;
02286
02287
02288
02289
02290
02291
02292
02293
02294 if (numberPasses == 1 && howOftenGlobalScan_ > 0 &&
02295 (numberNodes_%howOftenGlobalScan_) == 0)
02296 { int numberCuts = globalCuts_.sizeColCuts() ;
02297 int i;
02298 for ( i = 0 ; i < numberCuts ; i++)
02299 { const OsiColCut *thisCut = globalCuts_.colCutPtr(i) ;
02300 if (thisCut->violated(solution)>primalTolerance) {
02301 printf("Global cut added - violation %g\n",
02302 thisCut->violated(solution)) ;
02303 theseCuts.insert(*thisCut) ;
02304 }
02305 }
02306 numberCuts = globalCuts_.sizeRowCuts() ;
02307 for ( i = 0;i<numberCuts;i++) {
02308 const OsiRowCut * thisCut = globalCuts_.rowCutPtr(i) ;
02309 if (thisCut->violated(solution)>primalTolerance) {
02310 printf("Global cut added - violation %g\n",
02311 thisCut->violated(solution)) ;
02312 theseCuts.insert(*thisCut) ;
02313 }
02314 }
02315 }
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326 double * newSolution = new double [numberColumns] ;
02327 double heuristicValue = getCutoff() ;
02328 int found = -1;
02329
02330 for (int i = 0;i<numberCutGenerators_+numberHeuristics_;i++) {
02331 int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
02332 int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
02333 if (i<numberCutGenerators_) {
02334 bool mustResolve =
02335 generator_[i]->generateCuts(theseCuts,fullScan) ;
02336 if (mustResolve) {
02337 feasible = resolve() ;
02338 # ifdef SBB_DEBUG
02339 debugger = solver_->getRowCutDebugger() ;
02340 if (debugger)
02341 onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
02342 if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
02343 assert(feasible) ;
02344 # endif
02345 if (!feasible)
02346 break ;
02347 }
02348 } else {
02349
02350 double saveValue = heuristicValue ;
02351 int ifSol =
02352 heuristic_[i-numberCutGenerators_]->solution(heuristicValue,
02353 newSolution,
02354 theseCuts) ;
02355 if (ifSol>0) {
02356
02357 found = i ;
02358 } else if (ifSol<0) {
02359 heuristicValue = saveValue ;
02360 }
02361 }
02362 int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
02363 int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
02364
02365 #ifdef SBB_DEBUG
02366 if (onOptimalPath) {
02367 int k ;
02368 for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
02369 OsiRowCut thisCut = theseCuts.rowCut(k) ;
02370 assert(!debugger->invalidCut(thisCut)) ;
02371 }
02372 }
02373 #endif
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386 int numberBefore =
02387 numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
02388 int numberAfter =
02389 numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
02390 if (numberAfter > maximumWhich) {
02391 maximumWhich = max(maximumWhich*2+100,numberAfter) ;
02392 int * temp = new int[2*maximumWhich] ;
02393 memcpy(temp,whichGenerator,numberBefore*sizeof(int)) ;
02394 delete [] whichGenerator ;
02395 whichGenerator = temp ;
02396 }
02397 int j ;
02398 if (fullScan) {
02399
02400 countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
02401 countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
02402 }
02403
02404 for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
02405 whichGenerator[numberBefore++] = i ;
02406 const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
02407 if (thisCut->globallyValid()) {
02408
02409 globalCuts_.insert(*thisCut) ;
02410 }
02411 }
02412 for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
02413 whichGenerator[numberBefore++] = i ;
02414 const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
02415 if (thisCut->globallyValid()) {
02416
02417 globalCuts_.insert(*thisCut) ;
02418 }
02419 }
02420 }
02421
02422
02423
02424
02425
02426
02427 if (found >= 0)
02428 { setBestSolution(SBB_ROUNDING,heuristicValue,newSolution) ; }
02429 delete [] newSolution ;
02430
02431 #if 0
02432
02433 theseCuts.printCuts() ;
02434 #endif
02435 int numberColumnCuts = theseCuts.sizeColCuts() ;
02436 int numberRowCuts = theseCuts.sizeRowCuts() ;
02437 violated = numberRowCuts + numberColumnCuts ;
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449
02450
02451
02452
02453 if (numberColumnCuts) {
02454
02455 #ifdef SBB_DEBUG
02456 double * oldLower = new double [numberColumns] ;
02457 double * oldUpper = new double [numberColumns] ;
02458 memcpy(oldLower,lower,numberColumns*sizeof(double)) ;
02459 memcpy(oldUpper,upper,numberColumns*sizeof(double)) ;
02460 #endif
02461
02462 double integerTolerance = getDblParam(SbbIntegerTolerance) ;
02463 for (int i = 0;i<numberColumnCuts;i++) {
02464 const OsiColCut * thisCut = theseCuts.colCutPtr(i) ;
02465 const CoinPackedVector & lbs = thisCut->lbs() ;
02466 const CoinPackedVector & ubs = thisCut->ubs() ;
02467 int j ;
02468 int n ;
02469 const int * which ;
02470 const double * values ;
02471 n = lbs.getNumElements() ;
02472 which = lbs.getIndices() ;
02473 values = lbs.getElements() ;
02474 for (j = 0;j<n;j++) {
02475 int iColumn = which[j] ;
02476 double value = solution[iColumn] ;
02477 #ifdef SBB_DEBUG
02478 printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn],
02479 solution[iColumn],oldUpper[iColumn],values[j]) ;
02480 #endif
02481 solver_->setColLower(iColumn,values[j]) ;
02482 if (value<values[j]-integerTolerance)
02483 violated = -1 ;
02484 if (values[j]>upper[iColumn]+integerTolerance) {
02485
02486 violated = -2 ;
02487 break ;
02488 }
02489 }
02490 n = ubs.getNumElements() ;
02491 which = ubs.getIndices() ;
02492 values = ubs.getElements() ;
02493 for (j = 0;j<n;j++) {
02494 int iColumn = which[j] ;
02495 double value = solution[iColumn] ;
02496 #ifdef SBB_DEBUG
02497 printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn],
02498 solution[iColumn],oldUpper[iColumn],values[j]) ;
02499 #endif
02500 solver_->setColUpper(iColumn,values[j]) ;
02501 if (value>values[j]+integerTolerance)
02502 violated = -1 ;
02503 if (values[j]<lower[iColumn]-integerTolerance) {
02504
02505 violated = -2 ;
02506 break ;
02507 }
02508 }
02509 }
02510 #ifdef SBB_DEBUG
02511 delete [] oldLower ;
02512 delete [] oldUpper ;
02513 #endif
02514 }
02515
02516
02517
02518
02519 if (violated == -2) {
02520
02521 feasible = false ;
02522 break ;
02523 }
02524
02525
02526
02527
02528
02529
02530
02531 int numberRowsNow = solver_->getNumRows() ;
02532 assert(numberRowsNow == numberRowsAtStart+lastNumberCuts) ;
02533 int numberToAdd = theseCuts.sizeRowCuts() ;
02534 numberNewCuts = lastNumberCuts+numberToAdd ;
02535
02536
02537
02538
02539 delete basis_ ;
02540 basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
02541 assert(basis_ != NULL);
02542 basis_->resize(numberRowsAtStart+numberNewCuts,numberColumns) ;
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562
02563
02564
02565
02566
02567
02568 if (numberRowCuts > 0 || numberColumnCuts > 0)
02569 { if (numberToAdd > 0)
02570 { int i ;
02571 OsiRowCut * addCuts = new OsiRowCut [numberToAdd] ;
02572 for (i = 0 ; i < numberToAdd ; i++)
02573 { addCuts[i] = theseCuts.rowCut(i) ; }
02574 solver_->applyRowCuts(numberToAdd,addCuts) ;
02575
02576 delete [] addCuts ;
02577 for (i = 0 ; i < numberToAdd ; i++)
02578 { cuts.insert(theseCuts.rowCut(i)) ; }
02579 for (i = 0 ; i < numberToAdd ; i++)
02580 { basis_->setArtifStatus(numberRowsNow+i,
02581 CoinWarmStartBasis::basic) ; }
02582 if (solver_->setWarmStart(basis_) == false)
02583 { throw CoinError("Fail setWarmStart() after cut installation.",
02584 "solveWithCuts","SbbModel") ; } }
02585 feasible = resolve() ;
02586 # ifdef SBB_DEBUG
02587 printf("Obj value after cuts %g %d rows\n",solver_->getObjValue(),
02588 solver_->getNumRows()) ;
02589 if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
02590 assert(feasible) ;
02591 # endif
02592 }
02593
02594
02595
02596 else
02597 { numberTries = 0 ; }
02598
02599
02600
02601
02602
02603
02604
02605
02606
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618
02619
02620
02621
02622 if (feasible)
02623 { int cutIterations = solver_->getIterationCount() ;
02624 takeOffCuts(cuts,whichGenerator,numberOldActiveCuts,numberNewCuts,true) ;
02625 if (solver_->isDualObjectiveLimitReached())
02626 { feasible = false ;
02627 # ifdef SBB_DEBUG
02628 double z = solver_->getObjValue() ;
02629 double cut = getCutoff() ;
02630 printf("Lost feasibility by %g in takeOffCuts; z = %g, cutoff = %g\n",
02631 z-cut,z,cut) ;
02632 # endif
02633 }
02634 if (feasible)
02635 { numberRowsAtStart = numberOldActiveCuts+numberRowsAtContinuous_ ;
02636 lastNumberCuts = numberNewCuts ;
02637 if (direction*solver_->getObjValue() < lastObjective+minimumDrop &&
02638 numberPasses >= 3)
02639 { numberTries = 0 ; }
02640 if (numberRowCuts+numberColumnCuts == 0 || cutIterations == 0)
02641 { break ; }
02642 if (numberTries > 0)
02643 { reducedCostFix() ;
02644 lastObjective = direction*solver_->getObjValue() ;
02645 lower = solver_->getColLower() ;
02646 upper = solver_->getColUpper() ;
02647 solution = solver_->getColSolution() ; } } }
02648
02649
02650
02651
02652
02653
02654
02655 if (!feasible)
02656 { int i ;
02657 for (i = 0;i<currentNumberCuts_;i++) {
02658
02659 if (addedCuts_[i]) {
02660 if (!addedCuts_[i]->decrement())
02661 delete addedCuts_[i] ;
02662 addedCuts_[i] = NULL ;
02663 }
02664 }
02665 numberTries = 0 ;
02666 }
02667 } while (numberTries) ;
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686 if (feasible&&fullScan&&numberCutGenerators_) {
02687
02688 int i ;
02689 double thisObjective = solver_->getObjValue()*direction ;
02690 double totalCuts = 0.0 ;
02691 for (i = 0;i<numberCutGenerators_;i++)
02692 totalCuts += countRowCuts[i] + 5.0*countColumnCuts[i] ;
02693 if (!numberNodes_)
02694 handler_->message(SBB_ROOT,messages_)
02695 <<numberNewCuts
02696 <<startObjective<<thisObjective
02697 <<numberPasses
02698 <<CoinMessageEol ;
02699 int * count = new int[numberCutGenerators_] ;
02700 memset(count,0,numberCutGenerators_*sizeof(int)) ;
02701 for (i = 0;i<numberNewCuts;i++)
02702 count[whichGenerator[i]]++ ;
02703 double small = (0.5* totalCuts) /
02704 ((double) numberCutGenerators_) ;
02705 for (i = 0;i<numberCutGenerators_;i++) {
02706 int howOften = generator_[i]->howOften() ;
02707 if (howOften<-99)
02708 continue ;
02709 if (howOften<0||howOften >= 1000000) {
02710
02711 double thisCuts = countRowCuts[i] + 5.0*countColumnCuts[i] ;
02712 if (!thisCuts||howOften == -99) {
02713 if (howOften == -99)
02714 howOften = -100 ;
02715 else
02716 howOften = 1000000+SCANCUTS;
02717 } else if (thisCuts<small) {
02718 int k = (int) sqrt(small/thisCuts) ;
02719 howOften = k+1000000 ;
02720 } else {
02721 howOften = 1+1000000 ;
02722 }
02723 }
02724 generator_[i]->setHowOften(howOften) ;
02725 int newFrequency = generator_[i]->howOften()%1000000 ;
02726 if (handler_->logLevel()>1||!numberNodes_)
02727 handler_->message(SBB_GENERATOR,messages_)
02728 <<i
02729 <<generator_[i]->cutGeneratorName()
02730 <<countRowCuts[i]<<count[i]
02731 <<countColumnCuts[i]
02732 <<newFrequency
02733 <<CoinMessageEol ;
02734 }
02735 delete [] count ;
02736 }
02737
02738 delete [] countRowCuts ;
02739 delete [] countColumnCuts ;
02740
02741
02742 #ifdef CHECK_CUT_COUNTS
02743 if (feasible)
02744 { delete basis_ ;
02745 basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
02746 printf("solveWithCuts: Number of rows at end (only active cuts) %d\n",
02747 numberRowsAtContinuous_+numberNewCuts+numberOldActiveCuts) ;
02748 basis_->print() ; }
02749 #endif
02750 #ifdef SBB_DEBUG
02751 if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
02752 assert(feasible) ;
02753 #endif
02754
02755 return feasible ; }
02756
02757
02758
02759
02760
02761
02762
02763
02764
02765
02766
02767
02768
02769
02770
02771
02772
02773
02774
02775
02776
02777
02778
02779
02780
02781
02782
02783
02784
02785
02786
02787
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797
02798
02799
02800
02801 void
02802 SbbModel::takeOffCuts (OsiCuts &newCuts, int *whichGenerator,
02803 int &numberOldActiveCuts, int &numberNewCuts,
02804 bool allowResolve)
02805
02806 {
02807 int firstOldCut = numberRowsAtContinuous_ ;
02808 int totalNumberCuts = numberNewCuts+numberOldActiveCuts ;
02809 int *solverCutIndices = new int[totalNumberCuts] ;
02810 int *newCutIndices = new int[numberNewCuts] ;
02811 const CoinWarmStartBasis* ws ;
02812 CoinWarmStartBasis::Status status ;
02813 bool needPurge = true ;
02814
02815
02816
02817
02818
02819 while (needPurge)
02820 { int numberNewToDelete = 0 ;
02821 int numberOldToDelete = 0 ;
02822 int i ;
02823 ws = dynamic_cast<const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
02824
02825
02826
02827
02828
02829
02830
02831 int oldCutIndex = 0 ;
02832 for (i = 0 ; i < numberOldActiveCuts ; i++)
02833 { status = ws->getArtifStatus(i+firstOldCut) ;
02834 while (!addedCuts_[oldCutIndex]) oldCutIndex++ ;
02835 assert(oldCutIndex < currentNumberCuts_) ;
02836 if (status == CoinWarmStartBasis::basic)
02837 { solverCutIndices[numberOldToDelete++] = i+firstOldCut ;
02838 if (addedCuts_[oldCutIndex]->decrement() == 0)
02839 delete addedCuts_[oldCutIndex] ;
02840 addedCuts_[oldCutIndex] = NULL ;
02841 oldCutIndex++ ; }
02842 else
02843 { oldCutIndex++ ; } }
02844
02845
02846
02847
02848
02849
02850
02851 int firstNewCut = firstOldCut+numberOldActiveCuts ;
02852 int k = 0 ;
02853 for (i = 0 ; i < numberNewCuts ; i++)
02854 { status = ws->getArtifStatus(i+firstNewCut) ;
02855 if (status == CoinWarmStartBasis::basic)
02856 { solverCutIndices[numberNewToDelete+numberOldToDelete] = i+firstNewCut ;
02857 newCutIndices[numberNewToDelete++] = i ; }
02858 else
02859 {
02860 whichGenerator[k++] = whichGenerator[i] ; } }
02861 for (i = numberNewToDelete-1 ; i >= 0 ; i--)
02862 { int iCut = newCutIndices[i] ;
02863 newCuts.eraseRowCut(iCut) ; }
02864
02865
02866
02867
02868
02869
02870 if (numberNewToDelete+numberOldToDelete > 0)
02871 { solver_->deleteRows(numberNewToDelete+numberOldToDelete,
02872 solverCutIndices) ;
02873 numberNewCuts -= numberNewToDelete ;
02874 numberOldActiveCuts -= numberOldToDelete ;
02875 # ifdef SBB_DEBUG
02876 std::cout << "takeOffCuts: purged " << numberOldToDelete << "+"
02877 << numberNewToDelete << " cuts." << std::endl ;
02878 # endif
02879 if (allowResolve)
02880 { solver_->resolve() ;
02881 if (solver_->getIterationCount() == 0)
02882 { needPurge = false ; }
02883 # ifdef SBB_DEBUG
02884 else
02885 { std::cout << "Repeating purging loop. "
02886 << solver_->getIterationCount() << " iters."
02887 << std::endl ; }
02888 # endif
02889 }
02890 else
02891 { needPurge = false ; } }
02892 else
02893 { needPurge = false ; } }
02894
02895
02896
02897 delete ws ;
02898 delete [] solverCutIndices ;
02899 delete [] newCutIndices ;
02900 }
02901
02902
02903
02904 bool
02905 SbbModel::resolve()
02906 {
02907
02908 int iRow;
02909 int numberRows = solver_->getNumRows();
02910 const double * rowLower = solver_->getRowLower();
02911 const double * rowUpper = solver_->getRowUpper();
02912 bool feasible=true;
02913 for (iRow= numberRowsAtContinuous_;iRow<numberRows;iRow++) {
02914 if (rowLower[iRow]>rowUpper[iRow]+1.0e-8)
02915 feasible=false;
02916 }
02917
02918
02919
02920
02921
02922 if (feasible)
02923 { solver_->resolve() ;
02924 numberIterations_ += getIterationCount() ;
02925 feasible = (solver_->isProvenOptimal() &&
02926 !solver_->isDualObjectiveLimitReached()) ; }
02927 return feasible ; }
02928
02929
02930
02931
02932
02933
02934
02935
02936 SbbModel *
02937 SbbModel::findCliques(bool makeEquality,
02938 int atLeastThisMany, int lessThanThis, int defaultValue)
02939 {
02940
02941 assert(numberObjects_==numberIntegers_);
02942 CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
02943 int numberRows = solver_->getNumRows();
02944 int numberColumns = solver_->getNumCols();
02945
02946
02947 int numberSlacks=0;
02948 int * rows = new int[numberRows];
02949 double * element =new double[numberRows];
02950
02951 int iRow;
02952
02953 findIntegers(true);
02954 numberObjects_=numberIntegers_;
02955
02956 int numberCliques=0;
02957 SbbObject ** object = new SbbObject * [numberRows];
02958 int * which = new int[numberIntegers_];
02959 char * type = new char[numberIntegers_];
02960 int * lookup = new int[numberColumns];
02961 int i;
02962 for (i=0;i<numberColumns;i++)
02963 lookup[i]=-1;
02964 for (i=0;i<numberIntegers_;i++)
02965 lookup[integerVariable_[i]]=i;
02966
02967
02968 int totalP1=0,totalM1=0;
02969 int numberBig=0,totalBig=0;
02970 int numberFixed=0;
02971
02972
02973 const double * elementByRow = matrixByRow.getElements();
02974 const int * column = matrixByRow.getIndices();
02975 const int * rowStart = matrixByRow.getVectorStarts();
02976 const int * rowLength = matrixByRow.getVectorLengths();
02977
02978
02979 const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
02980
02981 const double * lower = getColLower();
02982 const double * upper = getColUpper();
02983 const double * rowLower = getRowLower();
02984 const double * rowUpper = getRowUpper();
02985
02986 for (iRow=0;iRow<numberRows;iRow++) {
02987 int numberP1=0, numberM1=0;
02988 int j;
02989 double upperValue=rowUpper[iRow];
02990 double lowerValue=rowLower[iRow];
02991 bool good=true;
02992 int slack = -1;
02993 for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
02994 int iColumn = column[j];
02995 int iInteger=lookup[iColumn];
02996 if (upper[iColumn]-lower[iColumn]<1.0e-8) {
02997
02998 upperValue -= lower[iColumn]*elementByRow[j];
02999 lowerValue -= lower[iColumn]*elementByRow[j];
03000 continue;
03001 } else if (upper[iColumn]!=1.0||lower[iColumn]!=0.0) {
03002 good = false;
03003 break;
03004 } else if (iInteger<0) {
03005 good = false;
03006 break;
03007 } else {
03008 if (columnLength[iColumn]==1)
03009 slack = iInteger;
03010 }
03011 if (fabs(elementByRow[j])!=1.0) {
03012 good=false;
03013 break;
03014 } else if (elementByRow[j]>0.0) {
03015 which[numberP1++]=iInteger;
03016 } else {
03017 numberM1++;
03018 which[numberIntegers_-numberM1]=iInteger;
03019 }
03020 }
03021 int iUpper = (int) floor(upperValue+1.0e-5);
03022 int iLower = (int) ceil(lowerValue-1.0e-5);
03023 int state=0;
03024 if (upperValue<1.0e6) {
03025 if (iUpper==1-numberM1)
03026 state=1;
03027 else if (iUpper==-numberM1)
03028 state=2;
03029 else if (iUpper<-numberM1)
03030 state=3;
03031 }
03032 if (!state&&lowerValue>-1.0e6) {
03033 if (-iLower==1-numberP1)
03034 state=-1;
03035 else if (-iLower==-numberP1)
03036 state=-2;
03037 else if (-iLower<-numberP1)
03038 state=-3;
03039 }
03040 if (good&&state) {
03041 if (abs(state)==3) {
03042
03043 numberObjects_=-1;
03044 break;
03045 } else if (abs(state)==2) {
03046
03047 numberFixed += numberP1+numberM1;
03048 if (state>0) {
03049
03050 for (i=0;i<numberP1;i++)
03051 solver_->setColUpper(integerVariable_[which[i]],0.0);
03052 for (i=0;i<numberM1;i++)
03053 solver_->setColLower(integerVariable_[which[numberIntegers_-i-1]],
03054 1.0);
03055 } else {
03056
03057 for (i=0;i<numberP1;i++)
03058 solver_->setColLower(integerVariable_[which[i]],1.0);
03059 for (i=0;i<numberM1;i++)
03060 solver_->setColUpper(integerVariable_[which[numberIntegers_-i-1]],
03061 0.0);
03062 }
03063 } else {
03064 int length = numberP1+numberM1;
03065 if (length >= atLeastThisMany&&length<lessThanThis) {
03066
03067 bool addOne=false;
03068 int objectType;
03069 if (iLower==iUpper) {
03070 objectType=1;
03071 } else {
03072 if (makeEquality) {
03073 objectType=1;
03074 element[numberSlacks]=state;
03075 rows[numberSlacks++]=iRow;
03076 addOne=true;
03077 } else {
03078 objectType=0;
03079 }
03080 }
03081 if (state>0) {
03082 totalP1 += numberP1;
03083 totalM1 += numberM1;
03084 for (i=0;i<numberP1;i++)
03085 type[i]=1;
03086 for (i=0;i<numberM1;i++) {
03087 which[numberP1]=which[numberIntegers_-i-1];
03088 type[numberP1++]=0;
03089 }
03090 } else {
03091 totalP1 += numberM1;
03092 totalM1 += numberP1;
03093 for (i=0;i<numberP1;i++)
03094 type[i]=0;
03095 for (i=0;i<numberM1;i++) {
03096 which[numberP1]=which[numberIntegers_-i-1];
03097 type[numberP1++]=1;
03098 }
03099 }
03100 if (addOne) {
03101
03102 which[numberP1]=numberIntegers_+numberSlacks-1;
03103 slack = numberP1;
03104 type[numberP1++]=1;
03105 } else if (slack >= 0) {
03106 for (i=0;i<numberP1;i++) {
03107 if (which[i]==slack) {
03108 slack=i;
03109 }
03110 }
03111 }
03112 object[numberCliques++] = new SbbClique(this,objectType,numberP1,
03113 which,type,
03114 1000000+numberCliques,slack);
03115 } else if (numberP1+numberM1 >= lessThanThis) {
03116
03117 numberBig++;
03118 totalBig += numberP1+numberM1;
03119 }
03120 }
03121 }
03122 }
03123 delete [] which;
03124 delete [] type;
03125 delete [] lookup;
03126 if (numberCliques<0) {
03127 printf("*** Problem infeasible\n");
03128 } else {
03129 if (numberCliques)
03130 printf("%d cliques of average size %g found, %d P1, %d M1\n",
03131 numberCliques,
03132 ((double)(totalP1+totalM1))/((double) numberCliques),
03133 totalP1,totalM1);
03134 else
03135 printf("No cliques found\n");
03136 if (numberBig)
03137 printf("%d large cliques ( >= %d) found, total %d\n",
03138 numberBig,lessThanThis,totalBig);
03139 if (numberFixed)
03140 printf("%d variables fixed\n",numberFixed);
03141 }
03142 if (numberCliques>0&&numberSlacks&&makeEquality) {
03143 printf("adding %d integer slacks\n",numberSlacks);
03144
03145 int * temp = new int[numberIntegers_+numberSlacks];
03146 memcpy(temp,integerVariable_,numberIntegers_*sizeof(int));
03147
03148 SbbModel * newModel = new SbbModel(*this);
03149 OsiSolverInterface * newSolver = newModel->solver();
03150 for (i=0;i<numberSlacks;i++) {
03151 temp[i+numberIntegers_]=i+numberColumns;
03152 int iRow = rows[i];
03153 double value = element[i];
03154 double lowerValue = 0.0;
03155 double upperValue = 1.0;
03156 double objValue = 0.0;
03157 CoinPackedVector column(1,&iRow,&value);
03158 newSolver->addCol(column,lowerValue,upperValue,objValue);
03159
03160 newSolver->setInteger(numberColumns+i);
03161 if (value >0)
03162 newSolver->setRowLower(iRow,rowUpper[iRow]);
03163 else
03164 newSolver->setRowUpper(iRow,rowLower[iRow]);
03165 }
03166
03167 for (i=0;i<newModel->numberObjects_;i++)
03168 delete newModel->object_[i];
03169 newModel->numberObjects_ = 0;
03170 delete [] newModel->object_;
03171 newModel->object_=NULL;
03172 newModel->findIntegers(true);
03173 if (priority_) {
03174
03175 delete [] newModel->priority_;
03176 newModel->priority_ = new int[newModel->numberIntegers_+numberCliques];
03177 memcpy(newModel->priority_,priority_,numberIntegers_*sizeof(int));
03178 for (i=numberIntegers_;i<newModel->numberIntegers_+numberCliques;i++)
03179 newModel->priority_[i]=defaultValue;
03180 }
03181 if (originalColumns_) {
03182
03183 delete [] newModel->originalColumns_;
03184 newModel->originalColumns_ = new int[numberColumns+numberSlacks];
03185 memcpy(newModel->originalColumns_,originalColumns_,numberColumns*sizeof(int));
03186
03187 for (i=numberColumns;i<numberColumns+numberSlacks;i++)
03188 newModel->originalColumns_[i]=-1;
03189 }
03190 delete [] rows;
03191 delete [] element;
03192 newModel->addObjects(numberCliques,object);
03193 for (;i<numberCliques;i++)
03194 delete object[i];
03195 delete [] object;
03196 newModel->synchronizeModel();
03197 return newModel;
03198 } else {
03199 if (numberCliques>0) {
03200 addObjects(numberCliques,object);
03201 for (;i<numberCliques;i++)
03202 delete object[i];
03203 synchronizeModel();
03204 }
03205 delete [] object;
03206 if (priority_) {
03207
03208 int * temp = new int[numberIntegers_+numberCliques];
03209 memcpy(temp,priority_,numberIntegers_*sizeof(int));
03210 delete [] priority_;
03211 priority_=temp;
03212 for (i=numberIntegers_;i<numberIntegers_+numberCliques;i++)
03213 priority_[i]=defaultValue;
03214 }
03215 delete [] rows;
03216 delete [] element;
03217 return this;
03218 }
03219 }
03220
03221
03222
03223
03224
03225
03226
03227
03228
03229
03230 void
03231 SbbModel::passInPriorities (const int * priorities,
03232 bool ifObject, int defaultValue)
03233 {
03234 findIntegers(false);
03235 int i;
03236 if (!priority_) {
03237 priority_ = new int[numberObjects_];
03238 for (i=0;i<numberObjects_;i++)
03239 priority_[i]=defaultValue;
03240 }
03241 if (priorities) {
03242 if (ifObject)
03243 memcpy(priority_+numberIntegers_,priorities,
03244 (numberObjects_-numberIntegers_)*sizeof(int));
03245 else
03246 memcpy(priority_,priorities,numberIntegers_*sizeof(int));
03247 }
03248 }
03249
03250
03251 void
03252 SbbModel::deleteObjects()
03253 {
03254 delete [] priority_;
03255 priority_=NULL;
03256 int i;
03257 for (i=0;i<numberObjects_;i++)
03258 delete object_[i];
03259 delete [] object_;
03260 object_ = NULL;
03261 numberObjects_=0;
03262 findIntegers(true);
03263 }
03264
03269 void SbbModel::synchronizeModel()
03270 {
03271 int i;
03272 for (i=0;i<numberHeuristics_;i++)
03273 heuristic_[i]->setModel(this);
03274 for (i=0;i<numberObjects_;i++)
03275 object_[i]->setModel(this);
03276 for (i=0;i<numberCutGenerators_;i++)
03277 generator_[i]->refreshModel(this);
03278 }
03279
03280
03281
03299 void
03300 SbbModel::findIntegers(bool startAgain)
03301 {
03302 assert(solver_);
03303
03304
03305
03306
03307 if (numberIntegers_&&!startAgain&&object_)
03308 return;
03309
03310
03311
03312
03313 delete [] integerVariable_;
03314 numberIntegers_=0;
03315 int numberColumns = getNumCols();
03316 int iColumn;
03317 for (iColumn=0;iColumn<numberColumns;iColumn++) {
03318 if (isInteger(iColumn))
03319 numberIntegers_++;
03320 }
03321
03322
03323
03324
03325 if (numberIntegers_) {
03326 if (!object_) {
03327 object_ = new SbbObject * [numberIntegers_];
03328 memset(object_,0,numberIntegers_*sizeof(SbbObject *));
03329 numberObjects_=numberIntegers_;
03330 }
03331 integerVariable_ = new int [numberIntegers_];
03332
03333
03334
03335
03336
03337 numberIntegers_=0;
03338 for (iColumn=0;iColumn<numberColumns;iColumn++) {
03339 if(isInteger(iColumn)) {
03340 delete object_[numberIntegers_];
03341 object_[numberIntegers_] =
03342 new SbbSimpleInteger(this,numberIntegers_,iColumn);
03343 integerVariable_[numberIntegers_++]=iColumn;
03344 }
03345 }
03346 } else {
03347 handler_->message(SBB_NOINT,messages_) << CoinMessageEol ;
03348 return;
03349 }
03350 }
03351
03352
03353
03354 void
03355 SbbModel::addObjects(int numberObjects, SbbObject ** objects)
03356 {
03357 int newNumberObjects= numberObjects+numberObjects_;
03358 SbbObject ** temp = new SbbObject * [newNumberObjects];
03359 int i;
03360 for (i=0;i<numberObjects_;i++)
03361 temp[i]=object_[i];
03362 for (;i<newNumberObjects;i++) {
03363 temp[i]=objects[i-numberObjects_]->clone();
03364 temp[i]->setModel(this);
03365 }
03366 delete [] object_;
03367 object_ = temp;
03368 numberObjects_ = newNumberObjects;
03369 }
03370
03388 void SbbModel::setCutoff (double value)
03389
03390 { double tol = 0 ;
03391 int fathomStrict = getIntParam(SbbFathomDiscipline) ;
03392 double direction = solver_->getObjSense();
03393 if (fathomStrict == 1)
03394 { solver_->getDblParam(OsiDualTolerance,tol) ;
03395 tol = tol*(1+fabs(value)) ;
03396
03397 value += tol ; }
03398
03399 solver_->setDblParam(OsiDualObjectiveLimit,value*direction); }
03400
03401
03402
03403
03404
03405
03406
03407
03408
03409
03410
03411
03412
03413
03414
03415
03416
03417 double
03418 SbbModel::checkSolution (double cutoff, const double *solution,
03419 bool fixVariables)
03420
03421 { int numberColumns = solver_->getNumCols();
03422
03423
03424
03425
03426
03427
03428
03429
03430
03431
03432 OsiSolverInterface * saveSolver = solver_;
03433 solver_ = continuousSolver_;
03434
03435 solver_->setColSolution(solution);
03436
03437 memcpy(currentSolution_,solver_->getColSolution(),
03438 numberColumns*sizeof(double));
03439
03440
03441
03442 double * saveUpper = new double[numberColumns];
03443 double * saveLower = new double[numberColumns];
03444 memcpy(saveUpper,getColUpper(),numberColumns*sizeof(double));
03445 memcpy(saveLower,getColLower(),numberColumns*sizeof(double));
03446
03447
03448
03449
03450
03451
03452 int i;
03453 for (i=0;i<numberObjects_;i++)
03454 object_[i]->feasibleRegion();
03455
03456
03457
03458
03459
03460
03461
03462
03463
03464
03465
03466
03467
03468
03469
03470
03471
03472
03473 CoinWarmStartBasis slack;
03474 solver_->setWarmStart(&slack);
03475
03476
03477
03478
03479
03480
03481 solver_->initialSolve();
03482
03483 if (!solver_->isProvenOptimal())
03484 { std::cout << "checkSolution infeas! Retrying wihout scaling."
03485 << std::endl ;
03486 bool saveTakeHint;
03487 OsiHintStrength saveStrength;
03488 bool savePrintHint;
03489 bool gotHint = (solver_->getHintParam(OsiDoReducePrint,savePrintHint,saveStrength));
03490 gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
03491 solver_->setHintParam(OsiDoScale,false,OsiHintTry);
03492 solver_->setHintParam(OsiDoReducePrint,false,OsiHintTry) ;
03493 solver_->initialSolve();
03494 solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
03495 solver_->setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ;
03496 }
03497 assert(solver_->isProvenOptimal());
03498 double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
03499
03500
03501
03502
03503
03504
03505
03506
03507
03508
03509
03510
03511 if (solver_->isProvenOptimal() && objectiveValue <= cutoff)
03512 { solution = solver_->getColSolution() ;
03513 memcpy(currentSolution_,solution,numberColumns*sizeof(double)) ;
03514
03515 const double * rowLower = solver_->getRowLower() ;
03516 const double * rowUpper = solver_->getRowUpper() ;
03517 int numberRows = solver_->getNumRows() ;
03518 double *rowActivity = new double[numberRows] ;
03519 memset(rowActivity,0,numberRows*sizeof(double)) ;
03520
03521 double integerTolerance = getIntegerTolerance() ;
03522
03523 int iColumn;
03524 for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
03525 { double value = currentSolution_[iColumn] ;
03526 value = max(value, saveLower[iColumn]) ;
03527 value = min(value, saveUpper[iColumn]) ;
03528 if (solver_->isInteger(iColumn))
03529 assert(fabs(value-currentSolution_[iColumn]) <= integerTolerance) ;
03530 currentSolution_[iColumn] = value ; }
03531
03532 solver_->getMatrixByCol()->times(currentSolution_,rowActivity) ;
03533 double primalTolerance ;
03534 solver_->getDblParam(OsiPrimalTolerance,primalTolerance) ;
03535 double largestInfeasibility =0.0;
03536 for (i=0 ; i < numberRows ; i++) {
03537 largestInfeasibility = max(largestInfeasibility,
03538 rowLower[i]-rowActivity[i]);
03539 largestInfeasibility = max(largestInfeasibility,
03540 rowActivity[i]-rowUpper[i]);
03541 }
03542 if (largestInfeasibility>100.0*primalTolerance)
03543 handler_->message(SBB_NOTFEAS3, messages_)
03544 << largestInfeasibility << CoinMessageEol ;
03545
03546 delete [] rowActivity ; }
03547 else
03548 { objectiveValue=1.0e50 ; }
03549
03550
03551
03552
03553
03554 if (!fixVariables)
03555 { for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++)
03556 { solver_->setColLower(iColumn,saveLower[iColumn]) ;
03557 solver_->setColUpper(iColumn,saveUpper[iColumn]) ; } }
03558 delete [] saveLower;
03559 delete [] saveUpper;
03560
03561
03562
03563
03564 solver_=saveSolver;
03565 return objectiveValue;
03566 }
03567
03568
03569
03570
03571
03572
03573
03574
03575
03576
03577
03578
03579
03580
03581
03582
03583
03584
03585
03586 void
03587 SbbModel::setBestSolution (SBB_Message how,
03588 double & objectiveValue, const double *solution,
03589 bool fixVariables)
03590
03591 { double cutoff = getCutoff() ;
03592
03593
03594
03595
03596 objectiveValue = checkSolution(cutoff,solution,fixVariables);
03597 if (objectiveValue > cutoff)
03598 { if (objectiveValue>1.0e30)
03599 handler_->message(SBB_NOTFEAS1, messages_) << CoinMessageEol ;
03600 else
03601 handler_->message(SBB_NOTFEAS2, messages_)
03602 << objectiveValue << cutoff << CoinMessageEol ; }
03603
03604
03605
03606
03607
03608 else
03609 { bestObjective_ = objectiveValue;
03610 int numberColumns = solver_->getNumCols();
03611 if (!bestSolution_)
03612 bestSolution_ = new double[numberColumns];
03613 memcpy(bestSolution_,currentSolution_,numberColumns*sizeof(double));
03614
03615 cutoff = bestObjective_-dblParam_[SbbCutoffIncrement];
03616
03617
03618
03619 setCutoff(cutoff);
03620
03621 if (how==SBB_ROUNDING)
03622 numberHeuristicSolutions_++;
03623 numberSolutions_++;
03624
03625 handler_->message(how,messages_)
03626 <<bestObjective_<<numberIterations_
03627 <<numberNodes_
03628 <<CoinMessageEol;
03629
03630
03631
03632
03633
03634
03635 OsiCuts theseCuts;
03636 int i;
03637 for (i=0;i<numberCutGenerators_;i++) {
03638 if (generator_[i]->atSolution())
03639 generator_[i]->generateCuts(theseCuts,true);
03640 }
03641 int numberCuts = theseCuts.sizeColCuts();
03642 for (i=0;i<numberCuts;i++) {
03643 const OsiColCut * thisCut = theseCuts.colCutPtr(i);
03644 if (thisCut->globallyValid()) {
03645
03646 globalCuts_.insert(*thisCut);
03647 }
03648 }
03649 numberCuts = theseCuts.sizeRowCuts();
03650 for (i=0;i<numberCuts;i++) {
03651 const OsiRowCut * thisCut = globalCuts_.rowCutPtr(i);
03652 if (thisCut->globallyValid()) {
03653
03654 globalCuts_.insert(*thisCut);
03655 }
03656 }
03657 }
03658
03659 return ; }
03660
03661
03662
03663
03664
03665
03666
03667
03668
03669
03670 bool
03671 SbbModel::feasibleSolution(int & numberIntegerInfeasibilities,
03672 int & numberObjectInfeasibilities) const
03673 {
03674 int numberUnsatisfied=0;
03675 double sumUnsatisfied=0.0;
03676 int preferredWay;
03677 double otherWay;
03678 double integerTolerance = getIntegerTolerance();
03679 int j;
03680
03681 memcpy(currentSolution_,solver_->getColSolution(),
03682 solver_->getNumCols()*sizeof(double));
03683 for (j=0;j<numberIntegers_;j++) {
03684 const SbbObject * object = object_[j];
03685 double infeasibility = object->infeasibility(preferredWay,otherWay);
03686 if (infeasibility>integerTolerance) {
03687 numberUnsatisfied++;
03688 sumUnsatisfied += infeasibility;
03689 }
03690 }
03691 numberIntegerInfeasibilities = numberUnsatisfied;
03692 for (;j<numberObjects_;j++) {
03693 const SbbObject * object = object_[j];
03694 double infeasibility = object->infeasibility(preferredWay,otherWay);
03695 if (infeasibility>integerTolerance) {
03696 numberUnsatisfied++;
03697 sumUnsatisfied += infeasibility;
03698 }
03699 }
03700 numberObjectInfeasibilities = numberUnsatisfied-numberIntegerInfeasibilities;
03701 return (!numberUnsatisfied);
03702 }
03703
03704
03705
03706
03707
03708
03709
03710 bool
03711 SbbModel::tightenVubs(int type, bool allowMultipleBinary, double useCutoff)
03712 {
03713
03714 CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
03715 int numberRows = solver_->getNumRows();
03716 int numberColumns = solver_->getNumCols();
03717
03718 int iRow,iColumn;
03719
03720
03721
03722 const int * column = matrixByRow.getIndices();
03723 const int * rowStart = matrixByRow.getVectorStarts();
03724 const int * rowLength = matrixByRow.getVectorLengths();
03725
03726 const double * colUpper = solver_->getColUpper();
03727 const double * colLower = solver_->getColLower();
03728
03729
03730
03731 const double * objective = solver_->getObjCoefficients();
03732
03733 const double * colsol = solver_->getColSolution();
03734
03735 int numberVub=0;
03736 int * continuous = new int[numberColumns];
03737 if (type >= 0) {
03738 double * sort = new double[numberColumns];
03739 for (iRow=0;iRow<numberRows;iRow++) {
03740 int j;
03741 int numberBinary=0;
03742 int numberUnsatisfiedBinary=0;
03743 int numberContinuous=0;
03744 int iCont=-1;
03745 double weight=1.0e30;
03746 for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
03747 int iColumn = column[j];
03748 if (colUpper[iColumn]-colLower[iColumn]>1.0e-8) {
03749 if (solver_->isFreeBinary(iColumn)) {
03750 numberBinary++;
03751
03752
03753
03754
03755 if (colsol[iColumn]>colLower[iColumn]+1.0e-6&&
03756 colsol[iColumn]<colUpper[iColumn]-1.0e-6) {
03757 numberUnsatisfiedBinary++;
03758 weight = min(weight,fabs(objective[iColumn]));
03759 }
03760 } else {
03761 numberContinuous++;
03762 iCont=iColumn;
03763 }
03764 }
03765 }
03766 if (numberContinuous==1&&numberBinary) {
03767 if (numberBinary==1||allowMultipleBinary) {
03768
03769 if (!numberUnsatisfiedBinary)
03770 weight=-1.0;
03771 sort[numberVub]=-weight;
03772 continuous[numberVub++] = iCont;
03773 }
03774 }
03775 }
03776 if (type>0) {
03777
03778 CoinSort_2(sort,sort+numberVub,continuous);
03779 numberVub = min(numberVub,type);
03780 }
03781 delete [] sort;
03782 } else {
03783 for (iColumn=0;iColumn<numberColumns;iColumn++)
03784 continuous[iColumn]=iColumn;
03785 numberVub=numberColumns;
03786 }
03787 bool feasible = tightenVubs(numberVub,continuous,useCutoff);
03788 delete [] continuous;
03789
03790 return feasible;
03791 }
03792
03793 bool
03794 SbbModel::tightenVubs(int numberSolves, const int * which,
03795 double useCutoff)
03796 {
03797
03798 int numberColumns = solver_->getNumCols();
03799
03800 int iColumn;
03801
03802 OsiSolverInterface * solver = solver_;
03803 double saveCutoff = getCutoff() ;
03804
03805 double * objective = new double[numberColumns];
03806 memcpy(objective,solver_->getObjCoefficients(),numberColumns*sizeof(double));
03807 double direction = solver_->getObjSense();
03808
03809
03810 if (useCutoff<1.0e30) {
03811
03812 solver = solver_->clone();
03813 CoinPackedVector newRow;
03814 for (iColumn=0;iColumn<numberColumns;iColumn++) {
03815 solver->setObjCoeff(iColumn,0.0);
03816 if (objective[iColumn])
03817 newRow.insert(iColumn,direction * objective[iColumn]);
03818
03819 }
03820 solver->addRow(newRow,-DBL_MAX,useCutoff);
03821
03822 delete [] objective;
03823 objective=NULL;
03824 }
03825 setCutoff(DBL_MAX);
03826
03827
03828 bool * vub = new bool [numberColumns];
03829 int iVub;
03830
03831
03832 for (iColumn=0;iColumn<numberColumns;iColumn++)
03833 vub[iColumn]=false;
03834 for (iVub=0;iVub<numberSolves;iVub++)
03835 vub[which[iVub]]=true;
03836 OsiCuts cuts;
03837
03838 CglProbing* generator = NULL;
03839 int iGen;
03840 for (iGen=0;iGen<numberCutGenerators_;iGen++) {
03841 generator = dynamic_cast<CglProbing*>(generator_[iGen]->generator());
03842 if (generator)
03843 break;
03844 }
03845 int numberFixed=0;
03846 int numberTightened=0;
03847 int numberFixedByProbing=0;
03848 int numberTightenedByProbing=0;
03849 int printFrequency = (numberSolves+19)/20;
03850 int save[4];
03851 if (generator) {
03852
03853 save[0]=generator->getMaxPass();
03854 save[1]=generator->getMaxProbe();
03855 save[2]=generator->getMaxLook();
03856 save[3]=generator->rowCuts();
03857 generator->setMaxPass(1);
03858 generator->setMaxProbe(10);
03859 generator->setMaxLook(50);
03860 generator->setRowCuts(0);
03861
03862
03863 generator->generateCutsAndModify(*solver,cuts);
03864 const double * tightLower = generator->tightLower();
03865 const double * lower = solver->getColLower();
03866 const double * tightUpper = generator->tightUpper();
03867 const double * upper = solver->getColUpper();
03868 for (iColumn=0;iColumn<numberColumns;iColumn++) {
03869 double newUpper = tightUpper[iColumn];
03870 double newLower = tightLower[iColumn];
03871 if (newUpper<upper[iColumn]-1.0e-8*(fabs(upper[iColumn])+1)||
03872 newLower>lower[iColumn]+1.0e-8*(fabs(lower[iColumn])+1)) {
03873 if (newUpper<newLower) {
03874 fprintf(stderr,"Problem is infeasible\n");
03875 return false;
03876 }
03877 if (newUpper==newLower) {
03878 numberFixed++;
03879 numberFixedByProbing++;
03880 solver->setColLower(iColumn,newLower);
03881 solver->setColUpper(iColumn,newUpper);
03882 } else if (vub[iColumn]) {
03883 numberTightened++;
03884 numberTightenedByProbing++;
03885 if (!solver->isInteger(iColumn)) {
03886
03887 newLower=max(lower[iColumn],
03888 newLower
03889 -1.0e-5*(fabs(lower[iColumn])+1));
03890 newUpper=min(upper[iColumn],
03891 newUpper
03892 +1.0e-5*(fabs(upper[iColumn])+1));
03893 }
03894 solver->setColLower(iColumn,newLower);
03895 solver->setColUpper(iColumn,newUpper);
03896 }
03897 }
03898 }
03899 }
03900 CoinWarmStart * ws = solver->getWarmStart();
03901 double * solution = new double [numberColumns];
03902 memcpy(solution,solver->getColSolution(),numberColumns*sizeof(double));
03903 for (iColumn=0;iColumn<numberColumns;iColumn++)
03904 solver->setObjCoeff(iColumn,0.0);
03905
03906 for (iVub=0;iVub<numberSolves;iVub++) {
03907 iColumn = which[iVub];
03908 int iTry;
03909 for (iTry=0;iTry<2;iTry++) {
03910 double saveUpper = solver->getColUpper()[iColumn];
03911 double saveLower = solver->getColLower()[iColumn];
03912 double value;
03913 if (iTry==1) {
03914
03915 solver->setObjCoeff(iColumn,-1.0);
03916 } else {
03917
03918 solver->setObjCoeff(iColumn,1.0);
03919 }
03920 solver->initialSolve();
03921 value = solver->getColSolution()[iColumn];
03922 bool change=false;
03923 if (iTry==1) {
03924 if (value<saveUpper-1.0e-4) {
03925 if (solver->isInteger(iColumn)) {
03926 value = floor(value+0.00001);
03927 } else {
03928
03929 value=min(saveUpper,value+1.0e-5*(fabs(saveUpper)+1));
03930 }
03931 if (value-saveLower<1.0e-7)
03932 value = saveLower;
03933 solver->setColUpper(iColumn,value);
03934 saveUpper=value;
03935 change=true;
03936 }
03937 } else {
03938 if (value>saveLower+1.0e-4) {
03939 if (solver->isInteger(iColumn)) {
03940 value = ceil(value-0.00001);
03941 } else {
03942
03943 value=max(saveLower,value-1.0e-5*(fabs(saveLower)+1));
03944 }
03945 if (saveUpper-value<1.0e-7)
03946 value = saveUpper;
03947 solver->setColLower(iColumn,value);
03948 saveLower=value;
03949 change=true;
03950 }
03951 }
03952 solver->setObjCoeff(iColumn,0.0);
03953 if (change) {
03954 if (saveUpper==saveLower)
03955 numberFixed++;
03956 else
03957 numberTightened++;
03958 int saveFixed=numberFixed;
03959
03960 int jColumn;
03961 if (generator) {
03962
03963 cuts = OsiCuts();
03964 generator->generateCutsAndModify(*solver,cuts);
03965 const double * tightLower = generator->tightLower();
03966 const double * lower = solver->getColLower();
03967 const double * tightUpper = generator->tightUpper();
03968 const double * upper = solver->getColUpper();
03969 for (jColumn=0;jColumn<numberColumns;jColumn++) {
03970 double newUpper = tightUpper[jColumn];
03971 double newLower = tightLower[jColumn];
03972 if (newUpper<upper[jColumn]-1.0e-8*(fabs(upper[jColumn])+1)||
03973 newLower>lower[jColumn]+1.0e-8*(fabs(lower[jColumn])+1)) {
03974 if (newUpper<newLower) {
03975 fprintf(stderr,"Problem is infeasible\n");
03976 return false;
03977 }
03978 if (newUpper==newLower) {
03979 numberFixed++;
03980 numberFixedByProbing++;
03981 solver->setColLower(jColumn,newLower);
03982 solver->setColUpper(jColumn,newUpper);
03983 } else if (vub[jColumn]) {
03984 numberTightened++;
03985 numberTightenedByProbing++;
03986 if (!solver->isInteger(jColumn)) {
03987
03988 newLower=max(lower[jColumn],
03989 newLower
03990 -1.0e-5*(fabs(lower[jColumn])+1));
03991 newUpper=min(upper[jColumn],
03992 newUpper
03993 +1.0e-5*(fabs(upper[jColumn])+1));
03994 }
03995 solver->setColLower(jColumn,newLower);
03996 solver->setColUpper(jColumn,newUpper);
03997 }
03998 }
03999 }
04000 }
04001 if (numberFixed>saveFixed) {
04002
04003
04004 if (objective) {
04005 for (jColumn=0;jColumn<numberColumns;jColumn++)
04006 solver->setObjCoeff(jColumn,objective[jColumn]);
04007 }
04008 solver->setColSolution(solution);
04009 solver->setWarmStart(ws);
04010 solver->resolve();
04011 if (!solver->isProvenOptimal()) {
04012 fprintf(stderr,"Problem is infeasible\n");
04013 return false;
04014 }
04015 delete ws;
04016 ws = solver->getWarmStart();
04017 memcpy(solution,solver->getColSolution(),
04018 numberColumns*sizeof(double));
04019 for (jColumn=0;jColumn<numberColumns;jColumn++)
04020 solver->setObjCoeff(jColumn,0.0);
04021 }
04022 }
04023 solver->setColSolution(solution);
04024 solver->setWarmStart(ws);
04025 }
04026 if (iVub%printFrequency==0)
04027 handler_->message(SBB_VUB_PASS,messages_)
04028 <<iVub+1<<numberFixed<<numberTightened
04029 <<CoinMessageEol;
04030 }
04031 handler_->message(SBB_VUB_END,messages_)
04032 <<numberFixed<<numberTightened
04033 <<CoinMessageEol;
04034 delete ws;
04035 delete [] solution;
04036
04037 if (objective) {
04038 for (iColumn=0;iColumn<numberColumns;iColumn++)
04039 solver_->setObjCoeff(iColumn,objective[iColumn]);
04040 delete [] objective;
04041 }
04042 delete [] vub;
04043 if (generator) {
04044
04045
04046
04047 if (generator_[iGen]->howOften()==-1&&
04048 (numberFixedByProbing+numberTightenedByProbing)*5>
04049 (numberFixed+numberTightened))
04050 generator_[iGen]->setHowOften(1000000+1);
04051 generator->setMaxPass(save[0]);
04052 generator->setMaxProbe(save[1]);
04053 generator->setMaxLook(save[2]);
04054 generator->setRowCuts(save[3]);
04055 }
04056
04057 if (solver!=solver_) {
04058
04059 const double * lower = solver->getColLower();
04060 const double * upper = solver->getColUpper();
04061 const double * lowerOrig = solver_->getColLower();
04062 const double * upperOrig = solver_->getColUpper();
04063 for (iColumn=0;iColumn<numberColumns;iColumn++) {
04064 solver_->setColLower(iColumn,max(lower[iColumn],lowerOrig[iColumn]));
04065 solver_->setColUpper(iColumn,min(upper[iColumn],upperOrig[iColumn]));
04066 }
04067 delete solver;
04068 }
04069 setCutoff(saveCutoff);
04070 return true;
04071 }
04072
04073
04074
04075
04076
04077 SbbModel *
04078 SbbModel::integerPresolve(bool weak)
04079 {
04080 status_ = 0;
04081
04082 bool feasible = resolve();
04083
04084 SbbModel * newModel = NULL;
04085 if (feasible) {
04086
04087
04088 newModel = new SbbModel(*this);
04089 newModel->messageHandler()->setLogLevel(messageHandler()->logLevel());
04090
04091 feasible = newModel->integerPresolveThisModel(solver_,weak);
04092 }
04093 if (!feasible) {
04094 handler_->message(SBB_INFEAS,messages_)
04095 <<CoinMessageEol;
04096 status_ = 1;
04097 delete newModel;
04098 return NULL;
04099 } else {
04100 newModel->synchronizeModel();
04101 return newModel;
04102 }
04103 }
04104
04105
04106
04107 bool
04108 SbbModel::integerPresolveThisModel(OsiSolverInterface * originalSolver,
04109 bool weak)
04110 {
04111 status_ = 0;
04112
04113 bool feasible = resolve();
04114
04115 bestObjective_=1.0e50;
04116 numberSolutions_=0;
04117 numberHeuristicSolutions_=0;
04118 double cutoff = getCutoff() ;
04119 double direction = solver_->getObjSense();
04120 if (cutoff < 1.0e20&&direction<0.0)
04121 messageHandler()->message(SBB_CUTOFF_WARNING1,
04122 messages())
04123 << cutoff << -cutoff << CoinMessageEol ;
04124 if (cutoff > bestObjective_)
04125 cutoff = bestObjective_ ;
04126 setCutoff(cutoff) ;
04127 int iColumn;
04128 int numberColumns = getNumCols();
04129 int originalNumberColumns = numberColumns;
04130 int numberPasses=0;
04131 synchronizeModel();
04132
04133 delete continuousSolver_;
04134 continuousSolver_ = solver_;
04135
04136 OsiSolverInterface * cleanModel = originalSolver->clone();
04137 #ifdef SBB_DEBUG
04138 std::string problemName;
04139 cleanModel->getStrParam(OsiProbName,problemName);
04140 printf("Problem name - %s\n",problemName.c_str());
04141 cleanModel->activateRowCutDebugger(problemName.c_str());
04142 const OsiRowCutDebugger * debugger = cleanModel->getRowCutDebugger();
04143 #endif
04144
04145
04146 int * original = new int[numberColumns];
04147
04148
04149 double * originalLower = new double[numberColumns];
04150 double * originalUpper = new double[numberColumns];
04151 {
04152 const double * lower = getColLower();
04153 const double * upper = getColUpper();
04154 for (iColumn=0;iColumn<numberColumns;iColumn++) {
04155 original[iColumn]=iColumn;
04156 originalLower[iColumn] = lower[iColumn];
04157 originalUpper[iColumn] = upper[iColumn];
04158 }
04159 }
04160 findIntegers(true);
04161
04162 int * originalPriority = NULL;
04163 int * originalIntegers = new int[numberIntegers_];
04164 int originalNumberIntegers = numberIntegers_;
04165 memcpy(originalIntegers,integerVariable_,numberIntegers_*sizeof(int));
04166
04167 if (priority_) {
04168 originalPriority = new int[numberIntegers_];
04169 memcpy(originalPriority,priority_,numberIntegers_*sizeof(int));
04170 delete [] priority_;
04171 priority_=NULL;
04172 }
04173 int todo=20;
04174 if (weak)
04175 todo=1;
04176 while (numberPasses<todo) {
04177
04178 numberPasses++;
04179 numberSolutions_=0;
04180
04181 bool doIntegerPresolve=(numberPasses!=20);
04182
04183
04184
04185 {
04186 const double * objective = cleanModel->getObjCoefficients();
04187 const double * lower = cleanModel->getColLower();
04188 const double * upper = cleanModel->getColUpper();
04189 double maximumCost=0.0;
04190 bool possibleMultiple=true;
04191 int numberChanged=0;
04192 for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
04193 if (originalUpper[iColumn]>originalLower[iColumn]) {
04194 if( cleanModel->isInteger(iColumn)) {
04195 maximumCost = max(maximumCost,fabs(objective[iColumn]));
04196 } else if (objective[iColumn]) {
04197 possibleMultiple=false;
04198 }
04199 }
04200 if (originalUpper[iColumn]<upper[iColumn]) {
04201 #ifdef SBB_DEBUG
04202 printf("Changing upper bound on %d from %g to %g\n",
04203 iColumn,upper[iColumn],originalUpper[iColumn]);
04204 #endif
04205 cleanModel->setColUpper(iColumn,originalUpper[iColumn]);
04206 numberChanged++;
04207 }
04208 if (originalLower[iColumn]>lower[iColumn]) {
04209 #ifdef SBB_DEBUG
04210 printf("Changing lower bound on %d from %g to %g\n",
04211 iColumn,lower[iColumn],originalLower[iColumn]);
04212 #endif
04213 cleanModel->setColLower(iColumn,originalLower[iColumn]);
04214 numberChanged++;
04215 }
04216 }
04217
04218 if (numberPasses==1)
04219 numberChanged += 1;
04220 if (possibleMultiple) {
04221 int increment=0;
04222 double multiplier = 2520.0;
04223 while (10.0*multiplier*maximumCost<1.0e8)
04224 multiplier *= 10.0;
04225 for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
04226 if (originalUpper[iColumn]>originalLower[iColumn]) {
04227 if( isInteger(iColumn)&&objective[iColumn]) {
04228 double value = fabs(objective[iColumn])*multiplier;
04229 int nearest = (int) floor(value+0.5);
04230 if (fabs(value-floor(value+0.5))>1.0e-8) {
04231 increment=0;
04232 break;
04233 } else if (!increment) {
04234
04235 increment=nearest;
04236 } else {
04237 increment = gcd(increment,nearest);
04238 }
04239 }
04240 }
04241 }
04242 if (increment) {
04243 double value = increment;
04244 value /= multiplier;
04245 if (value*0.999>dblParam_[SbbCutoffIncrement]) {
04246 messageHandler()->message(SBB_INTEGERINCREMENT,messages())
04247 <<value
04248 <<CoinMessageEol;
04249 dblParam_[SbbCutoffIncrement]=value*0.999;
04250 }
04251 }
04252 }
04253 if (!numberChanged) {
04254 doIntegerPresolve=false;
04255 }
04256 }
04257 #ifdef SBB_DEBUG
04258 if (debugger)
04259 assert(debugger->onOptimalPath(*cleanModel));
04260 #endif
04261 #ifdef COIN_USE_CLP
04262
04263 OsiClpSolverInterface * clpSolver
04264 = dynamic_cast<OsiClpSolverInterface *> (cleanModel);
04265 assert (clpSolver);
04266 ClpSimplex * clp = clpSolver->getModelPtr();
04267 ClpPresolve pinfo;
04268 ClpSimplex * model2 = pinfo.presolvedModel(*clp,1.0e-8);
04269 if (!model2) {
04270
04271 feasible=false;
04272 } else {
04273
04274 const int * originalColumns = pinfo.originalColumns();
04275
04276 OsiClpSolverInterface * temp = new OsiClpSolverInterface(model2);
04277 numberColumns = temp->getNumCols();
04278 for (iColumn=0;iColumn<originalNumberColumns;iColumn++)
04279 original[iColumn]=-1;
04280 for (iColumn=0;iColumn<numberColumns;iColumn++)
04281 original[originalColumns[iColumn]]=iColumn;
04282
04283 temp->copyParameters(*solver_);
04284 delete solver_;
04285 solver_ = temp;
04286 setCutoff(cutoff);
04287 deleteObjects();
04288 synchronizeModel();
04289
04290 continuousSolver_ = solver_;
04291 feasible=resolve();
04292 if (!feasible||!doIntegerPresolve||weak) break;
04293
04294 int found=-1;
04295 int iHeuristic;
04296 double * newSolution = new double [numberColumns];
04297 double heuristicValue=getCutoff();
04298 for (iHeuristic=0;iHeuristic<numberHeuristics_;iHeuristic++) {
04299 double saveValue=heuristicValue;
04300 int ifSol = heuristic_[iHeuristic]->solution(heuristicValue,
04301 newSolution);
04302 if (ifSol>0) {
04303
04304 found=iHeuristic;
04305 } else if (ifSol<0) {
04306 heuristicValue = saveValue;
04307 }
04308 }
04309 if (found >= 0) {
04310
04311 int numberColumns = getNumCols() ;
04312 if (!currentSolution_)
04313 currentSolution_ = new double[numberColumns] ;
04314
04315 setBestSolution(SBB_ROUNDING,heuristicValue,
04316 newSolution);
04317
04318 cutoff = getCutoff();
04319 }
04320 delete [] newSolution;
04321
04322 int maximumWhich=1000;
04323 int * whichGenerator = new int[maximumWhich];
04324
04325 numberRowsAtContinuous_ = getNumRows();
04326 maximumNumberCuts_=0;
04327 currentNumberCuts_=0;
04328 delete [] addedCuts_;
04329 addedCuts_ = NULL;
04330
04331
04332 maximumDepth_=10;
04333 delete [] walkback_;
04334 walkback_ = new SbbNodeInfo * [maximumDepth_];
04335
04336 OsiCuts cuts;
04337 int numberOldActiveCuts=0;
04338 int numberNewCuts = 0;
04339 feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
04340 NULL,numberOldActiveCuts,numberNewCuts,
04341 maximumWhich, whichGenerator);
04342 currentNumberCuts_=numberNewCuts;
04343 delete [] whichGenerator;
04344 delete [] walkback_;
04345 walkback_ = NULL;
04346 delete [] addedCuts_;
04347 addedCuts_=NULL;
04348 if (feasible) {
04349
04350
04351 const double * lower = solver_->getColLower();
04352 const double * upper = solver_->getColUpper();
04353 int i;
04354 for (i=0;i<originalNumberIntegers;i++) {
04355 iColumn = originalIntegers[i];
04356 int jColumn = original[iColumn];
04357 if (jColumn >= 0) {
04358 if (upper[jColumn]<originalUpper[iColumn])
04359 originalUpper[iColumn] = upper[jColumn];
04360 if (lower[jColumn]>originalLower[iColumn])
04361 originalLower[iColumn] = lower[jColumn];
04362 }
04363 }
04364 }
04365 }
04366 #endif
04367 if (!feasible||!doIntegerPresolve) {
04368 break;
04369 }
04370 }
04371
04372 delete cleanModel;
04373
04374 if (originalPriority) {
04375 priority_ = new int[numberIntegers_];
04376 int i;
04377 int number=0;
04378
04379 for (i=0;i<originalNumberIntegers;i++) {
04380 iColumn = originalIntegers[i];
04381 int newColumn=original[iColumn];
04382 if (newColumn >= 0)
04383 priority_[number++]=originalPriority[i];
04384 }
04385 assert (number==numberIntegers_);
04386 delete [] originalPriority;
04387 }
04388 delete [] originalIntegers;
04389 numberColumns = getNumCols();
04390 delete [] originalColumns_;
04391 originalColumns_ = new int[numberColumns];
04392 numberColumns=0;
04393 for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
04394 int jColumn = original[iColumn];
04395 if (jColumn >= 0)
04396 originalColumns_[numberColumns++]=iColumn;
04397 }
04398 delete [] original;
04399 delete [] originalLower;
04400 delete [] originalUpper;
04401
04402 deleteObjects();
04403 synchronizeModel();
04404 continuousSolver_=NULL;
04405 currentNumberCuts_=0;
04406 return feasible;
04407 }
04408
04409 void
04410 SbbModel::originalModel(SbbModel * presolvedModel,bool weak)
04411 {
04412 solver_->copyParameters(*(presolvedModel->solver_));
04413 bestObjective_ = presolvedModel->bestObjective_;
04414 delete [] bestSolution_;
04415 findIntegers(true);
04416 if (presolvedModel->bestSolution_) {
04417 int numberColumns = getNumCols();
04418 int numberOtherColumns = presolvedModel->getNumCols();
04419 bestSolution_ = new double[numberColumns];
04420
04421 int * back = new int[numberColumns];
04422 int i;
04423 for (i=0;i<numberColumns;i++)
04424 back[i]=-1;
04425 for (i=0;i<numberOtherColumns;i++)
04426 back[presolvedModel->originalColumns_[i]]=i;
04427 int iColumn;
04428
04429 double * otherSolution = presolvedModel->bestSolution_;
04430 for (i=0;i<numberIntegers_;i++) {
04431 iColumn = integerVariable_[i];
04432 int jColumn = back[iColumn];
04433 if (jColumn >= 0) {
04434 double value=floor(otherSolution[jColumn]+0.5);
04435 solver_->setColLower(iColumn,value);
04436 solver_->setColUpper(iColumn,value);
04437 }
04438 }
04439 #if 0
04440
04441
04442 OsiClpSolverInterface * clpSolver
04443 = dynamic_cast<OsiClpSolverInterface *> (solver_);
04444 assert (clpSolver);
04445 ClpSimplex * clp = clpSolver->getModelPtr();
04446 Presolve pinfo;
04447 ClpSimplex * model2 = pinfo.presolvedModel(*clp,1.0e-8);
04448 model2->primal(1);
04449 pinfo.postsolve(true);
04450 const double * solution = solver_->getColSolution();
04451 for (i=0;i<numberIntegers_;i++) {
04452 iColumn = integerVariable_[i];
04453 double value=floor(solution[iColumn]+0.5);
04454 solver_->setColLower(iColumn,value);
04455 solver_->setColUpper(iColumn,value);
04456 }
04457 #else
04458 if (!weak) {
04459
04460 int save = numberCutGenerators_;
04461 numberCutGenerators_=0;
04462 branchAndBound();
04463 numberCutGenerators_=save;
04464 }
04465 #endif
04466
04467 resolve();
04468
04469 int numberIntegerInfeasibilities;
04470 int numberObjectInfeasibilities;
04471 if (!currentSolution_)
04472 currentSolution_ = new double[numberColumns] ;
04473 assert(feasibleSolution(numberIntegerInfeasibilities,
04474 numberObjectInfeasibilities));
04475 } else {
04476 bestSolution_=NULL;
04477 }
04478 numberSolutions_=presolvedModel->numberSolutions_;
04479 numberHeuristicSolutions_=presolvedModel->numberHeuristicSolutions_;
04480 numberNodes_ = presolvedModel->numberNodes_;
04481 numberIterations_ = presolvedModel->numberIterations_;
04482 status_ = presolvedModel->status_;
04483 synchronizeModel();
04484 }
04485
04486 void
04487 SbbModel::passInMessageHandler(CoinMessageHandler * handler)
04488 {
04489 if (defaultHandler_) {
04490 delete handler_;
04491 handler_ = NULL;
04492 }
04493 defaultHandler_=false;
04494 handler_=handler;
04495 }
04496