source: trunk/CbcTreeLocal.cpp @ 190

Last change on this file since 190 was 190, checked in by forrest, 14 years ago

local search

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.3 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CbcModel.hpp"
5#include "CbcNode.hpp"
6#include "CbcTreeLocal.hpp"
7#include "CoinPackedMatrix.hpp"
8#include "CoinTime.hpp"
9#include <cassert>
10#if 0
11// gdb doesn't always put breakpoints in this virtual function
12// just stick xxxxxx() where you want to start
13static void xxxxxx()
14{
15  printf("break\n");
16}
17#endif
18CbcTreeLocal::CbcTreeLocal()
19  : localNode_(NULL),
20    bestSolution_(NULL),
21    savedSolution_(NULL),
22    saveNumberSolutions_(0),
23    model_(NULL),
24    originalLower_(NULL),
25    originalUpper_(NULL),
26    range_(0),
27    typeCuts_(-1),
28    maxDiversification_(0),
29    diversification_(0),
30    nextStrong_(false),
31    rhs_(0.0),
32    savedGap_(0.0),
33    bestCutoff_(0.0),
34    timeLimit_(0),
35    startTime_(0),
36    nodeLimit_(0),
37    startNode_(-1),
38    searchType_(-1),
39    refine_(false)
40{
41 
42}
43/* Constructor with solution.
44   range is upper bound on difference from given solution.
45   maxDiversification is maximum number of diversifications to try
46   timeLimit is seconds in subTree
47   nodeLimit is nodes in subTree
48*/
49CbcTreeLocal::CbcTreeLocal(CbcModel * model,const double * solution ,
50                                  int range,int typeCuts,int maxDiversification,
51                                  int timeLimit, int nodeLimit,bool refine)
52  : localNode_(NULL),
53    bestSolution_(NULL),
54    savedSolution_(NULL),
55    saveNumberSolutions_(0),
56    model_(model),
57    originalLower_(NULL),
58    originalUpper_(NULL),
59    range_(range),
60    typeCuts_(typeCuts),
61    maxDiversification_(maxDiversification),
62    diversification_(0),
63    nextStrong_(false),
64    rhs_(0.0),
65    savedGap_(0.0),
66    bestCutoff_(0.0),
67    timeLimit_(timeLimit),
68    startTime_(0),
69    nodeLimit_(nodeLimit),
70    startNode_(-1),
71    searchType_(-1),
72    refine_(refine)
73{
74
75  OsiSolverInterface * solver = model_->solver();
76  const double * lower = solver->getColLower();
77  const double * upper = solver->getColUpper();
78  //const double * solution = solver->getColSolution();
79  //const double * objective = solver->getObjCoefficients();
80  double primalTolerance;
81  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
82
83  // Get increment
84  model_->analyzeObjective();
85
86  bestCutoff_ = model_->getCutoff();
87  // save current gap
88  savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
89
90  // make sure integers found
91  model_->findIntegers(false);
92  int numberIntegers = model_->numberIntegers();
93  const int * integerVariable = model_->integerVariable();
94  int i;
95  double direction = solver->getObjSense();
96  double newSolutionValue=1.0e50;
97  if (solution) {
98    // copy solution
99    solver->setColSolution(solution);
100    newSolutionValue = direction*solver->getObjValue();
101  }
102  originalLower_=new double [numberIntegers];
103  originalUpper_=new double [numberIntegers];
104  bool all01=true;
105  int number01=0;
106  for (i=0;i<numberIntegers;i++) {
107    int iColumn = integerVariable[i];
108    originalLower_[i]=lower[iColumn];
109    originalUpper_[i]=upper[iColumn];
110    if (upper[iColumn]-lower[iColumn]>1.5)
111      all01=false;
112    else if (upper[iColumn]-lower[iColumn]==1.0)
113      number01++;
114  }
115  if (all01 && !typeCuts_)
116    typeCuts_=1; // may as well so we don't have to deal with refine
117  if (!number01&&!typeCuts_) {
118    printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
119    typeCuts_=-1;
120  } else {
121    std::string type;
122    if (all01) {
123      printf("%d 0-1 variables normal local  cuts\n",
124             number01);
125    } else if (typeCuts_) {
126      printf("%d 0-1 variables, %d other - general integer local cuts\n",
127             number01,numberIntegers-number01);
128    } else {
129      printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
130             number01,numberIntegers-number01);
131    }
132    printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
133           maxDiversification_,range_,timeLimit_,nodeLimit_);
134  }
135  int numberColumns = model_->getNumCols();
136  savedSolution_ = new double [numberColumns];
137  memset(savedSolution_,0,numberColumns*sizeof(double));
138  if (solution) {
139    rhs_=range_;
140    // Check feasible
141    int goodSolution=createCut(solution,cut_);
142    if (goodSolution>=0) {
143      for (i=0;i<numberIntegers;i++) {
144        int iColumn = integerVariable[i];
145        double value = floor(solution[iColumn]+0.5);
146        // fix so setBestSolution will work
147        solver->setColLower(iColumn,value);
148        solver->setColUpper(iColumn,value);
149      }
150      model_->reserveCurrentSolution();
151      // Create cut and get total gap
152      if (newSolutionValue<bestCutoff_) {
153        model_->setBestSolution(CBC_ROUNDING,newSolutionValue,solution,
154                                false);
155        bestCutoff_=newSolutionValue;
156        // save as best solution
157        memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
158      }
159      for (i=0;i<numberIntegers;i++) {
160        int iColumn = integerVariable[i];
161        // restore bounds
162        solver->setColLower(iColumn,originalLower_[i]);
163        solver->setColUpper(iColumn,originalUpper_[i]);
164      }
165      // make sure can't stop on gap
166      model_->setDblParam(CbcModel::CbcAllowableGap,-1.0e50);
167    } else {
168      model_=NULL;
169    }
170  } else {
171    // no solution
172    rhs_=1.0e50;
173    // make sure can't stop on gap
174    model_->setDblParam(CbcModel::CbcAllowableGap,-1.0e50);
175  }
176}
177CbcTreeLocal::~CbcTreeLocal()
178{
179  delete [] originalLower_;
180  delete [] originalUpper_;
181  delete [] bestSolution_;
182  delete [] savedSolution_;
183  delete localNode_;
184}
185// Copy constructor
186CbcTreeLocal::CbcTreeLocal ( const CbcTreeLocal & rhs)
187  :CbcTree(rhs),
188   saveNumberSolutions_(rhs.saveNumberSolutions_),
189   model_(rhs.model_),
190   range_(rhs.range_),
191   typeCuts_(rhs.typeCuts_),
192   maxDiversification_(rhs.maxDiversification_),
193   diversification_(rhs.diversification_),
194   nextStrong_(rhs.nextStrong_),
195   rhs_(rhs.rhs_),
196   savedGap_(rhs.savedGap_),
197   bestCutoff_(rhs.bestCutoff_),
198   timeLimit_(rhs.timeLimit_),
199   startTime_(rhs.startTime_),
200   nodeLimit_(rhs.nodeLimit_),
201   startNode_(rhs.startNode_),
202   searchType_(rhs.searchType_),
203   refine_(rhs.refine_)
204{ 
205  cut_=rhs.cut_;
206  fixedCut_=rhs.fixedCut_;
207  if (rhs.localNode_)
208    localNode_=new CbcNode(*rhs.localNode_);
209  else
210    localNode_=NULL;
211  if (rhs.originalLower_) {
212    int numberIntegers = model_->numberIntegers();
213    originalLower_=new double [numberIntegers];
214    memcpy(originalLower_,rhs.originalLower_,numberIntegers*sizeof(double));
215    originalUpper_=new double [numberIntegers];
216    memcpy(originalUpper_,rhs.originalUpper_,numberIntegers*sizeof(double));
217  } else {
218    originalLower_=NULL;
219    originalUpper_=NULL;
220  }
221  if (rhs.bestSolution_) {
222    int numberColumns = model_->getNumCols();
223    bestSolution_ = new double [numberColumns];
224    memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
225  } else {
226    bestSolution_=NULL;
227  }
228  if (rhs.savedSolution_) {
229    int numberColumns = model_->getNumCols();
230    savedSolution_ = new double [numberColumns];
231    memcpy(savedSolution_,rhs.savedSolution_,numberColumns*sizeof(double));
232  } else {
233    savedSolution_=NULL;
234  }
235}
236//----------------------------------------------------------------
237// Assignment operator
238//-------------------------------------------------------------------
239CbcTreeLocal &
240CbcTreeLocal::operator=(const CbcTreeLocal& rhs)
241{
242  if (this != &rhs) {
243    CbcTree::operator=(rhs);
244    saveNumberSolutions_ = rhs.saveNumberSolutions_;
245    cut_=rhs.cut_;
246    fixedCut_=rhs.fixedCut_;
247    delete localNode_;
248    if (rhs.localNode_)
249      localNode_=new CbcNode(*rhs.localNode_);
250    else
251      localNode_=NULL;
252    model_ = rhs.model_;
253    range_ = rhs.range_;
254    typeCuts_ = rhs.typeCuts_;
255    maxDiversification_ = rhs.maxDiversification_;
256    diversification_ = rhs.diversification_;
257    nextStrong_ = rhs.nextStrong_;
258    rhs_ = rhs.rhs_;
259    savedGap_ = rhs.savedGap_;
260    bestCutoff_ = rhs.bestCutoff_;
261    timeLimit_ = rhs.timeLimit_;
262    startTime_ = rhs.startTime_;
263    nodeLimit_ = rhs.nodeLimit_;
264    startNode_ = rhs.startNode_;
265    searchType_ = rhs.searchType_;
266    refine_ = rhs.refine_;
267    delete [] originalLower_;
268    delete [] originalUpper_;
269    if (rhs.originalLower_) {
270      int numberIntegers = model_->numberIntegers();
271      originalLower_=new double [numberIntegers];
272      memcpy(originalLower_,rhs.originalLower_,numberIntegers*sizeof(double));
273      originalUpper_=new double [numberIntegers];
274      memcpy(originalUpper_,rhs.originalUpper_,numberIntegers*sizeof(double));
275    } else {
276      originalLower_=NULL;
277      originalUpper_=NULL;
278    }
279    delete [] bestSolution_;
280    if (rhs.bestSolution_) {
281      int numberColumns = model_->getNumCols();
282      bestSolution_ = new double [numberColumns];
283      memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
284    } else {
285      bestSolution_=NULL;
286    }
287    delete [] savedSolution_;
288    if (rhs.savedSolution_) {
289      int numberColumns = model_->getNumCols();
290      savedSolution_ = new double [numberColumns];
291      memcpy(savedSolution_,rhs.savedSolution_,numberColumns*sizeof(double));
292    } else {
293      savedSolution_=NULL;
294    }
295  }
296  return *this;
297}
298// Clone
299CbcTree *
300CbcTreeLocal::clone() const
301{
302  return new CbcTreeLocal(*this);
303}
304// Pass in solution (so can be used after heuristic)
305void 
306CbcTreeLocal::passInSolution(const double * solution, double solutionValue)
307{
308  int numberColumns = model_->getNumCols();
309  delete [] savedSolution_;
310  savedSolution_ = new double [numberColumns];
311  memcpy(savedSolution_,solution,numberColumns*sizeof(double));
312  rhs_=range_;
313  // Check feasible
314  int goodSolution=createCut(solution,cut_);
315  if (goodSolution>=0) {
316    bestCutoff_=solutionValue;
317  } else {
318    model_=NULL;
319  }
320}
321// Return the top node of the heap
322CbcNode * 
323CbcTreeLocal::top() {
324#ifdef CBC_DEBUG
325  int smallest=9999999;
326  int largest=-1;
327  double smallestD=1.0e30;
328  double largestD=-1.0e30;
329  int n=nodes_.size();
330  for (int i=0;i<n;i++) {
331    int nn=nodes_[i]->nodeNumber();
332    double dd = nodes_[i]->objectiveValue();
333    largest=max(largest,nn);
334    smallest=min(smallest,nn);
335    largestD=max(largestD,dd);
336    smallestD=min(smallestD,dd);
337  }
338  printf("smallest %d, largest %d, top %d\n",smallest,largest,nodes_.front()->nodeNumber());
339  printf("smallestD %g, largestD %g, top %g\n",smallestD,largestD,nodes_.front()->objectiveValue());
340#endif
341  return nodes_.front();
342}
343
344// Add a node to the heap
345void 
346CbcTreeLocal::push(CbcNode * x) {
347  if (typeCuts_>=0&&!nodes_.size()&&searchType_<0) {
348    startNode_=model_->getNodeCount();
349    // save copy of node
350    localNode_ = new CbcNode(*x);
351
352    if (cut_.row().getNumElements()) {
353      // Add to global cuts
354      // we came in with solution
355      model_->globalCuts()->insert(cut_);
356      printf("initial cut - rhs %g %g\n",
357             cut_.lb(),cut_.ub());
358      searchType_=1;
359    } else {
360      // stop on first solution
361      searchType_=0;
362    }
363    startTime_ = (int) CoinCpuTime();
364    saveNumberSolutions_ = model_->getSolutionCount();
365  }
366  nodes_.push_back(x);
367#ifdef CBC_DEBUG
368  printf("pushing node onto heap %d %x %x\n",x->nodeNumber(),x,x->nodeInfo());
369#endif
370  push_heap(nodes_.begin(), nodes_.end(), comparison_);
371}
372
373// Remove the top node from the heap
374void 
375CbcTreeLocal::pop() {
376  pop_heap(nodes_.begin(), nodes_.end(), comparison_);
377  nodes_.pop_back();
378}
379// Test if empty - does work if so
380bool 
381CbcTreeLocal::empty()   
382{
383  if (typeCuts_<0) 
384    return !nodes_.size();
385  /* state -
386     0 iterating
387     1 subtree finished optimal solution for subtree found
388     2 subtree finished and no solution found
389     3 subtree exiting and solution found
390     4 subtree exiting and no solution found
391  */
392  int state=0;
393  assert (searchType_!=2);
394  if (searchType_) {
395    if (CoinCpuTime()-startTime_>timeLimit_||model_->getNodeCount()-startNode_>=nodeLimit_) {
396      state=4;
397    }
398  } else {
399    if (model_->getSolutionCount()>saveNumberSolutions_) {
400      state=4;
401    }
402  }
403  if (!nodes_.size()) 
404    state=2;
405  if (!state) {
406    return false;
407  }
408  // Finished this phase
409  int numberColumns = model_->getNumCols();
410  if (model_->getSolutionCount()>saveNumberSolutions_) {
411    if(model_->getCutoff()<bestCutoff_) {
412      // Save solution
413      if (!bestSolution_) 
414        bestSolution_ = new double [numberColumns];
415      memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
416      bestCutoff_=model_->getCutoff();
417    }
418    state--;
419  } 
420  // get rid of all nodes (safe even if already done)
421  double bestPossibleObjective;
422  cleanTree(model_,-COIN_DBL_MAX,bestPossibleObjective);
423
424  double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
425  printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
426         state,
427         model_->getNodeCount()-startNode_,
428         (int) CoinCpuTime()-startTime_,
429         model_->getCutoff()+increment,bestCutoff_+increment,rhs_);
430  saveNumberSolutions_ = model_->getSolutionCount();
431  bool finished=false;
432  bool lastTry=false;
433  switch (state) {
434  case 1:
435    // solution found and subtree exhausted
436    if (rhs_>1.0e30) {
437      finished=true;
438    } else {
439      // find global cut and reverse
440      reverseCut(1);
441      searchType_=1; // first false
442      rhs_ = range_; // reset range
443      nextStrong_=false;
444
445      // save best solution in this subtree
446      memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
447    }
448    break;
449  case 2:
450    // solution not found and subtree exhausted
451    if (rhs_>1.0e30) {
452      finished=true;
453    } else {
454      // find global cut and reverse
455      reverseCut(2);
456      searchType_=1; // first false
457      if (diversification_<maxDiversification_) {
458        if (nextStrong_) {
459          diversification_++;
460          model_->setCutoff(1.0e50);
461          searchType_=0;
462        }
463        nextStrong_=true;
464        rhs_ += range_/2;
465      } else {
466        // This will be last try (may hit max time)
467        lastTry=true;
468        if (!maxDiversification_)
469          typeCuts_=-1; // make sure can't start again
470        model_->setCutoff(bestCutoff_);
471        printf("Exiting local search with current set of cuts\n");
472        rhs_=1.0e100;
473        // Can now stop on gap
474        model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
475      }
476    }
477    break;
478  case 3:
479    // solution found and subtree not exhausted
480    if (rhs_<1.0e30) {
481      if (searchType_) {
482        if (!typeCuts_&&refine_&&searchType_==1) {
483          // We need to check we have best solution given these 0-1 values
484          OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
485          CbcModel * subModel = model_->subTreeModel(subSolver);
486          CbcTree normalTree;
487          subModel->passInTreeHandler(normalTree);
488          int numberIntegers = model_->numberIntegers();
489          const int * integerVariable = model_->integerVariable();
490          const double * solution = model_->bestSolution();
491          int i;
492          int numberColumns = model_->getNumCols();
493          for (i=0;i<numberIntegers;i++) {
494            int iColumn = integerVariable[i];
495            double value = floor(solution[iColumn]+0.5);
496            if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
497              continue; // skip as not 0-1
498            if (originalLower_[i]==originalUpper_[i])
499              continue;
500            subSolver->setColLower(iColumn,value);
501            subSolver->setColUpper(iColumn,value);
502          }
503          subSolver->initialSolve();
504          // We can copy cutoff
505          // But adjust
506          subModel->setCutoff(model_->getCutoff()+model_->getDblParam(CbcModel::CbcCutoffIncrement)+1.0e-6);
507          subModel->setSolutionCount(0);
508          assert (subModel->isProvenOptimal());
509          if (!subModel->typePresolve()) {
510            subModel->branchAndBound();
511            if (subModel->status()) {
512              model_->incrementSubTreeStopped(); 
513            }
514            //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
515            //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
516            double newCutoff = subModel->getMinimizationObjValue()-
517              subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
518            if (subModel->getSolutionCount()) {
519              if (!subModel->status())
520                assert (subModel->isProvenOptimal());
521              memcpy(model_->bestSolution(),subModel->bestSolution(),
522                     numberColumns*sizeof(double));
523              model_->setCutoff(newCutoff);
524            }
525          } else if (subModel->typePresolve()==1) {
526            CbcModel * model2 = subModel->integerPresolve(true);
527            if (model2) {
528              // Do complete search
529              model2->branchAndBound();
530              // get back solution
531              subModel->originalModel(model2,false);
532              if (model2->status()) {
533                model_->incrementSubTreeStopped(); 
534              } 
535              double newCutoff = model2->getMinimizationObjValue()-
536                model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
537              if (model2->getSolutionCount()) {
538                if (!model2->status())
539                  assert (model2->isProvenOptimal());
540                memcpy(model_->bestSolution(),subModel->bestSolution(),
541                       numberColumns*sizeof(double));
542                model_->setCutoff(newCutoff);
543              }
544              delete model2;
545            } else {
546              // infeasible - could just be - due to cutoff
547            }
548          } else {
549            // too dangerous at present
550            assert (subModel->typePresolve()!=2);
551          }
552          if(model_->getCutoff()<bestCutoff_) {
553            // Save solution
554            if (!bestSolution_) 
555              bestSolution_ = new double [numberColumns];
556            memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
557            bestCutoff_=model_->getCutoff();
558          }
559          delete subModel;
560        }
561        // we have done search to make sure best general solution
562        searchType_=1;
563        // Reverse cut weakly
564        reverseCut(3,rhs_);
565      } else {
566        searchType_=1;
567        // delete last cut
568        deleteCut(cut_);
569      }
570    } else {
571      searchType_=1;
572    }
573    // save best solution in this subtree
574    memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
575    nextStrong_=false;
576    rhs_=range_;
577    break;
578  case 4:
579    // solution not found and subtree not exhausted
580    if (maxDiversification_) {
581      if (nextStrong_) {
582        // Reverse cut weakly
583        reverseCut(4,rhs_);
584        model_->setCutoff(1.0e50);
585        diversification_++;
586        searchType_=0;
587      } else {
588        // delete last cut
589        deleteCut(cut_);
590        searchType_=1;
591      }
592      nextStrong_=true;
593      rhs_ += range_/2;
594    } else {
595      // special case when using as heuristic
596      // This will be last try (may hit max time0
597      lastTry=true;
598      model_->setCutoff(bestCutoff_);
599      printf("Exiting local search with current set of cuts\n");
600      rhs_=1.0e100;
601      // Can now stop on gap
602      model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
603      typeCuts_ =-1;
604    }
605    break;
606  }
607  if (rhs_<1.0e30||lastTry) {
608    int goodSolution=createCut(savedSolution_,cut_);
609    if (goodSolution>=0) {
610      // Add to global cuts
611      model_->globalCuts()->insert(cut_);
612      OsiCuts * global = model_->globalCuts();
613      int n = global->sizeRowCuts();
614      OsiRowCut * rowCut = global->rowCutPtr(n-1);
615      printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
616             n,rowCut->lb(),rowCut->ub(),rhs_,diversification_);
617      for (int i=0;i<n;i++) {
618        rowCut = global->rowCutPtr(i);
619        printf("%d - rhs %g %g\n",
620               i,rowCut->lb(),rowCut->ub());
621      }
622    }
623    // put back node
624    startTime_ = (int) CoinCpuTime();
625    startNode_=model_->getNodeCount();
626    if (localNode_) {
627      // save copy of node
628      CbcNode * localNode2 = new CbcNode(*localNode_);
629      // But localNode2 now owns cuts so swap
630      //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
631      //   localNode_,localNode_->nodeInfo());
632      nodes_.push_back(localNode_);
633      localNode_=localNode2;
634      make_heap(nodes_.begin(), nodes_.end(), comparison_);
635    }
636  }
637  return finished;
638}
639// We may have got an intelligent tree so give it one more chance
640void 
641CbcTreeLocal::endSearch() 
642{
643  if (typeCuts_>=0) {
644    // copy best solution to model
645    int numberColumns = model_->getNumCols();
646    if (bestSolution_&&bestCutoff_<model_->getCutoff()) {
647      memcpy(model_->bestSolution(),bestSolution_,numberColumns*sizeof(double));
648      model_->setCutoff(bestCutoff_);
649      // recompute objective value
650      const double * objCoef = model_->getObjCoefficients();
651      double objOffset=0.0;
652      model_->continuousSolver()->getDblParam(OsiObjOffset,objOffset);
653 
654      // Compute dot product of objCoef and colSol and then adjust by offset
655      double objValue = -objOffset;
656      for ( int i=0 ; i<numberColumns ; i++ )
657        objValue += objCoef[i]*bestSolution_[i];
658      model_->setMinimizationObjValue(objValue);
659    }
660    // Can now stop on gap
661    model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
662  }
663}
664// Create cut
665int
666CbcTreeLocal::createCut(const double * solution,OsiRowCut & rowCut)
667{
668  if (rhs_>1.0e20)
669    return -1;
670  OsiSolverInterface * solver = model_->solver();
671  const double * rowLower = solver->getRowLower();
672  const double * rowUpper = solver->getRowUpper();
673  //const double * solution = solver->getColSolution();
674  //const double * objective = solver->getObjCoefficients();
675  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
676  double primalTolerance;
677  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
678  // relax
679  primalTolerance *= 1000.0;
680
681  int numberRows = model_->getNumRows();
682
683  int numberIntegers = model_->numberIntegers();
684  const int * integerVariable = model_->integerVariable();
685  int i;
686
687  // Check feasible
688
689  double * rowActivity = new double[numberRows];
690  memset(rowActivity,0,numberRows*sizeof(double));
691  solver->getMatrixByCol()->times(solution,rowActivity) ;
692  int goodSolution=0;
693  // check was feasible
694  for (i=0;i<numberRows;i++) {
695    if(rowActivity[i]<rowLower[i]-primalTolerance) {
696      goodSolution=-1;
697    } else if(rowActivity[i]>rowUpper[i]+primalTolerance) {
698      goodSolution=-1;
699    }
700  }
701  delete [] rowActivity;
702  for (i=0;i<numberIntegers;i++) {
703    int iColumn = integerVariable[i];
704    double value=solution[iColumn];
705    if (fabs(floor(value+0.5)-value)>integerTolerance) {
706      goodSolution=-1;
707    }
708  }
709  // zap cut
710  if (goodSolution==0) {
711    // Create cut and get total gap
712    CoinPackedVector cut;
713    double rhs = rhs_;
714    double maxValue=0.0;
715    for (i=0;i<numberIntegers;i++) {
716      int iColumn = integerVariable[i];
717      double value = floor(solution[iColumn]+0.5);
718      if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
719        continue; // skip as not 0-1
720      if (originalLower_[i]==originalUpper_[i])
721        continue;
722      double mu =1.0/(originalUpper_[i]-originalLower_[i]);
723      if (value==originalLower_[i]) {
724        rhs += mu*originalLower_[i];
725        cut.insert(iColumn,1.0);
726        maxValue += originalUpper_[i];
727      } else if (value==originalUpper_[i]) {
728        rhs -= mu*originalUpper_[i];
729        cut.insert(iColumn,-1.0);
730        maxValue += originalLower_[i];
731      }
732    }
733    if (maxValue<rhs-primalTolerance) {
734      printf("slack cut\n");
735      goodSolution=1;
736    }
737    rowCut.setRow(cut);
738    rowCut.setLb(-COIN_DBL_MAX);
739    rowCut.setUb(rhs);   
740    rowCut.setGloballyValid();
741    printf("Cut size: %i Cut rhs: %g\n",cut.getNumElements(), rhs);
742#ifdef CBC_DEBUG
743    int k;
744    for (k=0; k<cut.getNumElements(); k++){
745      printf("%i    %g ",cut.getIndices()[k], cut.getElements()[k]);
746      if ((k+1)%5==0)
747        printf("\n");
748    }
749    if (k%5!=0)
750      printf("\n");
751#endif
752    return goodSolution;
753  } else {
754    printf("Not a good solution\n");
755    return -1;
756  }
757}
758// Other side of last cut branch
759void 
760CbcTreeLocal::reverseCut(int state, double bias)
761{
762  // find global cut
763  OsiCuts * global = model_->globalCuts();
764  int n = global->sizeRowCuts();
765  int i;
766  OsiRowCut * rowCut = NULL;
767  for ( i=0;i<n;i++) {
768    rowCut = global->rowCutPtr(i);
769    if (cut_==*rowCut) {
770      break;
771    }
772  }
773  if (!rowCut) {
774    // must have got here in odd way e.g. strong branching
775    return;
776  }
777  // get smallest element
778  double smallest=COIN_DBL_MAX;
779  CoinPackedVector row = cut_.row();
780  for (int k=0; k<row.getNumElements(); k++)
781    smallest = min(smallest,fabs(row.getElements()[k]));
782  if (!typeCuts_&&!refine_) {
783    // Reverse cut very very weakly
784    if (state>2)
785      smallest=0.0;
786  }
787  // replace by other way
788  printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
789         i,n,rowCut->lb(),rowCut->ub());
790  rowCut->setLb(rowCut->ub()+smallest-bias);
791  rowCut->setUb(COIN_DBL_MAX);
792  printf("new rhs %g %g, bias %g smallest %g ",
793         rowCut->lb(),rowCut->ub(),bias,smallest);
794}
795// Delete last cut branch
796void 
797CbcTreeLocal::deleteCut(OsiRowCut & cut)
798{
799  // find global cut
800  OsiCuts * global = model_->globalCuts();
801  int n = global->sizeRowCuts();
802  int i;
803  OsiRowCut * rowCut = NULL;
804  for ( i=0;i<n;i++) {
805    rowCut = global->rowCutPtr(i);
806    if (cut==*rowCut) {
807      break;
808    }
809  }
810  assert (i<n);
811  // delete last cut
812  printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
813         i,n,rowCut->lb(),rowCut->ub());
814  global->eraseRowCut(i);
815}
816
817
818
Note: See TracBrowser for help on using the repository browser.