source: trunk/CbcTreeLocal.cpp @ 195

Last change on this file since 195 was 195, checked in by forrest, 14 years ago

stuff

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