source: trunk/CbcTreeLocal.cpp @ 186

Last change on this file since 186 was 186, checked in by forrest, 16 years ago

heuristics

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