source: trunk/CbcTreeLocal.cpp @ 188

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

local tree

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.8 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CbcModel.hpp"
5#include "CbcNode.hpp"
6#include "CbcTreeLocal.hpp"
7#include "CoinPackedMatrix.hpp"
8#include "CoinTime.hpp"
9#include <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          // But adjust
504          subModel->setCutoff(model_->getCutoff()+model_->getDblParam(CbcModel::CbcCutoffIncrement)+1.0e-6);
505          subModel->setSolutionCount(0);
506          assert (subModel->isProvenOptimal());
507          if (!subModel->typePresolve()) {
508            subModel->branchAndBound();
509            if (subModel->status()) {
510              model_->incrementSubTreeStopped(); 
511            }
512            //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
513            //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
514            double newCutoff = subModel->getMinimizationObjValue()-
515              subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
516            if (subModel->getSolutionCount()) {
517              if (!subModel->status())
518                assert (subModel->isProvenOptimal());
519              memcpy(model_->bestSolution(),subModel->bestSolution(),
520                     numberColumns*sizeof(double));
521              model_->setCutoff(newCutoff);
522            }
523          } else if (subModel->typePresolve()==1) {
524            CbcModel * model2 = subModel->integerPresolve(true);
525            if (model2) {
526              // Do complete search
527              model2->branchAndBound();
528              // get back solution
529              subModel->originalModel(model2,false);
530              if (model2->status()) {
531                model_->incrementSubTreeStopped(); 
532              } 
533              double newCutoff = model2->getMinimizationObjValue()-
534                model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
535              if (model2->getSolutionCount()) {
536                if (!model2->status())
537                  assert (model2->isProvenOptimal());
538                memcpy(model_->bestSolution(),subModel->bestSolution(),
539                       numberColumns*sizeof(double));
540                model_->setCutoff(newCutoff);
541              }
542              delete model2;
543            } else {
544              // infeasible - could just be - due to cutoff
545            }
546          } else {
547            // too dangerous at present
548            assert (subModel->typePresolve()!=2);
549          }
550          if(model_->getCutoff()<bestCutoff_) {
551            // Save solution
552            if (!bestSolution_) 
553              bestSolution_ = new double [numberColumns];
554            memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
555            bestCutoff_=model_->getCutoff();
556          }
557          delete subModel;
558        }
559        // we have done search to make sure best general solution
560        searchType_=1;
561        // Reverse cut weakly
562        reverseCut(3,rhs_);
563      } else {
564        searchType_=1;
565        // delete last cut
566        deleteCut(cut_);
567      }
568    } else {
569      searchType_=1;
570    }
571    // save best solution in this subtree
572    memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
573    nextStrong_=false;
574    rhs_=range_;
575    break;
576  case 4:
577    // solution not found and subtree not exhausted
578    if (nextStrong_) {
579      // Reverse cut weakly
580      reverseCut(4,rhs_);
581      model_->setCutoff(1.0e50);
582      diversification_++;
583      searchType_=0;
584    } else {
585      // delete last cut
586      deleteCut(cut_);
587      searchType_=1;
588    }
589    nextStrong_=true;
590    rhs_ += range_/2;
591    break;
592  }
593  if (rhs_<1.0e30||lastTry) {
594    int goodSolution=createCut(savedSolution_,cut_);
595    if (goodSolution>=0) {
596      // Add to global cuts
597      model_->globalCuts()->insert(cut_);
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    if (localNode_) {
613      // save copy of node
614      CbcNode * localNode2 = new CbcNode(*localNode_);
615      // But localNode2 now owns cuts so swap
616      //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
617      //   localNode_,localNode_->nodeInfo());
618      nodes_.push_back(localNode_);
619      localNode_=localNode2;
620      make_heap(nodes_.begin(), nodes_.end(), comparison_);
621    }
622  }
623  return finished;
624}
625// We may have got an intelligent tree so give it one more chance
626void 
627CbcTreeLocal::endSearch() 
628{
629  if (typeCuts_>=0) {
630    // copy best solution to model
631    int numberColumns = model_->getNumCols();
632    if (bestSolution_&&bestCutoff_<model_->getCutoff()) {
633      memcpy(model_->bestSolution(),bestSolution_,numberColumns*sizeof(double));
634      model_->setCutoff(bestCutoff_);
635      // recompute objective value
636      const double * objCoef = model_->getObjCoefficients();
637      double objOffset=0.0;
638      model_->continuousSolver()->getDblParam(OsiObjOffset,objOffset);
639 
640      // Compute dot product of objCoef and colSol and then adjust by offset
641      double objValue = -objOffset;
642      for ( int i=0 ; i<numberColumns ; i++ )
643        objValue += objCoef[i]*bestSolution_[i];
644      model_->setMinimizationObjValue(objValue);
645    }
646    // Can now stop on gap
647    model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
648  }
649}
650// Create cut
651int
652CbcTreeLocal::createCut(const double * solution,OsiRowCut & rowCut)
653{
654  OsiSolverInterface * solver = model_->solver();
655  const double * rowLower = solver->getRowLower();
656  const double * rowUpper = solver->getRowUpper();
657  //const double * solution = solver->getColSolution();
658  //const double * objective = solver->getObjCoefficients();
659  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
660  double primalTolerance;
661  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
662  // relax
663  primalTolerance *= 1000.0;
664
665  int numberRows = model_->getNumRows();
666
667  int numberIntegers = model_->numberIntegers();
668  const int * integerVariable = model_->integerVariable();
669  int i;
670
671  // Check feasible
672
673  double * rowActivity = new double[numberRows];
674  memset(rowActivity,0,numberRows*sizeof(double));
675  solver->getMatrixByCol()->times(solution,rowActivity) ;
676  int goodSolution=0;
677  // check was feasible
678  for (i=0;i<numberRows;i++) {
679    if(rowActivity[i]<rowLower[i]-primalTolerance) {
680      goodSolution=-1;
681    } else if(rowActivity[i]>rowUpper[i]+primalTolerance) {
682      goodSolution=-1;
683    }
684  }
685  delete [] rowActivity;
686  for (i=0;i<numberIntegers;i++) {
687    int iColumn = integerVariable[i];
688    double value=solution[iColumn];
689    if (fabs(floor(value+0.5)-value)>integerTolerance) {
690      goodSolution=-1;
691    }
692  }
693  // zap cut
694  if (goodSolution==0) {
695    // Create cut and get total gap
696    CoinPackedVector cut;
697    double rhs = rhs_;
698    double maxValue=0.0;
699    for (i=0;i<numberIntegers;i++) {
700      int iColumn = integerVariable[i];
701      double value = floor(solution[iColumn]+0.5);
702      if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
703        continue; // skip as not 0-1
704      if (originalLower_[i]==originalUpper_[i])
705        continue;
706      double mu =1.0/(originalUpper_[i]-originalLower_[i]);
707      if (value==originalLower_[i]) {
708        rhs += mu*originalLower_[i];
709        cut.insert(iColumn,1.0);
710        maxValue += originalUpper_[i];
711      } else if (value==originalUpper_[i]) {
712        rhs -= mu*originalUpper_[i];
713        cut.insert(iColumn,-1.0);
714        maxValue += originalLower_[i];
715      }
716    }
717    if (maxValue<rhs-primalTolerance) {
718      printf("slack cut\n");
719      goodSolution=1;
720    }
721    rowCut.setRow(cut);
722    rowCut.setLb(-COIN_DBL_MAX);
723    rowCut.setUb(rhs);   
724    rowCut.setGloballyValid();
725    printf("Cut size: %i Cut rhs: %g\n",cut.getNumElements(), rhs);
726#ifdef CBC_DEBUG
727    int k;
728    for (k=0; k<cut.getNumElements(); k++){
729      printf("%i    %g ",cut.getIndices()[k], cut.getElements()[k]);
730      if ((k+1)%5==0)
731        printf("\n");
732    }
733    if (k%5!=0)
734      printf("\n");
735#endif
736    return goodSolution;
737  } else {
738    printf("Not a good solution\n");
739    return -1;
740  }
741}
742// Other side of last cut branch
743void 
744CbcTreeLocal::reverseCut(int state, double bias)
745{
746  // find global cut
747  OsiCuts * global = model_->globalCuts();
748  int n = global->sizeRowCuts();
749  int i;
750  OsiRowCut * rowCut = NULL;
751  for ( i=0;i<n;i++) {
752    rowCut = global->rowCutPtr(i);
753    if (cut_==*rowCut) {
754      break;
755    }
756  }
757  if (!rowCut) {
758    // must have got here in odd way e.g. strong branching
759    return;
760  }
761  // get smallest element
762  double smallest=COIN_DBL_MAX;
763  CoinPackedVector row = cut_.row();
764  for (int k=0; k<row.getNumElements(); k++)
765    smallest = min(smallest,fabs(row.getElements()[k]));
766  if (!typeCuts_&&!refine_) {
767    // Reverse cut very very weakly
768    if (state>2)
769      smallest=0.0;
770  }
771  // replace by other way
772  printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
773         i,n,rowCut->lb(),rowCut->ub());
774  rowCut->setLb(rowCut->ub()+smallest-bias);
775  rowCut->setUb(COIN_DBL_MAX);
776  printf("new rhs %g %g, bias %g smallest %g ",
777         rowCut->lb(),rowCut->ub(),bias,smallest);
778}
779// Delete last cut branch
780void 
781CbcTreeLocal::deleteCut(OsiRowCut & cut)
782{
783  // find global cut
784  OsiCuts * global = model_->globalCuts();
785  int n = global->sizeRowCuts();
786  int i;
787  OsiRowCut * rowCut = NULL;
788  for ( i=0;i<n;i++) {
789    rowCut = global->rowCutPtr(i);
790    if (cut==*rowCut) {
791      break;
792    }
793  }
794  assert (i<n);
795  // delete last cut
796  printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
797         i,n,rowCut->lb(),rowCut->ub());
798  global->eraseRowCut(i);
799}
800
801
802
Note: See TracBrowser for help on using the repository browser.