source: trunk/Cbc/src/CbcTreeLocal.cpp @ 862

Last change on this file since 862 was 765, checked in by andreasw, 12 years ago

merging changes from Bug Squashing Party Aug 2007 to regular trunk

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