source: stable/2.4/Cbc/src/CbcTreeLocal.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 53.4 KB
Line 
1/* $Id: CbcTreeLocal.cpp 1271 2009-11-05 15:57:25Z forrest $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CbcModel.hpp"
6#include "CbcNode.hpp"
7#include "CbcTreeLocal.hpp"
8#include "CoinPackedMatrix.hpp"
9#include "CoinTime.hpp"
10#include "OsiRowCutDebugger.hpp"
11#include <cassert>
12#if 0
13// gdb doesn't always put breakpoints in this virtual function
14// just stick xxxxxx() where you want to start
15static void xxxxxx()
16{
17  printf("break\n");
18}
19#endif
20CbcTreeLocal::CbcTreeLocal()
21  : localNode_(NULL),
22    bestSolution_(NULL),
23    savedSolution_(NULL),
24    saveNumberSolutions_(0),
25    model_(NULL),
26    originalLower_(NULL),
27    originalUpper_(NULL),
28    range_(0),
29    typeCuts_(-1),
30    maxDiversification_(0),
31    diversification_(0),
32    nextStrong_(false),
33    rhs_(0.0),
34    savedGap_(0.0),
35    bestCutoff_(0.0),
36    timeLimit_(0),
37    startTime_(0),
38    nodeLimit_(0),
39    startNode_(-1),
40    searchType_(-1),
41    refine_(false)
42{
43 
44}
45/* Constructor with solution.
46   range is upper bound on difference from given solution.
47   maxDiversification is maximum number of diversifications to try
48   timeLimit is seconds in subTree
49   nodeLimit is nodes in subTree
50*/
51CbcTreeLocal::CbcTreeLocal(CbcModel * model,const double * solution ,
52                                  int range,int typeCuts,int maxDiversification,
53                                  int timeLimit, int nodeLimit,bool refine)
54  : localNode_(NULL),
55    bestSolution_(NULL),
56    savedSolution_(NULL),
57    saveNumberSolutions_(0),
58    model_(model),
59    originalLower_(NULL),
60    originalUpper_(NULL),
61    range_(range),
62    typeCuts_(typeCuts),
63    maxDiversification_(maxDiversification),
64    diversification_(0),
65    nextStrong_(false),
66    rhs_(0.0),
67    savedGap_(0.0),
68    bestCutoff_(0.0),
69    timeLimit_(timeLimit),
70    startTime_(0),
71    nodeLimit_(nodeLimit),
72    startNode_(-1),
73    searchType_(-1),
74    refine_(refine)
75{
76
77  OsiSolverInterface * solver = model_->solver();
78  const double * lower = solver->getColLower();
79  const double * upper = solver->getColUpper();
80  //const double * solution = solver->getColSolution();
81  //const double * objective = solver->getObjCoefficients();
82  double primalTolerance;
83  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
84
85  // Get increment
86  model_->analyzeObjective();
87
88  { 
89    // needed to sync cutoffs
90    double value ;
91    solver->getDblParam(OsiDualObjectiveLimit,value) ;
92    model_->setCutoff(value * solver->getObjSense());
93  }
94  bestCutoff_ = model_->getCutoff();
95  // save current gap
96  savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
97
98  // make sure integers found
99  model_->findIntegers(false);
100  int numberIntegers = model_->numberIntegers();
101  const int * integerVariable = model_->integerVariable();
102  int i;
103  double direction = solver->getObjSense();
104  double newSolutionValue=1.0e50;
105  if (solution) {
106    // copy solution
107    solver->setColSolution(solution);
108    newSolutionValue = direction*solver->getObjValue();
109  }
110  originalLower_=new double [numberIntegers];
111  originalUpper_=new double [numberIntegers];
112  bool all01=true;
113  int number01=0;
114  for (i=0;i<numberIntegers;i++) {
115    int iColumn = integerVariable[i];
116    originalLower_[i]=lower[iColumn];
117    originalUpper_[i]=upper[iColumn];
118    if (upper[iColumn]-lower[iColumn]>1.5)
119      all01=false;
120    else if (upper[iColumn]-lower[iColumn]==1.0)
121      number01++;
122  }
123  if (all01 && !typeCuts_)
124    typeCuts_=1; // may as well so we don't have to deal with refine
125  if (!number01&&!typeCuts_) {
126    if (model_->messageHandler()->logLevel()>0)
127      printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
128    typeCuts_=-1;
129  } else {
130    if (model_->messageHandler()->logLevel()>0) {
131      std::string type;
132      if (all01) {
133        printf("%d 0-1 variables normal local  cuts\n",
134               number01);
135      } else if (typeCuts_) {
136        printf("%d 0-1 variables, %d other - general integer local cuts\n",
137               number01,numberIntegers-number01);
138      } else {
139        printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
140               number01,numberIntegers-number01);
141      }
142      printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
143             maxDiversification_,range_,timeLimit_,nodeLimit_);
144    }
145  }
146  int numberColumns = model_->getNumCols();
147  savedSolution_ = new double [numberColumns];
148  memset(savedSolution_,0,numberColumns*sizeof(double));
149  if (solution) {
150    rhs_=range_;
151    // Check feasible
152    int goodSolution=createCut(solution,cut_);
153    if (goodSolution>=0) {
154      for (i=0;i<numberIntegers;i++) {
155        int iColumn = integerVariable[i];
156        double value = floor(solution[iColumn]+0.5);
157        // fix so setBestSolution will work
158        solver->setColLower(iColumn,value);
159        solver->setColUpper(iColumn,value);
160      }
161      model_->reserveCurrentSolution();
162      // Create cut and get total gap
163      if (newSolutionValue<bestCutoff_) {
164        model_->setBestSolution(CBC_ROUNDING,newSolutionValue,solution);
165        bestCutoff_=model_->getCutoff();
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_=CoinMin(solutionValue,model_->getCutoff());
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=CoinMax(largest,nn);
344    smallest=CoinMin(smallest,nn);
345    largestD=CoinMax(largestD,dd);
346    smallestD=CoinMin(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_ = static_cast<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           static_cast<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          // cut is valid so don't 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_ = static_cast<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 = CoinMin(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
882CbcTreeVariable::CbcTreeVariable()
883  : localNode_(NULL),
884    bestSolution_(NULL),
885    savedSolution_(NULL),
886    saveNumberSolutions_(0),
887    model_(NULL),
888    originalLower_(NULL),
889    originalUpper_(NULL),
890    range_(0),
891    typeCuts_(-1),
892    maxDiversification_(0),
893    diversification_(0),
894    nextStrong_(false),
895    rhs_(0.0),
896    savedGap_(0.0),
897    bestCutoff_(0.0),
898    timeLimit_(0),
899    startTime_(0),
900    nodeLimit_(0),
901    startNode_(-1),
902    searchType_(-1),
903    refine_(false)
904{
905 
906}
907/* Constructor with solution.
908   range is upper bound on difference from given solution.
909   maxDiversification is maximum number of diversifications to try
910   timeLimit is seconds in subTree
911   nodeLimit is nodes in subTree
912*/
913CbcTreeVariable::CbcTreeVariable(CbcModel * model,const double * solution ,
914                                  int range,int typeCuts,int maxDiversification,
915                                  int timeLimit, int nodeLimit,bool refine)
916  : localNode_(NULL),
917    bestSolution_(NULL),
918    savedSolution_(NULL),
919    saveNumberSolutions_(0),
920    model_(model),
921    originalLower_(NULL),
922    originalUpper_(NULL),
923    range_(range),
924    typeCuts_(typeCuts),
925    maxDiversification_(maxDiversification),
926    diversification_(0),
927    nextStrong_(false),
928    rhs_(0.0),
929    savedGap_(0.0),
930    bestCutoff_(0.0),
931    timeLimit_(timeLimit),
932    startTime_(0),
933    nodeLimit_(nodeLimit),
934    startNode_(-1),
935    searchType_(-1),
936    refine_(refine)
937{
938
939  OsiSolverInterface * solver = model_->solver();
940  const double * lower = solver->getColLower();
941  const double * upper = solver->getColUpper();
942  //const double * solution = solver->getColSolution();
943  //const double * objective = solver->getObjCoefficients();
944  double primalTolerance;
945  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
946
947  // Get increment
948  model_->analyzeObjective();
949
950  { 
951    // needed to sync cutoffs
952    double value ;
953    solver->getDblParam(OsiDualObjectiveLimit,value) ;
954    model_->setCutoff(value * solver->getObjSense());
955  }
956  bestCutoff_ = model_->getCutoff();
957  // save current gap
958  savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
959
960  // make sure integers found
961  model_->findIntegers(false);
962  int numberIntegers = model_->numberIntegers();
963  const int * integerVariable = model_->integerVariable();
964  int i;
965  double direction = solver->getObjSense();
966  double newSolutionValue=1.0e50;
967  if (solution) {
968    // copy solution
969    solver->setColSolution(solution);
970    newSolutionValue = direction*solver->getObjValue();
971  }
972  originalLower_=new double [numberIntegers];
973  originalUpper_=new double [numberIntegers];
974  bool all01=true;
975  int number01=0;
976  for (i=0;i<numberIntegers;i++) {
977    int iColumn = integerVariable[i];
978    originalLower_[i]=lower[iColumn];
979    originalUpper_[i]=upper[iColumn];
980    if (upper[iColumn]-lower[iColumn]>1.5)
981      all01=false;
982    else if (upper[iColumn]-lower[iColumn]==1.0)
983      number01++;
984  }
985  if (all01 && !typeCuts_)
986    typeCuts_=1; // may as well so we don't have to deal with refine
987  if (!number01&&!typeCuts_) {
988    if (model_->messageHandler()->logLevel()>0)
989      printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
990    typeCuts_=-1;
991  } else {
992    if (model_->messageHandler()->logLevel()>0) {
993      std::string type;
994      if (all01) {
995        printf("%d 0-1 variables normal local  cuts\n",
996               number01);
997      } else if (typeCuts_) {
998        printf("%d 0-1 variables, %d other - general integer local cuts\n",
999               number01,numberIntegers-number01);
1000      } else {
1001        printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
1002               number01,numberIntegers-number01);
1003      }
1004      printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
1005             maxDiversification_,range_,timeLimit_,nodeLimit_);
1006    }
1007  }
1008  int numberColumns = model_->getNumCols();
1009  savedSolution_ = new double [numberColumns];
1010  memset(savedSolution_,0,numberColumns*sizeof(double));
1011  if (solution) {
1012    rhs_=range_;
1013    // Check feasible
1014    int goodSolution=createCut(solution,cut_);
1015    if (goodSolution>=0) {
1016      for (i=0;i<numberIntegers;i++) {
1017        int iColumn = integerVariable[i];
1018        double value = floor(solution[iColumn]+0.5);
1019        // fix so setBestSolution will work
1020        solver->setColLower(iColumn,value);
1021        solver->setColUpper(iColumn,value);
1022      }
1023      model_->reserveCurrentSolution();
1024      // Create cut and get total gap
1025      if (newSolutionValue<bestCutoff_) {
1026        model_->setBestSolution(CBC_ROUNDING,newSolutionValue,solution);
1027        bestCutoff_=model_->getCutoff();
1028        // save as best solution
1029        memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
1030      }
1031      for (i=0;i<numberIntegers;i++) {
1032        int iColumn = integerVariable[i];
1033        // restore bounds
1034        solver->setColLower(iColumn,originalLower_[i]);
1035        solver->setColUpper(iColumn,originalUpper_[i]);
1036      }
1037      // make sure can't stop on gap
1038      model_->setDblParam(CbcModel::CbcAllowableGap,-1.0e50);
1039    } else {
1040      model_=NULL;
1041    }
1042  } else {
1043    // no solution
1044    rhs_=1.0e50;
1045    // make sure can't stop on gap
1046    model_->setDblParam(CbcModel::CbcAllowableGap,-1.0e50);
1047  }
1048}
1049CbcTreeVariable::~CbcTreeVariable()
1050{
1051  delete [] originalLower_;
1052  delete [] originalUpper_;
1053  delete [] bestSolution_;
1054  delete [] savedSolution_;
1055  delete localNode_;
1056}
1057// Copy constructor
1058CbcTreeVariable::CbcTreeVariable ( const CbcTreeVariable & rhs)
1059  :CbcTree(rhs),
1060   saveNumberSolutions_(rhs.saveNumberSolutions_),
1061   model_(rhs.model_),
1062   range_(rhs.range_),
1063   typeCuts_(rhs.typeCuts_),
1064   maxDiversification_(rhs.maxDiversification_),
1065   diversification_(rhs.diversification_),
1066   nextStrong_(rhs.nextStrong_),
1067   rhs_(rhs.rhs_),
1068   savedGap_(rhs.savedGap_),
1069   bestCutoff_(rhs.bestCutoff_),
1070   timeLimit_(rhs.timeLimit_),
1071   startTime_(rhs.startTime_),
1072   nodeLimit_(rhs.nodeLimit_),
1073   startNode_(rhs.startNode_),
1074   searchType_(rhs.searchType_),
1075   refine_(rhs.refine_)
1076{ 
1077  cut_=rhs.cut_;
1078  fixedCut_=rhs.fixedCut_;
1079  if (rhs.localNode_)
1080    localNode_=new CbcNode(*rhs.localNode_);
1081  else
1082    localNode_=NULL;
1083  if (rhs.originalLower_) {
1084    int numberIntegers = model_->numberIntegers();
1085    originalLower_=new double [numberIntegers];
1086    memcpy(originalLower_,rhs.originalLower_,numberIntegers*sizeof(double));
1087    originalUpper_=new double [numberIntegers];
1088    memcpy(originalUpper_,rhs.originalUpper_,numberIntegers*sizeof(double));
1089  } else {
1090    originalLower_=NULL;
1091    originalUpper_=NULL;
1092  }
1093  if (rhs.bestSolution_) {
1094    int numberColumns = model_->getNumCols();
1095    bestSolution_ = new double [numberColumns];
1096    memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
1097  } else {
1098    bestSolution_=NULL;
1099  }
1100  if (rhs.savedSolution_) {
1101    int numberColumns = model_->getNumCols();
1102    savedSolution_ = new double [numberColumns];
1103    memcpy(savedSolution_,rhs.savedSolution_,numberColumns*sizeof(double));
1104  } else {
1105    savedSolution_=NULL;
1106  }
1107}
1108//----------------------------------------------------------------
1109// Assignment operator
1110//-------------------------------------------------------------------
1111CbcTreeVariable &
1112CbcTreeVariable::operator=(const CbcTreeVariable& rhs)
1113{
1114  if (this != &rhs) {
1115    CbcTree::operator=(rhs);
1116    saveNumberSolutions_ = rhs.saveNumberSolutions_;
1117    cut_=rhs.cut_;
1118    fixedCut_=rhs.fixedCut_;
1119    delete localNode_;
1120    if (rhs.localNode_)
1121      localNode_=new CbcNode(*rhs.localNode_);
1122    else
1123      localNode_=NULL;
1124    model_ = rhs.model_;
1125    range_ = rhs.range_;
1126    typeCuts_ = rhs.typeCuts_;
1127    maxDiversification_ = rhs.maxDiversification_;
1128    diversification_ = rhs.diversification_;
1129    nextStrong_ = rhs.nextStrong_;
1130    rhs_ = rhs.rhs_;
1131    savedGap_ = rhs.savedGap_;
1132    bestCutoff_ = rhs.bestCutoff_;
1133    timeLimit_ = rhs.timeLimit_;
1134    startTime_ = rhs.startTime_;
1135    nodeLimit_ = rhs.nodeLimit_;
1136    startNode_ = rhs.startNode_;
1137    searchType_ = rhs.searchType_;
1138    refine_ = rhs.refine_;
1139    delete [] originalLower_;
1140    delete [] originalUpper_;
1141    if (rhs.originalLower_) {
1142      int numberIntegers = model_->numberIntegers();
1143      originalLower_=new double [numberIntegers];
1144      memcpy(originalLower_,rhs.originalLower_,numberIntegers*sizeof(double));
1145      originalUpper_=new double [numberIntegers];
1146      memcpy(originalUpper_,rhs.originalUpper_,numberIntegers*sizeof(double));
1147    } else {
1148      originalLower_=NULL;
1149      originalUpper_=NULL;
1150    }
1151    delete [] bestSolution_;
1152    if (rhs.bestSolution_) {
1153      int numberColumns = model_->getNumCols();
1154      bestSolution_ = new double [numberColumns];
1155      memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
1156    } else {
1157      bestSolution_=NULL;
1158    }
1159    delete [] savedSolution_;
1160    if (rhs.savedSolution_) {
1161      int numberColumns = model_->getNumCols();
1162      savedSolution_ = new double [numberColumns];
1163      memcpy(savedSolution_,rhs.savedSolution_,numberColumns*sizeof(double));
1164    } else {
1165      savedSolution_=NULL;
1166    }
1167  }
1168  return *this;
1169}
1170// Clone
1171CbcTree *
1172CbcTreeVariable::clone() const
1173{
1174  return new CbcTreeVariable(*this);
1175}
1176// Pass in solution (so can be used after heuristic)
1177void 
1178CbcTreeVariable::passInSolution(const double * solution, double solutionValue)
1179{
1180  int numberColumns = model_->getNumCols();
1181  delete [] savedSolution_;
1182  savedSolution_ = new double [numberColumns];
1183  memcpy(savedSolution_,solution,numberColumns*sizeof(double));
1184  rhs_=range_;
1185  // Check feasible
1186  int goodSolution=createCut(solution,cut_);
1187  if (goodSolution>=0) {
1188    bestCutoff_=CoinMin(solutionValue,model_->getCutoff());
1189  } else {
1190    model_=NULL;
1191  }
1192}
1193// Return the top node of the heap
1194CbcNode * 
1195CbcTreeVariable::top() const{
1196#ifdef CBC_DEBUG
1197  int smallest=9999999;
1198  int largest=-1;
1199  double smallestD=1.0e30;
1200  double largestD=-1.0e30;
1201  int n=nodes_.size();
1202  for (int i=0;i<n;i++) {
1203    int nn=nodes_[i]->nodeInfo()->nodeNumber();
1204    double dd = nodes_[i]->objectiveValue();
1205    largest=CoinMax(largest,nn);
1206    smallest=CoinMin(smallest,nn);
1207    largestD=CoinMax(largestD,dd);
1208    smallestD=CoinMin(smallestD,dd);
1209  }
1210  if (model_->messageHandler()->logLevel()>0) {
1211    printf("smallest %d, largest %d, top %d\n",smallest,largest,
1212           nodes_.front()->nodeInfo()->nodeNumber());
1213    printf("smallestD %g, largestD %g, top %g\n",smallestD,largestD,nodes_.front()->objectiveValue());
1214  }
1215#endif
1216  return nodes_.front();
1217}
1218
1219// Add a node to the heap
1220void 
1221CbcTreeVariable::push(CbcNode * x) {
1222  if (typeCuts_>=0&&!nodes_.size()&&searchType_<0) {
1223    startNode_=model_->getNodeCount();
1224    // save copy of node
1225    localNode_ = new CbcNode(*x);
1226
1227    if (cut_.row().getNumElements()) {
1228      // Add to global cuts
1229      // we came in with solution
1230      model_->globalCuts()->insert(cut_);
1231      if (model_->messageHandler()->logLevel()>0)
1232        printf("initial cut - rhs %g %g\n",
1233               cut_.lb(),cut_.ub());
1234      searchType_=1;
1235    } else {
1236      // stop on first solution
1237      searchType_=0;
1238    }
1239    startTime_ = static_cast<int> (CoinCpuTime());
1240    saveNumberSolutions_ = model_->getSolutionCount();
1241  }
1242  nodes_.push_back(x);
1243#ifdef CBC_DEBUG
1244  if (model_->messageHandler()->logLevel()>0)
1245    printf("pushing node onto heap %d %x %x\n",
1246           x->nodeInfo()->nodeNumber(),x,x->nodeInfo());
1247#endif
1248  std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
1249}
1250
1251// Remove the top node from the heap
1252void 
1253CbcTreeVariable::pop() {
1254  std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
1255  nodes_.pop_back();
1256}
1257// Test if empty - does work if so
1258bool 
1259CbcTreeVariable::empty() 
1260{
1261  if (typeCuts_<0) 
1262    return !nodes_.size();
1263  /* state -
1264     0 iterating
1265     1 subtree finished optimal solution for subtree found
1266     2 subtree finished and no solution found
1267     3 subtree exiting and solution found
1268     4 subtree exiting and no solution found
1269  */
1270  int state=0;
1271  assert (searchType_!=2);
1272  if (searchType_) {
1273    if (CoinCpuTime()-startTime_>timeLimit_||model_->getNodeCount()-startNode_>=nodeLimit_) {
1274      state=4;
1275    }
1276  } else {
1277    if (model_->getSolutionCount()>saveNumberSolutions_) {
1278      state=4;
1279    }
1280  }
1281  if (!nodes_.size()) 
1282    state=2;
1283  if (!state) {
1284    return false;
1285  }
1286  // Finished this phase
1287  int numberColumns = model_->getNumCols();
1288  if (model_->getSolutionCount()>saveNumberSolutions_) {
1289    if(model_->getCutoff()<bestCutoff_) {
1290      // Save solution
1291      if (!bestSolution_) 
1292        bestSolution_ = new double [numberColumns];
1293      memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
1294      bestCutoff_=model_->getCutoff();
1295    }
1296    state--;
1297  } 
1298  // get rid of all nodes (safe even if already done)
1299  double bestPossibleObjective;
1300  cleanTree(model_,-COIN_DBL_MAX,bestPossibleObjective);
1301
1302  double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
1303  if (model_->messageHandler()->logLevel()>0)
1304    printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
1305           state,
1306           model_->getNodeCount()-startNode_,
1307           static_cast<int> (CoinCpuTime())-startTime_,
1308           model_->getCutoff()+increment,bestCutoff_+increment,rhs_);
1309  saveNumberSolutions_ = model_->getSolutionCount();
1310  bool finished=false;
1311  bool lastTry=false;
1312  switch (state) {
1313  case 1:
1314    // solution found and subtree exhausted
1315    if (rhs_>1.0e30) {
1316      finished=true;
1317    } else {
1318      // find global cut and reverse
1319      reverseCut(1);
1320      searchType_=1; // first false
1321      rhs_ = range_; // reset range
1322      nextStrong_=false;
1323
1324      // save best solution in this subtree
1325      memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
1326    }
1327    break;
1328  case 2:
1329    // solution not found and subtree exhausted
1330    if (rhs_>1.0e30) {
1331      finished=true;
1332    } else {
1333      // find global cut and reverse
1334      reverseCut(2);
1335      searchType_=1; // first false
1336      if (diversification_<maxDiversification_) {
1337        if (nextStrong_) {
1338          diversification_++;
1339          // cut is valid so don't model_->setCutoff(1.0e50);
1340          searchType_=0;
1341        }
1342        nextStrong_=true;
1343        rhs_ += range_/2;
1344      } else {
1345        // This will be last try (may hit max time)
1346        lastTry=true;
1347        if (!maxDiversification_)
1348          typeCuts_=-1; // make sure can't start again
1349        model_->setCutoff(bestCutoff_);
1350        if (model_->messageHandler()->logLevel()>0)
1351          printf("Exiting local search with current set of cuts\n");
1352        rhs_=1.0e100;
1353        // Can now stop on gap
1354        model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
1355      }
1356    }
1357    break;
1358  case 3:
1359    // solution found and subtree not exhausted
1360    if (rhs_<1.0e30) {
1361      if (searchType_) {
1362        if (!typeCuts_&&refine_&&searchType_==1) {
1363          // We need to check we have best solution given these 0-1 values
1364          OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
1365          CbcModel * subModel = model_->subTreeModel(subSolver);
1366          CbcTree normalTree;
1367          subModel->passInTreeHandler(normalTree);
1368          int numberIntegers = model_->numberIntegers();
1369          const int * integerVariable = model_->integerVariable();
1370          const double * solution = model_->bestSolution();
1371          int i;
1372          int numberColumns = model_->getNumCols();
1373          for (i=0;i<numberIntegers;i++) {
1374            int iColumn = integerVariable[i];
1375            double value = floor(solution[iColumn]+0.5);
1376            if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
1377              continue; // skip as not 0-1
1378            if (originalLower_[i]==originalUpper_[i])
1379              continue;
1380            subSolver->setColLower(iColumn,value);
1381            subSolver->setColUpper(iColumn,value);
1382          }
1383          subSolver->initialSolve();
1384          // We can copy cutoff
1385          // But adjust
1386          subModel->setCutoff(model_->getCutoff()+model_->getDblParam(CbcModel::CbcCutoffIncrement)+1.0e-6);
1387          subModel->setSolutionCount(0);
1388          assert (subModel->isProvenOptimal());
1389          if (!subModel->typePresolve()) {
1390            subModel->branchAndBound();
1391            if (subModel->status()) {
1392              model_->incrementSubTreeStopped(); 
1393            }
1394            //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
1395            //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
1396            double newCutoff = subModel->getMinimizationObjValue()-
1397              subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
1398            if (subModel->getSolutionCount()) {
1399              if (!subModel->status())
1400                assert (subModel->isProvenOptimal());
1401              memcpy(model_->bestSolution(),subModel->bestSolution(),
1402                     numberColumns*sizeof(double));
1403              model_->setCutoff(newCutoff);
1404            }
1405          } else if (subModel->typePresolve()==1) {
1406            CbcModel * model2 = subModel->integerPresolve(true);
1407            if (model2) {
1408              // Do complete search
1409              model2->branchAndBound();
1410              // get back solution
1411              subModel->originalModel(model2,false);
1412              if (model2->status()) {
1413                model_->incrementSubTreeStopped(); 
1414              } 
1415              double newCutoff = model2->getMinimizationObjValue()-
1416                model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
1417              if (model2->getSolutionCount()) {
1418                if (!model2->status())
1419                  assert (model2->isProvenOptimal());
1420                memcpy(model_->bestSolution(),subModel->bestSolution(),
1421                       numberColumns*sizeof(double));
1422                model_->setCutoff(newCutoff);
1423              }
1424              delete model2;
1425            } else {
1426              // infeasible - could just be - due to cutoff
1427            }
1428          } else {
1429            // too dangerous at present
1430            assert (subModel->typePresolve()!=2);
1431          }
1432          if(model_->getCutoff()<bestCutoff_) {
1433            // Save solution
1434            if (!bestSolution_) 
1435              bestSolution_ = new double [numberColumns];
1436            memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
1437            bestCutoff_=model_->getCutoff();
1438          }
1439          delete subModel;
1440        }
1441        // we have done search to make sure best general solution
1442        searchType_=1;
1443        // Reverse cut weakly
1444        reverseCut(3,rhs_);
1445      } else {
1446        searchType_=1;
1447        // delete last cut
1448        deleteCut(cut_);
1449      }
1450    } else {
1451      searchType_=1;
1452    }
1453    // save best solution in this subtree
1454    memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
1455    nextStrong_=false;
1456    rhs_=range_;
1457    break;
1458  case 4:
1459    // solution not found and subtree not exhausted
1460    if (maxDiversification_) {
1461      if (nextStrong_) {
1462        // Reverse cut weakly
1463        reverseCut(4,rhs_);
1464        model_->setCutoff(1.0e50);
1465        diversification_++;
1466        searchType_=0;
1467      } else {
1468        // delete last cut
1469        deleteCut(cut_);
1470        searchType_=1;
1471      }
1472      nextStrong_=true;
1473      rhs_ += range_/2;
1474    } else {
1475      // special case when using as heuristic
1476      // Reverse cut weakly if lb -infinity
1477      reverseCut(4,rhs_);
1478      // This will be last try (may hit max time0
1479      lastTry=true;
1480      model_->setCutoff(bestCutoff_);
1481      if (model_->messageHandler()->logLevel()>0)
1482        printf("Exiting local search with current set of cuts\n");
1483      rhs_=1.0e100;
1484      // Can now stop on gap
1485      model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
1486      typeCuts_ =-1;
1487    }
1488    break;
1489  }
1490  if (rhs_<1.0e30||lastTry) {
1491    int goodSolution=createCut(savedSolution_,cut_);
1492    if (goodSolution>=0) {
1493      // Add to global cuts
1494      model_->globalCuts()->insert(cut_);
1495      OsiCuts * global = model_->globalCuts();
1496      int n = global->sizeRowCuts();
1497      OsiRowCut * rowCut = global->rowCutPtr(n-1);
1498      if (model_->messageHandler()->logLevel()>0)
1499        printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
1500               n,rowCut->lb(),rowCut->ub(),rhs_,diversification_);
1501      const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
1502      if (debugger) {
1503        if(debugger->invalidCut(*rowCut))
1504          printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1505      }
1506      for (int i=0;i<n;i++) {
1507        rowCut = global->rowCutPtr(i);
1508        if (model_->messageHandler()->logLevel()>0)
1509          printf("%d - rhs %g %g\n",
1510                 i,rowCut->lb(),rowCut->ub());
1511      }
1512    }
1513    // put back node
1514    startTime_ = static_cast<int> (CoinCpuTime());
1515    startNode_=model_->getNodeCount();
1516    if (localNode_) {
1517      // save copy of node
1518      CbcNode * localNode2 = new CbcNode(*localNode_);
1519      // But localNode2 now owns cuts so swap
1520      //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
1521      //   localNode_,localNode_->nodeInfo());
1522      nodes_.push_back(localNode_);
1523      localNode_=localNode2;
1524      std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
1525    }
1526  }
1527  return finished;
1528}
1529// We may have got an intelligent tree so give it one more chance
1530void 
1531CbcTreeVariable::endSearch() 
1532{
1533  if (typeCuts_>=0) {
1534    // copy best solution to model
1535    int numberColumns = model_->getNumCols();
1536    if (bestSolution_&&bestCutoff_<model_->getCutoff()) {
1537      memcpy(model_->bestSolution(),bestSolution_,numberColumns*sizeof(double));
1538      model_->setCutoff(bestCutoff_);
1539      // recompute objective value
1540      const double * objCoef = model_->getObjCoefficients();
1541      double objOffset=0.0;
1542      model_->continuousSolver()->getDblParam(OsiObjOffset,objOffset);
1543 
1544      // Compute dot product of objCoef and colSol and then adjust by offset
1545      double objValue = -objOffset;
1546      for ( int i=0 ; i<numberColumns ; i++ )
1547        objValue += objCoef[i]*bestSolution_[i];
1548      model_->setMinimizationObjValue(objValue);
1549    }
1550    // Can now stop on gap
1551    model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
1552  }
1553}
1554// Create cut
1555int
1556CbcTreeVariable::createCut(const double * solution,OsiRowCut & rowCut)
1557{
1558  if (rhs_>1.0e20)
1559    return -1;
1560  OsiSolverInterface * solver = model_->solver();
1561  const double * rowLower = solver->getRowLower();
1562  const double * rowUpper = solver->getRowUpper();
1563  //const double * solution = solver->getColSolution();
1564  //const double * objective = solver->getObjCoefficients();
1565  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1566  double primalTolerance;
1567  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
1568  // relax
1569  primalTolerance *= 1000.0;
1570
1571  int numberRows = model_->getNumRows();
1572
1573  int numberIntegers = model_->numberIntegers();
1574  const int * integerVariable = model_->integerVariable();
1575  int i;
1576
1577  // Check feasible
1578
1579  double * rowActivity = new double[numberRows];
1580  memset(rowActivity,0,numberRows*sizeof(double));
1581  solver->getMatrixByCol()->times(solution,rowActivity) ;
1582  int goodSolution=0;
1583  // check was feasible
1584  for (i=0;i<numberRows;i++) {
1585    if(rowActivity[i]<rowLower[i]-primalTolerance) {
1586      goodSolution=-1;
1587    } else if(rowActivity[i]>rowUpper[i]+primalTolerance) {
1588      goodSolution=-1;
1589    }
1590  }
1591  delete [] rowActivity;
1592  for (i=0;i<numberIntegers;i++) {
1593    int iColumn = integerVariable[i];
1594    double value=solution[iColumn];
1595    if (fabs(floor(value+0.5)-value)>integerTolerance) {
1596      goodSolution=-1;
1597    }
1598  }
1599  // zap cut
1600  if (goodSolution==0) {
1601    // Create cut and get total gap
1602    CoinPackedVector cut;
1603    double rhs = rhs_;
1604    double maxValue=0.0;
1605    for (i=0;i<numberIntegers;i++) {
1606      int iColumn = integerVariable[i];
1607      double value = floor(solution[iColumn]+0.5);
1608      if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
1609        continue; // skip as not 0-1
1610      if (originalLower_[i]==originalUpper_[i])
1611        continue;
1612      double mu =1.0/(originalUpper_[i]-originalLower_[i]);
1613      if (value==originalLower_[i]) {
1614        rhs += mu*originalLower_[i];
1615        cut.insert(iColumn,1.0);
1616        maxValue += originalUpper_[i];
1617      } else if (value==originalUpper_[i]) {
1618        rhs -= mu*originalUpper_[i];
1619        cut.insert(iColumn,-1.0);
1620        maxValue += originalLower_[i];
1621      }
1622    }
1623    if (maxValue<rhs-primalTolerance) {
1624      if (model_->messageHandler()->logLevel()>0)
1625        printf("slack cut\n");
1626      goodSolution=1;
1627    }
1628    rowCut.setRow(cut);
1629    rowCut.setLb(-COIN_DBL_MAX);
1630    rowCut.setUb(rhs);   
1631    rowCut.setGloballyValid();
1632    if (model_->messageHandler()->logLevel()>0)
1633    printf("Cut size: %i Cut rhs: %g\n",cut.getNumElements(), rhs);
1634#ifdef CBC_DEBUG
1635    if (model_->messageHandler()->logLevel()>0) {
1636      int k;
1637      for (k=0; k<cut.getNumElements(); k++){
1638        printf("%i    %g ",cut.getIndices()[k], cut.getElements()[k]);
1639        if ((k+1)%5==0)
1640          printf("\n");
1641      }
1642      if (k%5!=0)
1643        printf("\n");
1644    }
1645#endif
1646    return goodSolution;
1647  } else {
1648    if (model_->messageHandler()->logLevel()>0)
1649      printf("Not a good solution\n");
1650    return -1;
1651  }
1652}
1653// Other side of last cut branch
1654void 
1655CbcTreeVariable::reverseCut(int state, double bias)
1656{
1657  // find global cut
1658  OsiCuts * global = model_->globalCuts();
1659  int n = global->sizeRowCuts();
1660  int i;
1661  OsiRowCut * rowCut = NULL;
1662  for ( i=0;i<n;i++) {
1663    rowCut = global->rowCutPtr(i);
1664    if (cut_==*rowCut) {
1665      break;
1666    }
1667  }
1668  if (!rowCut) {
1669    // must have got here in odd way e.g. strong branching
1670    return;
1671  }
1672  if (rowCut->lb()>-1.0e10)
1673    return;
1674  // get smallest element
1675  double smallest=COIN_DBL_MAX;
1676  CoinPackedVector row = cut_.row();
1677  for (int k=0; k<row.getNumElements(); k++)
1678    smallest = CoinMin(smallest,fabs(row.getElements()[k]));
1679  if (!typeCuts_&&!refine_) {
1680    // Reverse cut very very weakly
1681    if (state>2)
1682      smallest=0.0;
1683  }
1684  // replace by other way
1685  if (model_->messageHandler()->logLevel()>0)
1686    printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
1687           i,n,rowCut->lb(),rowCut->ub());
1688  rowCut->setLb(rowCut->ub()+smallest-bias);
1689  rowCut->setUb(COIN_DBL_MAX);
1690  if (model_->messageHandler()->logLevel()>0)
1691    printf("new rhs %g %g, bias %g smallest %g ",
1692           rowCut->lb(),rowCut->ub(),bias,smallest);
1693  const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
1694  if (debugger) {
1695    if(debugger->invalidCut(*rowCut))
1696      printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1697  }
1698}
1699// Delete last cut branch
1700void 
1701CbcTreeVariable::deleteCut(OsiRowCut & cut)
1702{
1703  // find global cut
1704  OsiCuts * global = model_->globalCuts();
1705  int n = global->sizeRowCuts();
1706  int i;
1707  OsiRowCut * rowCut = NULL;
1708  for ( i=0;i<n;i++) {
1709    rowCut = global->rowCutPtr(i);
1710    if (cut==*rowCut) {
1711      break;
1712    }
1713  }
1714  assert (i<n);
1715  // delete last cut
1716  if (model_->messageHandler()->logLevel()>0)
1717    printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
1718           i,n,rowCut->lb(),rowCut->ub());
1719  global->eraseRowCut(i);
1720}
1721// Create C++ lines to get to current state
1722void 
1723CbcTreeVariable::generateCpp( FILE * fp) 
1724{
1725  CbcTreeVariable other;
1726  fprintf(fp,"0#include \"CbcTreeVariable.hpp\"\n");
1727  fprintf(fp,"5  CbcTreeVariable variableTree(cbcModel,NULL);\n");
1728  if (range_!=other.range_)
1729    fprintf(fp,"5  variableTree.setRange(%d);\n",range_);
1730  if (typeCuts_!=other.typeCuts_)
1731    fprintf(fp,"5  variableTree.setTypeCuts(%d);\n",typeCuts_);
1732  if (maxDiversification_!=other.maxDiversification_)
1733    fprintf(fp,"5  variableTree.setMaxDiversification(%d);\n",maxDiversification_);
1734  if (timeLimit_!=other.timeLimit_)
1735    fprintf(fp,"5  variableTree.setTimeLimit(%d);\n",timeLimit_);
1736  if (nodeLimit_!=other.nodeLimit_)
1737    fprintf(fp,"5  variableTree.setNodeLimit(%d);\n",nodeLimit_);
1738  if (refine_!=other.refine_)
1739    fprintf(fp,"5  variableTree.setRefine(%s);\n",refine_ ? "true" : "false");
1740  fprintf(fp,"5  cbcModel->passInTreeHandler(variableTree);\n");
1741}
1742
1743
1744
Note: See TracBrowser for help on using the repository browser.