source: branches/devel/Cbc/src/CbcTreeLocal.cpp @ 419

Last change on this file since 419 was 419, checked in by forrest, 13 years ago

possible bug in tree

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