source: trunk/CbcTreeLocal.cpp @ 185

Last change on this file since 185 was 185, checked in by forrest, 16 years ago

move CbcTreeLocal?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.0 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// Return the top node of the heap
305CbcNode * 
306CbcTreeLocal::top() {
307#ifdef CBC_DEBUG
308  int smallest=9999999;
309  int largest=-1;
310  double smallestD=1.0e30;
311  double largestD=-1.0e30;
312  int n=nodes_.size();
313  for (int i=0;i<n;i++) {
314    int nn=nodes_[i]->nodeNumber();
315    double dd = nodes_[i]->objectiveValue();
316    largest=max(largest,nn);
317    smallest=min(smallest,nn);
318    largestD=max(largestD,dd);
319    smallestD=min(smallestD,dd);
320  }
321  printf("smallest %d, largest %d, top %d\n",smallest,largest,nodes_.front()->nodeNumber());
322  printf("smallestD %g, largestD %g, top %g\n",smallestD,largestD,nodes_.front()->objectiveValue());
323#endif
324  return nodes_.front();
325}
326
327// Add a node to the heap
328void 
329CbcTreeLocal::push(CbcNode * x) {
330  if (typeCuts_>=0&&!nodes_.size()&&searchType_<0) {
331    startNode_=model_->getNodeCount();
332    // save copy of node
333    localNode_ = new CbcNode(*x);
334
335    if (cut_.row().getNumElements()) {
336      // Add to global cuts
337      // we came in with solution
338      model_->globalCuts()->insert(cut_);
339      printf("initial cut - rhs %g %g\n",
340             cut_.lb(),cut_.ub());
341      searchType_=1;
342    } else {
343      // stop on first solution
344      searchType_=0;
345    }
346    startTime_ = (int) CoinCpuTime();
347    saveNumberSolutions_ = model_->getSolutionCount();
348  }
349  nodes_.push_back(x);
350#ifdef CBC_DEBUG
351  printf("pushing node onto heap %d %x %x\n",x->nodeNumber(),x,x->nodeInfo());
352#endif
353  push_heap(nodes_.begin(), nodes_.end(), comparison_);
354}
355
356// Remove the top node from the heap
357void 
358CbcTreeLocal::pop() {
359  pop_heap(nodes_.begin(), nodes_.end(), comparison_);
360  nodes_.pop_back();
361}
362// Test if empty - does work if so
363bool 
364CbcTreeLocal::empty()   
365{
366  if (typeCuts_<0) 
367    return !nodes_.size();
368  /* state -
369     0 iterating
370     1 subtree finished optimal solution for subtree found
371     2 subtree finished and no solution found
372     3 subtree exiting and solution found
373     4 subtree exiting and no solution found
374  */
375  int state=0;
376  assert (searchType_!=2);
377  if (searchType_) {
378    if (CoinCpuTime()-startTime_>timeLimit_||model_->getNodeCount()-startNode_>=nodeLimit_) {
379      state=4;
380    }
381  } else {
382    if (model_->getSolutionCount()>saveNumberSolutions_) {
383      state=4;
384    }
385  }
386  if (!nodes_.size()) 
387    state=2;
388  if (!state) {
389    return false;
390  }
391  // Finished this phase
392  int numberColumns = model_->getNumCols();
393  if (model_->getSolutionCount()>saveNumberSolutions_) {
394    if(model_->getCutoff()<bestCutoff_) {
395      // Save solution
396      if (!bestSolution_) 
397        bestSolution_ = new double [numberColumns];
398      memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
399      bestCutoff_=model_->getCutoff();
400    }
401    state--;
402  } 
403  // get rid of all nodes (safe even if already done)
404  double bestPossibleObjective;
405  cleanTree(model_,-COIN_DBL_MAX,bestPossibleObjective);
406
407  double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
408  printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
409         state,
410         model_->getNodeCount()-startNode_,
411         (int) CoinCpuTime()-startTime_,
412         model_->getCutoff()+increment,bestCutoff_+increment,rhs_);
413  saveNumberSolutions_ = model_->getSolutionCount();
414  bool finished=false;
415  bool lastTry=false;
416  switch (state) {
417  case 1:
418    // solution found and subtree exhausted
419    if (rhs_>1.0e30) {
420      finished=true;
421    } else {
422      // find global cut and reverse
423      reverseCut(1);
424      searchType_=1; // first false
425      rhs_ = range_; // reset range
426      nextStrong_=false;
427
428      // save best solution in this subtree
429      memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
430    }
431    break;
432  case 2:
433    // solution not found and subtree exhausted
434    if (rhs_>1.0e30) {
435      finished=true;
436    } else {
437      // find global cut and reverse
438      reverseCut(2);
439      searchType_=1; // first false
440      if (diversification_<maxDiversification_) {
441        if (nextStrong_) {
442          diversification_++;
443          model_->setCutoff(1.0e50);
444          searchType_=0;
445        }
446        nextStrong_=true;
447        rhs_ += range_/2;
448      } else {
449        // This will be last try (may hit max time0
450        lastTry=true;
451        model_->setCutoff(bestCutoff_);
452        printf("Exiting local search with current set of cuts\n");
453        rhs_=1.0e100;
454        // Can now stop on gap
455        model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
456      }
457    }
458    break;
459  case 3:
460    // solution found and subtree not exhausted
461    if (rhs_<1.0e30) {
462      if (searchType_) {
463        if (!typeCuts_&&refine_&&searchType_==1) {
464          // We need to check we have best solution given these 0-1 values
465          OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
466          CbcModel * subModel = model_->subTreeModel(subSolver);
467          CbcTree normalTree;
468          subModel->passInTreeHandler(normalTree);
469          int numberIntegers = model_->numberIntegers();
470          const int * integerVariable = model_->integerVariable();
471          const double * solution = model_->bestSolution();
472          int i;
473          int numberColumns = model_->getNumCols();
474          for (i=0;i<numberIntegers;i++) {
475            int iColumn = integerVariable[i];
476            double value = floor(solution[iColumn]+0.5);
477            if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
478              continue; // skip as not 0-1
479            if (originalLower_[i]==originalUpper_[i])
480              continue;
481            subSolver->setColLower(iColumn,value);
482            subSolver->setColUpper(iColumn,value);
483          }
484          subSolver->initialSolve();
485          // We can copy cutoff
486          subModel->setCutoff(model_->getCutoff());
487          subModel->setSolutionCount(0);
488          assert (subModel->isProvenOptimal());
489          if (!subModel->typePresolve()) {
490            subModel->branchAndBound();
491            if (subModel->status()) {
492              model_->incrementSubTreeStopped(); 
493            }
494            printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
495                   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
496            double newCutoff = subModel->getMinimizationObjValue()-
497              subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
498            if (subModel->getSolutionCount()) {
499              if (!subModel->status())
500                assert (subModel->isProvenOptimal());
501              memcpy(model_->bestSolution(),subModel->bestSolution(),
502                     numberColumns*sizeof(double));
503              model_->setCutoff(newCutoff);
504            }
505          } else if (subModel->typePresolve()==1) {
506            CbcModel * model2 = subModel->integerPresolve(true);
507            if (model2) {
508              // Do complete search
509              model2->branchAndBound();
510              // get back solution
511              subModel->originalModel(model2,false);
512              if (model2->status()) {
513                model_->incrementSubTreeStopped(); 
514              } 
515              double newCutoff = model2->getMinimizationObjValue()-
516                model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
517              if (model2->getSolutionCount()) {
518                if (!model2->status())
519                  assert (model2->isProvenOptimal());
520                memcpy(model_->bestSolution(),subModel->bestSolution(),
521                       numberColumns*sizeof(double));
522                model_->setCutoff(newCutoff);
523              }
524              delete model2;
525            } else {
526              // infeasible - could just be - due to cutoff
527            }
528          } else {
529            // too dangerous at present
530            assert (subModel->typePresolve()!=2);
531          }
532          if(model_->getCutoff()<bestCutoff_) {
533            // Save solution
534            if (!bestSolution_) 
535              bestSolution_ = new double [numberColumns];
536            memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
537            bestCutoff_=model_->getCutoff();
538          }
539          delete subModel;
540        }
541        // we have done search to make sure best general solution
542        searchType_=1;
543        // Reverse cut weakly
544        reverseCut(3,rhs_);
545      } else {
546        searchType_=1;
547        // delete last cut
548        deleteCut(cut_);
549      }
550    } else {
551      searchType_=1;
552    }
553    // save best solution in this subtree
554    memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
555    nextStrong_=false;
556    rhs_=range_;
557    break;
558  case 4:
559    // solution not found and subtree not exhausted
560    if (nextStrong_) {
561      // Reverse cut weakly
562      reverseCut(4,rhs_);
563      model_->setCutoff(1.0e50);
564      diversification_++;
565      searchType_=0;
566    } else {
567      // delete last cut
568      deleteCut(cut_);
569      searchType_=1;
570    }
571    nextStrong_=true;
572    rhs_ += range_/2;
573    break;
574  }
575  if (rhs_<1.0e30||lastTry) {
576    int goodSolution=createCut(savedSolution_,cut_);
577    assert(goodSolution>=0);
578    // Add to global cuts
579    model_->globalCuts()->insert(cut_);
580    {
581      OsiCuts * global = model_->globalCuts();
582      int n = global->sizeRowCuts();
583      OsiRowCut * rowCut = global->rowCutPtr(n-1);
584      printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
585             n,rowCut->lb(),rowCut->ub(),rhs_,diversification_);
586      for (int i=0;i<n;i++) {
587        rowCut = global->rowCutPtr(i);
588        printf("%d - rhs %g %g\n",
589               i,rowCut->lb(),rowCut->ub());
590      }
591    }
592    // put back node
593    startTime_ = (int) CoinCpuTime();
594    startNode_=model_->getNodeCount();
595    // save copy of node
596    CbcNode * localNode2 = new CbcNode(*localNode_);
597    // But localNode2 now owns cuts so swap
598    //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
599    //   localNode_,localNode_->nodeInfo());
600    nodes_.push_back(localNode_);
601    localNode_=localNode2;
602    make_heap(nodes_.begin(), nodes_.end(), comparison_);
603  }
604  return finished;
605}
606// We may have got an intelligent tree so give it one more chance
607void 
608CbcTreeLocal::endSearch() 
609{
610  if (typeCuts_>=0) {
611    // copy best solution to model
612    int numberColumns = model_->getNumCols();
613    if (bestSolution_&&bestCutoff_<model_->getCutoff()) {
614      memcpy(model_->bestSolution(),bestSolution_,numberColumns*sizeof(double));
615      model_->setCutoff(bestCutoff_);
616      // recompute objective value
617      const double * objCoef = model_->getObjCoefficients();
618      double objOffset=0.0;
619      model_->continuousSolver()->getDblParam(OsiObjOffset,objOffset);
620 
621      // Compute dot product of objCoef and colSol and then adjust by offset
622      double objValue = -objOffset;
623      for ( int i=0 ; i<numberColumns ; i++ )
624        objValue += objCoef[i]*bestSolution_[i];
625      model_->setMinimizationObjValue(objValue);
626    }
627    // Can now stop on gap
628    model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
629  }
630}
631// Create cut
632int
633CbcTreeLocal::createCut(const double * solution,OsiRowCut & rowCut)
634{
635  OsiSolverInterface * solver = model_->solver();
636  const double * rowLower = solver->getRowLower();
637  const double * rowUpper = solver->getRowUpper();
638  //const double * solution = solver->getColSolution();
639  //const double * objective = solver->getObjCoefficients();
640  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
641  double primalTolerance;
642  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
643
644  int numberRows = model_->getNumRows();
645
646  int numberIntegers = model_->numberIntegers();
647  const int * integerVariable = model_->integerVariable();
648  int i;
649
650  // Check feasible
651
652  double * rowActivity = new double[numberRows];
653  memset(rowActivity,0,numberRows*sizeof(double));
654  solver->getMatrixByCol()->times(solution,rowActivity) ;
655  int goodSolution=0;
656  // check was feasible
657  for (i=0;i<numberRows;i++) {
658    if(rowActivity[i]<rowLower[i]-primalTolerance) {
659      goodSolution=-1;
660    } else if(rowActivity[i]>rowUpper[i]+primalTolerance) {
661      goodSolution=-1;
662    }
663  }
664  delete [] rowActivity;
665  for (i=0;i<numberIntegers;i++) {
666    int iColumn = integerVariable[i];
667    double value=solution[iColumn];
668    if (fabs(floor(value+0.5)-value)>integerTolerance) {
669      goodSolution=-1;
670    }
671  }
672  // zap cut
673  if (goodSolution==0) {
674    // Create cut and get total gap
675    CoinPackedVector cut;
676    double rhs = rhs_;
677    double maxValue=0.0;
678    for (i=0;i<numberIntegers;i++) {
679      int iColumn = integerVariable[i];
680      double value = floor(solution[iColumn]+0.5);
681      if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
682        continue; // skip as not 0-1
683      if (originalLower_[i]==originalUpper_[i])
684        continue;
685      double mu =1.0/(originalUpper_[i]-originalLower_[i]);
686      if (value==originalLower_[i]) {
687        rhs += mu*originalLower_[i];
688        cut.insert(iColumn,1.0);
689        maxValue += originalUpper_[i];
690      } else if (value==originalUpper_[i]) {
691        rhs -= mu*originalUpper_[i];
692        cut.insert(iColumn,-1.0);
693        maxValue += originalLower_[i];
694      }
695    }
696    if (maxValue<rhs-primalTolerance) {
697      printf("slack cut\n");
698      goodSolution=1;
699    }
700    rowCut.setRow(cut);
701    rowCut.setLb(-COIN_DBL_MAX);
702    rowCut.setUb(rhs);   
703    rowCut.setGloballyValid();
704    printf("Cut size: %i Cut rhs: %g\n",cut.getNumElements(), rhs);
705#ifdef CBC_DEBUG
706    int k;
707    for (k=0; k<cut.getNumElements(); k++){
708      printf("%i    %g ",cut.getIndices()[k], cut.getElements()[k]);
709      if ((k+1)%5==0)
710        printf("\n");
711    }
712    if (k%5!=0)
713      printf("\n");
714#endif
715    return goodSolution;
716  } else {
717    printf("Not a good solution\n");
718    return -1;
719  }
720}
721// Other side of last cut branch
722void 
723CbcTreeLocal::reverseCut(int state, double bias)
724{
725  // find global cut
726  OsiCuts * global = model_->globalCuts();
727  int n = global->sizeRowCuts();
728  int i;
729  OsiRowCut * rowCut = NULL;
730  for ( i=0;i<n;i++) {
731    rowCut = global->rowCutPtr(i);
732    if (cut_==*rowCut) {
733      break;
734    }
735  }
736  assert (i<n);
737  // get smallest element
738  double smallest=COIN_DBL_MAX;
739  CoinPackedVector row = cut_.row();
740  for (int k=0; k<row.getNumElements(); k++)
741    smallest = min(smallest,fabs(row.getElements()[k]));
742  if (!typeCuts_&&!refine_) {
743    // Reverse cut very very weakly
744    if (state>2)
745      smallest=0.0;
746  }
747  // replace by other way
748  printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
749         i,n,rowCut->lb(),rowCut->ub());
750  rowCut->setLb(rowCut->ub()+smallest-bias);
751  rowCut->setUb(COIN_DBL_MAX);
752  printf("new rhs %g %g, bias %g smallest %g ",
753         rowCut->lb(),rowCut->ub(),bias,smallest);
754}
755// Delete last cut branch
756void 
757CbcTreeLocal::deleteCut(OsiRowCut & cut)
758{
759  // find global cut
760  OsiCuts * global = model_->globalCuts();
761  int n = global->sizeRowCuts();
762  int i;
763  OsiRowCut * rowCut = NULL;
764  for ( i=0;i<n;i++) {
765    rowCut = global->rowCutPtr(i);
766    if (cut==*rowCut) {
767      break;
768    }
769  }
770  assert (i<n);
771  // delete last cut
772  printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
773         i,n,rowCut->lb(),rowCut->ub());
774  global->eraseRowCut(i);
775}
776
777
778
Note: See TracBrowser for help on using the repository browser.