source: trunk/CbcTreeLocal.cpp @ 216

Last change on this file since 216 was 197, checked in by forrest, 15 years ago

stuff

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