source: trunk/Cbc/src/CbcTreeLocal.cpp

Last change on this file was 2467, checked in by unxusr, 7 months ago

spaces after angles

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 57.2 KB
Line 
1/* $Id: CbcTreeLocal.cpp 2467 2019-01-03 21:26:29Z unxusr $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CbcModel.hpp"
7#include "CbcNode.hpp"
8#include "CbcTreeLocal.hpp"
9#include "CoinPackedMatrix.hpp"
10#include "CoinTime.hpp"
11#include "OsiRowCutDebugger.hpp"
12#include <cassert>
13#ifdef JJF_ZERO
14// gdb doesn't always put breakpoints in this virtual function
15// just stick xxxxxx() where you want to start
16static void xxxxxx()
17{
18  printf("break\n");
19}
20#endif
21CbcTreeLocal::CbcTreeLocal()
22  : localNode_(NULL)
23  , bestSolution_(NULL)
24  , savedSolution_(NULL)
25  , saveNumberSolutions_(0)
26  , model_(NULL)
27  , originalLower_(NULL)
28  , originalUpper_(NULL)
29  , range_(0)
30  , typeCuts_(-1)
31  , maxDiversification_(0)
32  , diversification_(0)
33  , nextStrong_(false)
34  , rhs_(0.0)
35  , savedGap_(0.0)
36  , bestCutoff_(0.0)
37  , timeLimit_(0)
38  , startTime_(0)
39  , nodeLimit_(0)
40  , startNode_(-1)
41  , searchType_(-1)
42  , refine_(false)
43{
44}
45/* Constructor with solution.
46   range is upper bound on difference from given solution.
47   maxDiversification is maximum number of diversifications to try
48   timeLimit is seconds in subTree
49   nodeLimit is nodes in subTree
50*/
51CbcTreeLocal::CbcTreeLocal(CbcModel *model, const double *solution,
52  int range, int typeCuts, int maxDiversification,
53  int timeLimit, int nodeLimit, bool refine)
54  : localNode_(NULL)
55  , bestSolution_(NULL)
56  , savedSolution_(NULL)
57  , saveNumberSolutions_(0)
58  , model_(model)
59  , originalLower_(NULL)
60  , originalUpper_(NULL)
61  , range_(range)
62  , typeCuts_(typeCuts)
63  , maxDiversification_(maxDiversification)
64  , diversification_(0)
65  , nextStrong_(false)
66  , rhs_(0.0)
67  , savedGap_(0.0)
68  , bestCutoff_(0.0)
69  , timeLimit_(timeLimit)
70  , startTime_(0)
71  , nodeLimit_(nodeLimit)
72  , startNode_(-1)
73  , searchType_(-1)
74  , refine_(refine)
75{
76
77  OsiSolverInterface *solver = model_->solver();
78  const double *lower = solver->getColLower();
79  const double *upper = solver->getColUpper();
80  //const double * solution = solver->getColSolution();
81  //const double * objective = solver->getObjCoefficients();
82  double primalTolerance;
83  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
84
85  // Get increment
86  model_->analyzeObjective();
87
88  {
89    // needed to sync cutoffs
90    double value;
91    solver->getDblParam(OsiDualObjectiveLimit, value);
92    model_->setCutoff(value * solver->getObjSense());
93  }
94  bestCutoff_ = model_->getCutoff();
95  // save current gap
96  savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
97
98  // make sure integers found
99  model_->findIntegers(false);
100  int numberIntegers = model_->numberIntegers();
101  const int *integerVariable = model_->integerVariable();
102  int i;
103  double direction = solver->getObjSense();
104  double newSolutionValue = 1.0e50;
105  if (solution) {
106    // copy solution
107    solver->setColSolution(solution);
108    newSolutionValue = direction * solver->getObjValue();
109  }
110  originalLower_ = new double[numberIntegers];
111  originalUpper_ = new double[numberIntegers];
112  bool all01 = true;
113  int number01 = 0;
114  for (i = 0; i < numberIntegers; i++) {
115    int iColumn = integerVariable[i];
116    originalLower_[i] = lower[iColumn];
117    originalUpper_[i] = upper[iColumn];
118    if (upper[iColumn] - lower[iColumn] > 1.5)
119      all01 = false;
120    else if (upper[iColumn] - lower[iColumn] == 1.0)
121      number01++;
122  }
123  if (all01 && !typeCuts_)
124    typeCuts_ = 1; // may as well so we don't have to deal with refine
125  if (!number01 && !typeCuts_) {
126    if (model_->messageHandler()->logLevel() > 1)
127      printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
128    typeCuts_ = -1;
129  } else {
130    if (model_->messageHandler()->logLevel() > 1) {
131      std::string type;
132      if (all01) {
133        printf("%d 0-1 variables normal local  cuts\n",
134          number01);
135      } else if (typeCuts_) {
136        printf("%d 0-1 variables, %d other - general integer local cuts\n",
137          number01, numberIntegers - number01);
138      } else {
139        printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
140          number01, numberIntegers - number01);
141      }
142      printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
143        maxDiversification_, range_, timeLimit_, nodeLimit_);
144    }
145  }
146  int numberColumns = model_->getNumCols();
147  savedSolution_ = new double[numberColumns];
148  memset(savedSolution_, 0, numberColumns * sizeof(double));
149  if (solution) {
150    rhs_ = range_;
151    // Check feasible
152    int goodSolution = createCut(solution, cut_);
153    if (goodSolution >= 0) {
154      for (i = 0; i < numberIntegers; i++) {
155        int iColumn = integerVariable[i];
156        double value = floor(solution[iColumn] + 0.5);
157        // fix so setBestSolution will work
158        solver->setColLower(iColumn, value);
159        solver->setColUpper(iColumn, value);
160      }
161      model_->reserveCurrentSolution();
162      // Create cut and get total gap
163      if (newSolutionValue < bestCutoff_) {
164        model_->setBestSolution(CBC_ROUNDING, newSolutionValue, solution);
165        bestCutoff_ = model_->getCutoff();
166        // save as best solution
167        memcpy(savedSolution_, model_->bestSolution(), numberColumns * sizeof(double));
168      }
169      for (i = 0; i < numberIntegers; i++) {
170        int iColumn = integerVariable[i];
171        // restore bounds
172        solver->setColLower(iColumn, originalLower_[i]);
173        solver->setColUpper(iColumn, originalUpper_[i]);
174      }
175      // make sure can't stop on gap
176      model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
177    } else {
178      model_ = NULL;
179    }
180  } else {
181    // no solution
182    rhs_ = 1.0e50;
183    // make sure can't stop on gap
184    model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
185  }
186}
187CbcTreeLocal::~CbcTreeLocal()
188{
189  delete[] originalLower_;
190  delete[] originalUpper_;
191  delete[] bestSolution_;
192  delete[] savedSolution_;
193  delete localNode_;
194}
195// Copy constructor
196CbcTreeLocal::CbcTreeLocal(const CbcTreeLocal &rhs)
197  : CbcTree(rhs)
198  , saveNumberSolutions_(rhs.saveNumberSolutions_)
199  , model_(rhs.model_)
200  , range_(rhs.range_)
201  , typeCuts_(rhs.typeCuts_)
202  , maxDiversification_(rhs.maxDiversification_)
203  , diversification_(rhs.diversification_)
204  , nextStrong_(rhs.nextStrong_)
205  , rhs_(rhs.rhs_)
206  , savedGap_(rhs.savedGap_)
207  , bestCutoff_(rhs.bestCutoff_)
208  , timeLimit_(rhs.timeLimit_)
209  , startTime_(rhs.startTime_)
210  , nodeLimit_(rhs.nodeLimit_)
211  , startNode_(rhs.startNode_)
212  , searchType_(rhs.searchType_)
213  , refine_(rhs.refine_)
214{
215  cut_ = rhs.cut_;
216  fixedCut_ = rhs.fixedCut_;
217  if (rhs.localNode_)
218    localNode_ = new CbcNode(*rhs.localNode_);
219  else
220    localNode_ = NULL;
221  if (rhs.originalLower_) {
222    int numberIntegers = model_->numberIntegers();
223    originalLower_ = new double[numberIntegers];
224    memcpy(originalLower_, rhs.originalLower_, numberIntegers * sizeof(double));
225    originalUpper_ = new double[numberIntegers];
226    memcpy(originalUpper_, rhs.originalUpper_, numberIntegers * sizeof(double));
227  } else {
228    originalLower_ = NULL;
229    originalUpper_ = NULL;
230  }
231  if (rhs.bestSolution_) {
232    int numberColumns = model_->getNumCols();
233    bestSolution_ = new double[numberColumns];
234    memcpy(bestSolution_, rhs.bestSolution_, numberColumns * sizeof(double));
235  } else {
236    bestSolution_ = NULL;
237  }
238  if (rhs.savedSolution_) {
239    int numberColumns = model_->getNumCols();
240    savedSolution_ = new double[numberColumns];
241    memcpy(savedSolution_, rhs.savedSolution_, numberColumns * sizeof(double));
242  } else {
243    savedSolution_ = NULL;
244  }
245}
246//----------------------------------------------------------------
247// Assignment operator
248//-------------------------------------------------------------------
249CbcTreeLocal &
250CbcTreeLocal::operator=(const CbcTreeLocal &rhs)
251{
252  if (this != &rhs) {
253    CbcTree::operator=(rhs);
254    saveNumberSolutions_ = rhs.saveNumberSolutions_;
255    cut_ = rhs.cut_;
256    fixedCut_ = rhs.fixedCut_;
257    delete localNode_;
258    if (rhs.localNode_)
259      localNode_ = new CbcNode(*rhs.localNode_);
260    else
261      localNode_ = NULL;
262    model_ = rhs.model_;
263    range_ = rhs.range_;
264    typeCuts_ = rhs.typeCuts_;
265    maxDiversification_ = rhs.maxDiversification_;
266    diversification_ = rhs.diversification_;
267    nextStrong_ = rhs.nextStrong_;
268    rhs_ = rhs.rhs_;
269    savedGap_ = rhs.savedGap_;
270    bestCutoff_ = rhs.bestCutoff_;
271    timeLimit_ = rhs.timeLimit_;
272    startTime_ = rhs.startTime_;
273    nodeLimit_ = rhs.nodeLimit_;
274    startNode_ = rhs.startNode_;
275    searchType_ = rhs.searchType_;
276    refine_ = rhs.refine_;
277    delete[] originalLower_;
278    delete[] originalUpper_;
279    if (rhs.originalLower_) {
280      int numberIntegers = model_->numberIntegers();
281      originalLower_ = new double[numberIntegers];
282      memcpy(originalLower_, rhs.originalLower_, numberIntegers * sizeof(double));
283      originalUpper_ = new double[numberIntegers];
284      memcpy(originalUpper_, rhs.originalUpper_, numberIntegers * sizeof(double));
285    } else {
286      originalLower_ = NULL;
287      originalUpper_ = NULL;
288    }
289    delete[] bestSolution_;
290    if (rhs.bestSolution_) {
291      int numberColumns = model_->getNumCols();
292      bestSolution_ = new double[numberColumns];
293      memcpy(bestSolution_, rhs.bestSolution_, numberColumns * sizeof(double));
294    } else {
295      bestSolution_ = NULL;
296    }
297    delete[] savedSolution_;
298    if (rhs.savedSolution_) {
299      int numberColumns = model_->getNumCols();
300      savedSolution_ = new double[numberColumns];
301      memcpy(savedSolution_, rhs.savedSolution_, numberColumns * sizeof(double));
302    } else {
303      savedSolution_ = NULL;
304    }
305  }
306  return *this;
307}
308// Clone
309CbcTree *
310CbcTreeLocal::clone() const
311{
312  return new CbcTreeLocal(*this);
313}
314// Pass in solution (so can be used after heuristic)
315void CbcTreeLocal::passInSolution(const double *solution, double solutionValue)
316{
317  int numberColumns = model_->getNumCols();
318  delete[] savedSolution_;
319  savedSolution_ = new double[numberColumns];
320  memcpy(savedSolution_, solution, numberColumns * sizeof(double));
321  rhs_ = range_;
322  // Check feasible
323  int goodSolution = createCut(solution, cut_);
324  if (goodSolution >= 0) {
325    bestCutoff_ = CoinMin(solutionValue, model_->getCutoff());
326  } else {
327    model_ = NULL;
328  }
329}
330// Return the top node of the heap
331CbcNode *
332CbcTreeLocal::top() const
333{
334#ifdef CBC_DEBUG
335  int smallest = 9999999;
336  int largest = -1;
337  double smallestD = 1.0e30;
338  double largestD = -1.0e30;
339  int n = nodes_.size();
340  for (int i = 0; i < n; i++) {
341    int nn = nodes_[i]->nodeInfo()->nodeNumber();
342    double dd = nodes_[i]->objectiveValue();
343    largest = CoinMax(largest, nn);
344    smallest = CoinMin(smallest, nn);
345    largestD = CoinMax(largestD, dd);
346    smallestD = CoinMin(smallestD, dd);
347  }
348  if (model_->messageHandler()->logLevel() > 1) {
349    printf("smallest %d, largest %d, top %d\n", smallest, largest,
350      nodes_.front()->nodeInfo()->nodeNumber());
351    printf("smallestD %g, largestD %g, top %g\n", smallestD, largestD, nodes_.front()->objectiveValue());
352  }
353#endif
354  return nodes_.front();
355}
356
357// Add a node to the heap
358void CbcTreeLocal::push(CbcNode *x)
359{
360  if (typeCuts_ >= 0 && !nodes_.size() && searchType_ < 0) {
361    startNode_ = model_->getNodeCount();
362    // save copy of node
363    localNode_ = new CbcNode(*x);
364
365    if (cut_.row().getNumElements()) {
366      // Add to global cuts
367      // we came in with solution
368      model_->makeGlobalCut(cut_);
369      if (model_->messageHandler()->logLevel() > 1)
370        printf("initial cut - rhs %g %g\n",
371          cut_.lb(), cut_.ub());
372      searchType_ = 1;
373    } else {
374      // stop on first solution
375      searchType_ = 0;
376    }
377    startTime_ = static_cast< int >(CoinCpuTime());
378    saveNumberSolutions_ = model_->getSolutionCount();
379  }
380  nodes_.push_back(x);
381#ifdef CBC_DEBUG
382  if (model_->messageHandler()->logLevel() > 0)
383    printf("pushing node onto heap %d %x %x\n",
384      x->nodeInfo()->nodeNumber(), x, x->nodeInfo());
385#endif
386  std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
387}
388
389// Remove the top node from the heap
390void CbcTreeLocal::pop()
391{
392  std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
393  nodes_.pop_back();
394}
395// Test if empty - does work if so
396bool CbcTreeLocal::empty()
397{
398  if (typeCuts_ < 0)
399    return !nodes_.size();
400  /* state -
401       0 iterating
402       1 subtree finished optimal solution for subtree found
403       2 subtree finished and no solution found
404       3 subtree exiting and solution found
405       4 subtree exiting and no solution found
406    */
407  int state = 0;
408  assert(searchType_ != 2);
409  if (searchType_) {
410    if (CoinCpuTime() - startTime_ > timeLimit_ || model_->getNodeCount() - startNode_ >= nodeLimit_) {
411      state = 4;
412    }
413  } else {
414    if (model_->getSolutionCount() > saveNumberSolutions_) {
415      state = 4;
416    }
417  }
418  if (!nodes_.size())
419    state = 2;
420  if (!state) {
421    return false;
422  }
423  // Finished this phase
424  int numberColumns = model_->getNumCols();
425  if (model_->getSolutionCount() > saveNumberSolutions_) {
426    if (model_->getCutoff() < bestCutoff_) {
427      // Save solution
428      if (!bestSolution_)
429        bestSolution_ = new double[numberColumns];
430      memcpy(bestSolution_, model_->bestSolution(), numberColumns * sizeof(double));
431      bestCutoff_ = model_->getCutoff();
432    }
433    state--;
434  }
435  // get rid of all nodes (safe even if already done)
436  double bestPossibleObjective;
437  cleanTree(model_, -COIN_DBL_MAX, bestPossibleObjective);
438
439  double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement);
440  if (model_->messageHandler()->logLevel() > 1)
441    printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
442      state,
443      model_->getNodeCount() - startNode_,
444      static_cast< int >(CoinCpuTime()) - startTime_,
445      model_->getCutoff() + increment, bestCutoff_ + increment, rhs_);
446  saveNumberSolutions_ = model_->getSolutionCount();
447  bool finished = false;
448  bool lastTry = false;
449  switch (state) {
450  case 1:
451    // solution found and subtree exhausted
452    if (rhs_ > 1.0e30) {
453      finished = true;
454    } else {
455      // find global cut and reverse
456      reverseCut(1);
457      searchType_ = 1; // first false
458      rhs_ = range_; // reset range
459      nextStrong_ = false;
460
461      // save best solution in this subtree
462      memcpy(savedSolution_, model_->bestSolution(), numberColumns * sizeof(double));
463    }
464    break;
465  case 2:
466    // solution not found and subtree exhausted
467    if (rhs_ > 1.0e30) {
468      finished = true;
469    } else {
470      // find global cut and reverse
471      reverseCut(2);
472      searchType_ = 1; // first false
473      if (diversification_ < maxDiversification_) {
474        if (nextStrong_) {
475          diversification_++;
476          // cut is valid so don't model_->setCutoff(1.0e50);
477          searchType_ = 0;
478        }
479        nextStrong_ = true;
480        rhs_ += range_ / 2;
481      } else {
482        // This will be last try (may hit max time)
483        lastTry = true;
484        if (!maxDiversification_)
485          typeCuts_ = -1; // make sure can't start again
486        model_->setCutoff(bestCutoff_);
487        if (model_->messageHandler()->logLevel() > 1)
488          printf("Exiting local search with current set of cuts\n");
489        rhs_ = 1.0e100;
490        // Can now stop on gap
491        model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
492      }
493    }
494    break;
495  case 3:
496    // solution found and subtree not exhausted
497    if (rhs_ < 1.0e30) {
498      if (searchType_) {
499        if (!typeCuts_ && refine_ && searchType_ == 1) {
500          // We need to check we have best solution given these 0-1 values
501          OsiSolverInterface *subSolver = model_->continuousSolver()->clone();
502          CbcModel *subModel = model_->subTreeModel(subSolver);
503          CbcTree normalTree;
504          subModel->passInTreeHandler(normalTree);
505          int numberIntegers = model_->numberIntegers();
506          const int *integerVariable = model_->integerVariable();
507          const double *solution = model_->bestSolution();
508          int i;
509          int numberColumns = model_->getNumCols();
510          for (i = 0; i < numberIntegers; i++) {
511            int iColumn = integerVariable[i];
512            double value = floor(solution[iColumn] + 0.5);
513            if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
514              continue; // skip as not 0-1
515            if (originalLower_[i] == originalUpper_[i])
516              continue;
517            subSolver->setColLower(iColumn, value);
518            subSolver->setColUpper(iColumn, value);
519          }
520          subSolver->initialSolve();
521          // We can copy cutoff
522          // But adjust
523          subModel->setCutoff(model_->getCutoff() + model_->getDblParam(CbcModel::CbcCutoffIncrement) + 1.0e-6);
524          subModel->setSolutionCount(0);
525          assert(subModel->isProvenOptimal());
526          if (!subModel->typePresolve()) {
527            subModel->branchAndBound();
528            if (subModel->status()) {
529              model_->incrementSubTreeStopped();
530            }
531            //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
532            //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
533            double newCutoff = subModel->getMinimizationObjValue() - subModel->getDblParam(CbcModel::CbcCutoffIncrement);
534            if (subModel->getSolutionCount()) {
535              if (!subModel->status())
536                assert(subModel->isProvenOptimal());
537              memcpy(model_->bestSolution(), subModel->bestSolution(),
538                numberColumns * sizeof(double));
539              model_->setCutoff(newCutoff);
540            }
541          } else if (subModel->typePresolve() == 1) {
542            CbcModel *model2 = subModel->integerPresolve(true);
543            if (model2) {
544              // Do complete search
545              model2->branchAndBound();
546              // get back solution
547              subModel->originalModel(model2, false);
548              if (model2->status()) {
549                model_->incrementSubTreeStopped();
550              }
551              double newCutoff = model2->getMinimizationObjValue() - model2->getDblParam(CbcModel::CbcCutoffIncrement);
552              if (model2->getSolutionCount()) {
553                if (!model2->status())
554                  assert(model2->isProvenOptimal());
555                memcpy(model_->bestSolution(), subModel->bestSolution(),
556                  numberColumns * sizeof(double));
557                model_->setCutoff(newCutoff);
558              }
559              delete model2;
560            } else {
561              // infeasible - could just be - due to cutoff
562            }
563          } else {
564            // too dangerous at present
565            assert(subModel->typePresolve() != 2);
566          }
567          if (model_->getCutoff() < bestCutoff_) {
568            // Save solution
569            if (!bestSolution_)
570              bestSolution_ = new double[numberColumns];
571            memcpy(bestSolution_, model_->bestSolution(), numberColumns * sizeof(double));
572            bestCutoff_ = model_->getCutoff();
573          }
574          delete subModel;
575        }
576        // we have done search to make sure best general solution
577        searchType_ = 1;
578        // Reverse cut weakly
579        reverseCut(3, rhs_);
580      } else {
581        searchType_ = 1;
582        // delete last cut
583        deleteCut(cut_);
584      }
585    } else {
586      searchType_ = 1;
587    }
588    // save best solution in this subtree
589    memcpy(savedSolution_, model_->bestSolution(), numberColumns * sizeof(double));
590    nextStrong_ = false;
591    rhs_ = range_;
592    break;
593  case 4:
594    // solution not found and subtree not exhausted
595    if (maxDiversification_) {
596      if (nextStrong_) {
597        // Reverse cut weakly
598        reverseCut(4, rhs_);
599        model_->setCutoff(1.0e50);
600        diversification_++;
601        searchType_ = 0;
602      } else {
603        // delete last cut
604        deleteCut(cut_);
605        searchType_ = 1;
606      }
607      nextStrong_ = true;
608      rhs_ += range_ / 2;
609    } else {
610      // special case when using as heuristic
611      // Reverse cut weakly if lb -infinity
612      reverseCut(4, rhs_);
613      // This will be last try (may hit max time0
614      lastTry = true;
615      model_->setCutoff(bestCutoff_);
616      if (model_->messageHandler()->logLevel() > 1)
617        printf("Exiting local search with current set of cuts\n");
618      rhs_ = 1.0e100;
619      // Can now stop on gap
620      model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
621      typeCuts_ = -1;
622    }
623    break;
624  }
625  if (rhs_ < 1.0e30 || lastTry) {
626    int goodSolution = createCut(savedSolution_, cut_);
627    if (goodSolution >= 0) {
628      // Add to global cuts
629      model_->makeGlobalCut(cut_);
630      CbcRowCuts *global = model_->globalCuts();
631      int n = global->sizeRowCuts();
632      OsiRowCut *rowCut = global->rowCutPtr(n - 1);
633      if (model_->messageHandler()->logLevel() > 1)
634        printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
635          n, rowCut->lb(), rowCut->ub(), rhs_, diversification_);
636      const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways();
637      if (debugger) {
638        if (debugger->invalidCut(*rowCut))
639          printf("ZZZZTree Global cut - cuts off optimal solution!\n");
640      }
641      for (int i = 0; i < n; i++) {
642        rowCut = global->rowCutPtr(i);
643        if (model_->messageHandler()->logLevel() > 1)
644          printf("%d - rhs %g %g\n",
645            i, rowCut->lb(), rowCut->ub());
646      }
647    }
648    // put back node
649    startTime_ = static_cast< int >(CoinCpuTime());
650    startNode_ = model_->getNodeCount();
651    if (localNode_) {
652      // save copy of node
653      CbcNode *localNode2 = new CbcNode(*localNode_);
654      // But localNode2 now owns cuts so swap
655      //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
656      //   localNode_,localNode_->nodeInfo());
657      nodes_.push_back(localNode_);
658      localNode_ = localNode2;
659      std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
660    }
661  }
662  return finished;
663}
664// We may have got an intelligent tree so give it one more chance
665void CbcTreeLocal::endSearch()
666{
667  if (typeCuts_ >= 0) {
668    // copy best solution to model
669    int numberColumns = model_->getNumCols();
670    if (bestSolution_ && bestCutoff_ < model_->getCutoff()) {
671      memcpy(model_->bestSolution(), bestSolution_, numberColumns * sizeof(double));
672      model_->setCutoff(bestCutoff_);
673      // recompute objective value
674      const double *objCoef = model_->getObjCoefficients();
675      double objOffset = 0.0;
676      model_->continuousSolver()->getDblParam(OsiObjOffset, objOffset);
677
678      // Compute dot product of objCoef and colSol and then adjust by offset
679      double objValue = -objOffset;
680      for (int i = 0; i < numberColumns; i++)
681        objValue += objCoef[i] * bestSolution_[i];
682      model_->setMinimizationObjValue(objValue);
683    }
684    // Can now stop on gap
685    model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
686  }
687}
688// Create cut
689int CbcTreeLocal::createCut(const double *solution, OsiRowCut &rowCut)
690{
691  if (rhs_ > 1.0e20)
692    return -1;
693  OsiSolverInterface *solver = model_->solver();
694  const double *rowLower = solver->getRowLower();
695  const double *rowUpper = solver->getRowUpper();
696  //const double * solution = solver->getColSolution();
697  //const double * objective = solver->getObjCoefficients();
698  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
699  double primalTolerance;
700  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
701  // relax
702  primalTolerance *= 1000.0;
703
704  int numberRows = model_->getNumRows();
705
706  int numberIntegers = model_->numberIntegers();
707  const int *integerVariable = model_->integerVariable();
708  int i;
709
710  // Check feasible
711
712  double *rowActivity = new double[numberRows];
713  memset(rowActivity, 0, numberRows * sizeof(double));
714  solver->getMatrixByCol()->times(solution, rowActivity);
715  int goodSolution = 0;
716  // check was feasible
717  for (i = 0; i < numberRows; i++) {
718    if (rowActivity[i] < rowLower[i] - primalTolerance) {
719      goodSolution = -1;
720    } else if (rowActivity[i] > rowUpper[i] + primalTolerance) {
721      goodSolution = -1;
722    }
723  }
724  delete[] rowActivity;
725  for (i = 0; i < numberIntegers; i++) {
726    int iColumn = integerVariable[i];
727    double value = solution[iColumn];
728    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
729      goodSolution = -1;
730    }
731  }
732  // zap cut
733  if (goodSolution == 0) {
734    // Create cut and get total gap
735    CoinPackedVector cut;
736    double rhs = rhs_;
737    double maxValue = 0.0;
738    for (i = 0; i < numberIntegers; i++) {
739      int iColumn = integerVariable[i];
740      double value = floor(solution[iColumn] + 0.5);
741      /*
742              typeCuts_ == 0 restricts to binary, 1 allows general integer. But we're
743              still restricted to being up against a bound. Consider: the notion is that
744              the cut restricts us to a k-neighbourhood. For binary variables, this
745              amounts to k variables which change value. For general integer, we could
746              end up with a single variable sucking up all of k (hence mu --- the
747              variable must swing to its other bound to look like a movement of 1).  For
748              variables in the middle of a range, we're talking about fabs(sol<j> - x<j>).
749            */
750      if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
751        continue; // skip as not 0-1
752      if (originalLower_[i] == originalUpper_[i])
753        continue;
754      double mu = 1.0 / (originalUpper_[i] - originalLower_[i]);
755      if (value == originalLower_[i]) {
756        rhs += mu * originalLower_[i];
757        cut.insert(iColumn, 1.0);
758        maxValue += originalUpper_[i];
759      } else if (value == originalUpper_[i]) {
760        rhs -= mu * originalUpper_[i];
761        cut.insert(iColumn, -1.0);
762        maxValue += originalLower_[i];
763      }
764    }
765    if (maxValue < rhs - primalTolerance) {
766      if (model_->messageHandler()->logLevel() > 1)
767        printf("slack cut\n");
768      goodSolution = 1;
769    }
770    rowCut.setRow(cut);
771    rowCut.setLb(-COIN_DBL_MAX);
772    rowCut.setUb(rhs);
773    rowCut.setGloballyValid();
774    if (model_->messageHandler()->logLevel() > 1)
775      printf("Cut size: %i Cut rhs: %g\n", cut.getNumElements(), rhs);
776#ifdef CBC_DEBUG
777    if (model_->messageHandler()->logLevel() > 0) {
778      int k;
779      for (k = 0; k < cut.getNumElements(); k++) {
780        printf("%i    %g ", cut.getIndices()[k], cut.getElements()[k]);
781        if ((k + 1) % 5 == 0)
782          printf("\n");
783      }
784      if (k % 5 != 0)
785        printf("\n");
786    }
787#endif
788    return goodSolution;
789  } else {
790    if (model_->messageHandler()->logLevel() > 1)
791      printf("Not a good solution\n");
792    return -1;
793  }
794}
795// Other side of last cut branch
796void CbcTreeLocal::reverseCut(int state, double bias)
797{
798  // find global cut
799  CbcRowCuts *global = model_->globalCuts();
800  int n = global->sizeRowCuts();
801  int i;
802  OsiRowCut *rowCut = NULL;
803  for (i = 0; i < n; i++) {
804    rowCut = global->rowCutPtr(i);
805    if (cut_ == *rowCut) {
806      break;
807    }
808  }
809  if (!rowCut) {
810    // must have got here in odd way e.g. strong branching
811    return;
812  }
813  if (rowCut->lb() > -1.0e10)
814    return;
815  // get smallest element
816  double smallest = COIN_DBL_MAX;
817  CoinPackedVector row = cut_.row();
818  for (int k = 0; k < row.getNumElements(); k++)
819    smallest = CoinMin(smallest, fabs(row.getElements()[k]));
820  if (!typeCuts_ && !refine_) {
821    // Reverse cut very very weakly
822    if (state > 2)
823      smallest = 0.0;
824  }
825  // replace by other way
826  if (model_->messageHandler()->logLevel() > 1)
827    printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
828      i, n, rowCut->lb(), rowCut->ub());
829  rowCut->setLb(rowCut->ub() + smallest - bias);
830  rowCut->setUb(COIN_DBL_MAX);
831  if (model_->messageHandler()->logLevel() > 1)
832    printf("new rhs %g %g, bias %g smallest %g ",
833      rowCut->lb(), rowCut->ub(), bias, smallest);
834  const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways();
835  if (debugger) {
836    if (debugger->invalidCut(*rowCut))
837      printf("ZZZZTree Global cut - cuts off optimal solution!\n");
838  }
839}
840// Delete last cut branch
841void CbcTreeLocal::deleteCut(OsiRowCut &cut)
842{
843  // find global cut
844  CbcRowCuts *global = model_->globalCuts();
845  int n = global->sizeRowCuts();
846  int i;
847  OsiRowCut *rowCut = NULL;
848  for (i = 0; i < n; i++) {
849    rowCut = global->rowCutPtr(i);
850    if (cut == *rowCut) {
851      break;
852    }
853  }
854  assert(i < n);
855  // delete last cut
856  if (model_->messageHandler()->logLevel() > 1)
857    printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
858      i, n, rowCut->lb(), rowCut->ub());
859  global->eraseRowCut(i);
860}
861// Create C++ lines to get to current state
862void CbcTreeLocal::generateCpp(FILE *fp)
863{
864  CbcTreeLocal other;
865  fprintf(fp, "0#include \"CbcTreeLocal.hpp\"\n");
866  fprintf(fp, "5  CbcTreeLocal localTree(cbcModel,NULL);\n");
867  if (range_ != other.range_)
868    fprintf(fp, "5  localTree.setRange(%d);\n", range_);
869  if (typeCuts_ != other.typeCuts_)
870    fprintf(fp, "5  localTree.setTypeCuts(%d);\n", typeCuts_);
871  if (maxDiversification_ != other.maxDiversification_)
872    fprintf(fp, "5  localTree.setMaxDiversification(%d);\n", maxDiversification_);
873  if (timeLimit_ != other.timeLimit_)
874    fprintf(fp, "5  localTree.setTimeLimit(%d);\n", timeLimit_);
875  if (nodeLimit_ != other.nodeLimit_)
876    fprintf(fp, "5  localTree.setNodeLimit(%d);\n", nodeLimit_);
877  if (refine_ != other.refine_)
878    fprintf(fp, "5  localTree.setRefine(%s);\n", refine_ ? "true" : "false");
879  fprintf(fp, "5  cbcModel->passInTreeHandler(localTree);\n");
880}
881
882CbcTreeVariable::CbcTreeVariable()
883  : localNode_(NULL)
884  , bestSolution_(NULL)
885  , savedSolution_(NULL)
886  , saveNumberSolutions_(0)
887  , model_(NULL)
888  , originalLower_(NULL)
889  , originalUpper_(NULL)
890  , range_(0)
891  , typeCuts_(-1)
892  , maxDiversification_(0)
893  , diversification_(0)
894  , nextStrong_(false)
895  , rhs_(0.0)
896  , savedGap_(0.0)
897  , bestCutoff_(0.0)
898  , timeLimit_(0)
899  , startTime_(0)
900  , nodeLimit_(0)
901  , startNode_(-1)
902  , searchType_(-1)
903  , refine_(false)
904{
905}
906/* Constructor with solution.
907   range is upper bound on difference from given solution.
908   maxDiversification is maximum number of diversifications to try
909   timeLimit is seconds in subTree
910   nodeLimit is nodes in subTree
911*/
912CbcTreeVariable::CbcTreeVariable(CbcModel *model, const double *solution,
913  int range, int typeCuts, int maxDiversification,
914  int timeLimit, int nodeLimit, bool refine)
915  : localNode_(NULL)
916  , bestSolution_(NULL)
917  , savedSolution_(NULL)
918  , saveNumberSolutions_(0)
919  , model_(model)
920  , originalLower_(NULL)
921  , originalUpper_(NULL)
922  , range_(range)
923  , typeCuts_(typeCuts)
924  , maxDiversification_(maxDiversification)
925  , diversification_(0)
926  , nextStrong_(false)
927  , rhs_(0.0)
928  , savedGap_(0.0)
929  , bestCutoff_(0.0)
930  , timeLimit_(timeLimit)
931  , startTime_(0)
932  , nodeLimit_(nodeLimit)
933  , startNode_(-1)
934  , searchType_(-1)
935  , refine_(refine)
936{
937
938  OsiSolverInterface *solver = model_->solver();
939  const double *lower = solver->getColLower();
940  const double *upper = solver->getColUpper();
941  //const double * solution = solver->getColSolution();
942  //const double * objective = solver->getObjCoefficients();
943  double primalTolerance;
944  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
945
946  // Get increment
947  model_->analyzeObjective();
948
949  {
950    // needed to sync cutoffs
951    double value;
952    solver->getDblParam(OsiDualObjectiveLimit, value);
953    model_->setCutoff(value * solver->getObjSense());
954  }
955  bestCutoff_ = model_->getCutoff();
956  // save current gap
957  savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
958
959  // make sure integers found
960  model_->findIntegers(false);
961  int numberIntegers = model_->numberIntegers();
962  const int *integerVariable = model_->integerVariable();
963  int i;
964  double direction = solver->getObjSense();
965  double newSolutionValue = 1.0e50;
966  if (solution) {
967    // copy solution
968    solver->setColSolution(solution);
969    newSolutionValue = direction * solver->getObjValue();
970  }
971  originalLower_ = new double[numberIntegers];
972  originalUpper_ = new double[numberIntegers];
973  bool all01 = true;
974  int number01 = 0;
975  for (i = 0; i < numberIntegers; i++) {
976    int iColumn = integerVariable[i];
977    originalLower_[i] = lower[iColumn];
978    originalUpper_[i] = upper[iColumn];
979    if (upper[iColumn] - lower[iColumn] > 1.5)
980      all01 = false;
981    else if (upper[iColumn] - lower[iColumn] == 1.0)
982      number01++;
983  }
984  if (all01 && !typeCuts_)
985    typeCuts_ = 1; // may as well so we don't have to deal with refine
986  if (!number01 && !typeCuts_) {
987    if (model_->messageHandler()->logLevel() > 1)
988      printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
989    typeCuts_ = -1;
990  } else {
991    if (model_->messageHandler()->logLevel() > 1) {
992      std::string type;
993      if (all01) {
994        printf("%d 0-1 variables normal local  cuts\n",
995          number01);
996      } else if (typeCuts_) {
997        printf("%d 0-1 variables, %d other - general integer local cuts\n",
998          number01, numberIntegers - number01);
999      } else {
1000        printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
1001          number01, numberIntegers - number01);
1002      }
1003      printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
1004        maxDiversification_, range_, timeLimit_, nodeLimit_);
1005    }
1006  }
1007  int numberColumns = model_->getNumCols();
1008  savedSolution_ = new double[numberColumns];
1009  memset(savedSolution_, 0, numberColumns * sizeof(double));
1010  if (solution) {
1011    rhs_ = range_;
1012    // Check feasible
1013    int goodSolution = createCut(solution, cut_);
1014    if (goodSolution >= 0) {
1015      for (i = 0; i < numberIntegers; i++) {
1016        int iColumn = integerVariable[i];
1017        double value = floor(solution[iColumn] + 0.5);
1018        // fix so setBestSolution will work
1019        solver->setColLower(iColumn, value);
1020        solver->setColUpper(iColumn, value);
1021      }
1022      model_->reserveCurrentSolution();
1023      // Create cut and get total gap
1024      if (newSolutionValue < bestCutoff_) {
1025        model_->setBestSolution(CBC_ROUNDING, newSolutionValue, solution);
1026        bestCutoff_ = model_->getCutoff();
1027        // save as best solution
1028        memcpy(savedSolution_, model_->bestSolution(), numberColumns * sizeof(double));
1029      }
1030      for (i = 0; i < numberIntegers; i++) {
1031        int iColumn = integerVariable[i];
1032        // restore bounds
1033        solver->setColLower(iColumn, originalLower_[i]);
1034        solver->setColUpper(iColumn, originalUpper_[i]);
1035      }
1036      // make sure can't stop on gap
1037      model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
1038    } else {
1039      model_ = NULL;
1040    }
1041  } else {
1042    // no solution
1043    rhs_ = 1.0e50;
1044    // make sure can't stop on gap
1045    model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
1046  }
1047}
1048CbcTreeVariable::~CbcTreeVariable()
1049{
1050  delete[] originalLower_;
1051  delete[] originalUpper_;
1052  delete[] bestSolution_;
1053  delete[] savedSolution_;
1054  delete localNode_;
1055}
1056// Copy constructor
1057CbcTreeVariable::CbcTreeVariable(const CbcTreeVariable &rhs)
1058  : CbcTree(rhs)
1059  , saveNumberSolutions_(rhs.saveNumberSolutions_)
1060  , model_(rhs.model_)
1061  , range_(rhs.range_)
1062  , typeCuts_(rhs.typeCuts_)
1063  , maxDiversification_(rhs.maxDiversification_)
1064  , diversification_(rhs.diversification_)
1065  , nextStrong_(rhs.nextStrong_)
1066  , rhs_(rhs.rhs_)
1067  , savedGap_(rhs.savedGap_)
1068  , bestCutoff_(rhs.bestCutoff_)
1069  , timeLimit_(rhs.timeLimit_)
1070  , startTime_(rhs.startTime_)
1071  , nodeLimit_(rhs.nodeLimit_)
1072  , startNode_(rhs.startNode_)
1073  , searchType_(rhs.searchType_)
1074  , refine_(rhs.refine_)
1075{
1076  cut_ = rhs.cut_;
1077  fixedCut_ = rhs.fixedCut_;
1078  if (rhs.localNode_)
1079    localNode_ = new CbcNode(*rhs.localNode_);
1080  else
1081    localNode_ = NULL;
1082  if (rhs.originalLower_) {
1083    int numberIntegers = model_->numberIntegers();
1084    originalLower_ = new double[numberIntegers];
1085    memcpy(originalLower_, rhs.originalLower_, numberIntegers * sizeof(double));
1086    originalUpper_ = new double[numberIntegers];
1087    memcpy(originalUpper_, rhs.originalUpper_, numberIntegers * sizeof(double));
1088  } else {
1089    originalLower_ = NULL;
1090    originalUpper_ = NULL;
1091  }
1092  if (rhs.bestSolution_) {
1093    int numberColumns = model_->getNumCols();
1094    bestSolution_ = new double[numberColumns];
1095    memcpy(bestSolution_, rhs.bestSolution_, numberColumns * sizeof(double));
1096  } else {
1097    bestSolution_ = NULL;
1098  }
1099  if (rhs.savedSolution_) {
1100    int numberColumns = model_->getNumCols();
1101    savedSolution_ = new double[numberColumns];
1102    memcpy(savedSolution_, rhs.savedSolution_, numberColumns * sizeof(double));
1103  } else {
1104    savedSolution_ = NULL;
1105  }
1106}
1107//----------------------------------------------------------------
1108// Assignment operator
1109//-------------------------------------------------------------------
1110CbcTreeVariable &
1111CbcTreeVariable::operator=(const CbcTreeVariable &rhs)
1112{
1113  if (this != &rhs) {
1114    CbcTree::operator=(rhs);
1115    saveNumberSolutions_ = rhs.saveNumberSolutions_;
1116    cut_ = rhs.cut_;
1117    fixedCut_ = rhs.fixedCut_;
1118    delete localNode_;
1119    if (rhs.localNode_)
1120      localNode_ = new CbcNode(*rhs.localNode_);
1121    else
1122      localNode_ = NULL;
1123    model_ = rhs.model_;
1124    range_ = rhs.range_;
1125    typeCuts_ = rhs.typeCuts_;
1126    maxDiversification_ = rhs.maxDiversification_;
1127    diversification_ = rhs.diversification_;
1128    nextStrong_ = rhs.nextStrong_;
1129    rhs_ = rhs.rhs_;
1130    savedGap_ = rhs.savedGap_;
1131    bestCutoff_ = rhs.bestCutoff_;
1132    timeLimit_ = rhs.timeLimit_;
1133    startTime_ = rhs.startTime_;
1134    nodeLimit_ = rhs.nodeLimit_;
1135    startNode_ = rhs.startNode_;
1136    searchType_ = rhs.searchType_;
1137    refine_ = rhs.refine_;
1138    delete[] originalLower_;
1139    delete[] originalUpper_;
1140    if (rhs.originalLower_) {
1141      int numberIntegers = model_->numberIntegers();
1142      originalLower_ = new double[numberIntegers];
1143      memcpy(originalLower_, rhs.originalLower_, numberIntegers * sizeof(double));
1144      originalUpper_ = new double[numberIntegers];
1145      memcpy(originalUpper_, rhs.originalUpper_, numberIntegers * sizeof(double));
1146    } else {
1147      originalLower_ = NULL;
1148      originalUpper_ = NULL;
1149    }
1150    delete[] bestSolution_;
1151    if (rhs.bestSolution_) {
1152      int numberColumns = model_->getNumCols();
1153      bestSolution_ = new double[numberColumns];
1154      memcpy(bestSolution_, rhs.bestSolution_, numberColumns * sizeof(double));
1155    } else {
1156      bestSolution_ = NULL;
1157    }
1158    delete[] savedSolution_;
1159    if (rhs.savedSolution_) {
1160      int numberColumns = model_->getNumCols();
1161      savedSolution_ = new double[numberColumns];
1162      memcpy(savedSolution_, rhs.savedSolution_, numberColumns * sizeof(double));
1163    } else {
1164      savedSolution_ = NULL;
1165    }
1166  }
1167  return *this;
1168}
1169// Clone
1170CbcTree *
1171CbcTreeVariable::clone() const
1172{
1173  return new CbcTreeVariable(*this);
1174}
1175// Pass in solution (so can be used after heuristic)
1176void CbcTreeVariable::passInSolution(const double *solution, double solutionValue)
1177{
1178  int numberColumns = model_->getNumCols();
1179  delete[] savedSolution_;
1180  savedSolution_ = new double[numberColumns];
1181  memcpy(savedSolution_, solution, numberColumns * sizeof(double));
1182  rhs_ = range_;
1183  // Check feasible
1184  int goodSolution = createCut(solution, cut_);
1185  if (goodSolution >= 0) {
1186    bestCutoff_ = CoinMin(solutionValue, model_->getCutoff());
1187  } else {
1188    model_ = NULL;
1189  }
1190}
1191// Return the top node of the heap
1192CbcNode *
1193CbcTreeVariable::top() const
1194{
1195#ifdef CBC_DEBUG
1196  int smallest = 9999999;
1197  int largest = -1;
1198  double smallestD = 1.0e30;
1199  double largestD = -1.0e30;
1200  int n = nodes_.size();
1201  for (int i = 0; i < n; i++) {
1202    int nn = nodes_[i]->nodeInfo()->nodeNumber();
1203    double dd = nodes_[i]->objectiveValue();
1204    largest = CoinMax(largest, nn);
1205    smallest = CoinMin(smallest, nn);
1206    largestD = CoinMax(largestD, dd);
1207    smallestD = CoinMin(smallestD, dd);
1208  }
1209  if (model_->messageHandler()->logLevel() > 1) {
1210    printf("smallest %d, largest %d, top %d\n", smallest, largest,
1211      nodes_.front()->nodeInfo()->nodeNumber());
1212    printf("smallestD %g, largestD %g, top %g\n", smallestD, largestD, nodes_.front()->objectiveValue());
1213  }
1214#endif
1215  return nodes_.front();
1216}
1217
1218// Add a node to the heap
1219void CbcTreeVariable::push(CbcNode *x)
1220{
1221  if (typeCuts_ >= 0 && !nodes_.size() && searchType_ < 0) {
1222    startNode_ = model_->getNodeCount();
1223    // save copy of node
1224    localNode_ = new CbcNode(*x);
1225
1226    if (cut_.row().getNumElements()) {
1227      // Add to global cuts
1228      // we came in with solution
1229      model_->makeGlobalCut(cut_);
1230      if (model_->messageHandler()->logLevel() > 1)
1231        printf("initial cut - rhs %g %g\n",
1232          cut_.lb(), cut_.ub());
1233      searchType_ = 1;
1234    } else {
1235      // stop on first solution
1236      searchType_ = 0;
1237    }
1238    startTime_ = static_cast< int >(CoinCpuTime());
1239    saveNumberSolutions_ = model_->getSolutionCount();
1240  }
1241  nodes_.push_back(x);
1242#ifdef CBC_DEBUG
1243  if (model_->messageHandler()->logLevel() > 0)
1244    printf("pushing node onto heap %d %x %x\n",
1245      x->nodeInfo()->nodeNumber(), x, x->nodeInfo());
1246#endif
1247  std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
1248}
1249
1250// Remove the top node from the heap
1251void CbcTreeVariable::pop()
1252{
1253  std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
1254  nodes_.pop_back();
1255}
1256// Test if empty - does work if so
1257bool CbcTreeVariable::empty()
1258{
1259  if (typeCuts_ < 0)
1260    return !nodes_.size();
1261  /* state -
1262       0 iterating
1263       1 subtree finished optimal solution for subtree found
1264       2 subtree finished and no solution found
1265       3 subtree exiting and solution found
1266       4 subtree exiting and no solution found
1267    */
1268  int state = 0;
1269  assert(searchType_ != 2);
1270  if (searchType_) {
1271    if (CoinCpuTime() - startTime_ > timeLimit_ || model_->getNodeCount() - startNode_ >= nodeLimit_) {
1272      state = 4;
1273    }
1274  } else {
1275    if (model_->getSolutionCount() > saveNumberSolutions_) {
1276      state = 4;
1277    }
1278  }
1279  if (!nodes_.size())
1280    state = 2;
1281  if (!state) {
1282    return false;
1283  }
1284  // Finished this phase
1285  int numberColumns = model_->getNumCols();
1286  if (model_->getSolutionCount() > saveNumberSolutions_) {
1287    if (model_->getCutoff() < bestCutoff_) {
1288      // Save solution
1289      if (!bestSolution_)
1290        bestSolution_ = new double[numberColumns];
1291      memcpy(bestSolution_, model_->bestSolution(), numberColumns * sizeof(double));
1292      bestCutoff_ = model_->getCutoff();
1293    }
1294    state--;
1295  }
1296  // get rid of all nodes (safe even if already done)
1297  double bestPossibleObjective;
1298  cleanTree(model_, -COIN_DBL_MAX, bestPossibleObjective);
1299
1300  double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement);
1301  if (model_->messageHandler()->logLevel() > 1)
1302    printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
1303      state,
1304      model_->getNodeCount() - startNode_,
1305      static_cast< int >(CoinCpuTime()) - startTime_,
1306      model_->getCutoff() + increment, bestCutoff_ + increment, rhs_);
1307  saveNumberSolutions_ = model_->getSolutionCount();
1308  bool finished = false;
1309  bool lastTry = false;
1310  switch (state) {
1311  case 1:
1312    // solution found and subtree exhausted
1313    if (rhs_ > 1.0e30) {
1314      finished = true;
1315    } else {
1316      // find global cut and reverse
1317      reverseCut(1);
1318      searchType_ = 1; // first false
1319      rhs_ = range_; // reset range
1320      nextStrong_ = false;
1321
1322      // save best solution in this subtree
1323      memcpy(savedSolution_, model_->bestSolution(), numberColumns * sizeof(double));
1324    }
1325    break;
1326  case 2:
1327    // solution not found and subtree exhausted
1328    if (rhs_ > 1.0e30) {
1329      finished = true;
1330    } else {
1331      // find global cut and reverse
1332      reverseCut(2);
1333      searchType_ = 1; // first false
1334      if (diversification_ < maxDiversification_) {
1335        if (nextStrong_) {
1336          diversification_++;
1337          // cut is valid so don't model_->setCutoff(1.0e50);
1338          searchType_ = 0;
1339        }
1340        nextStrong_ = true;
1341        rhs_ += range_ / 2;
1342      } else {
1343        // This will be last try (may hit max time)
1344        lastTry = true;
1345        if (!maxDiversification_)
1346          typeCuts_ = -1; // make sure can't start again
1347        model_->setCutoff(bestCutoff_);
1348        if (model_->messageHandler()->logLevel() > 1)
1349          printf("Exiting local search with current set of cuts\n");
1350        rhs_ = 1.0e100;
1351        // Can now stop on gap
1352        model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1353      }
1354    }
1355    break;
1356  case 3:
1357    // solution found and subtree not exhausted
1358    if (rhs_ < 1.0e30) {
1359      if (searchType_) {
1360        if (!typeCuts_ && refine_ && searchType_ == 1) {
1361          // We need to check we have best solution given these 0-1 values
1362          OsiSolverInterface *subSolver = model_->continuousSolver()->clone();
1363          CbcModel *subModel = model_->subTreeModel(subSolver);
1364          CbcTree normalTree;
1365          subModel->passInTreeHandler(normalTree);
1366          int numberIntegers = model_->numberIntegers();
1367          const int *integerVariable = model_->integerVariable();
1368          const double *solution = model_->bestSolution();
1369          int i;
1370          int numberColumns = model_->getNumCols();
1371          for (i = 0; i < numberIntegers; i++) {
1372            int iColumn = integerVariable[i];
1373            double value = floor(solution[iColumn] + 0.5);
1374            if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
1375              continue; // skip as not 0-1
1376            if (originalLower_[i] == originalUpper_[i])
1377              continue;
1378            subSolver->setColLower(iColumn, value);
1379            subSolver->setColUpper(iColumn, value);
1380          }
1381          subSolver->initialSolve();
1382          // We can copy cutoff
1383          // But adjust
1384          subModel->setCutoff(model_->getCutoff() + model_->getDblParam(CbcModel::CbcCutoffIncrement) + 1.0e-6);
1385          subModel->setSolutionCount(0);
1386          assert(subModel->isProvenOptimal());
1387          if (!subModel->typePresolve()) {
1388            subModel->branchAndBound();
1389            if (subModel->status()) {
1390              model_->incrementSubTreeStopped();
1391            }
1392            //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
1393            //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
1394            double newCutoff = subModel->getMinimizationObjValue() - subModel->getDblParam(CbcModel::CbcCutoffIncrement);
1395            if (subModel->getSolutionCount()) {
1396              if (!subModel->status())
1397                assert(subModel->isProvenOptimal());
1398              memcpy(model_->bestSolution(), subModel->bestSolution(),
1399                numberColumns * sizeof(double));
1400              model_->setCutoff(newCutoff);
1401            }
1402          } else if (subModel->typePresolve() == 1) {
1403            CbcModel *model2 = subModel->integerPresolve(true);
1404            if (model2) {
1405              // Do complete search
1406              model2->branchAndBound();
1407              // get back solution
1408              subModel->originalModel(model2, false);
1409              if (model2->status()) {
1410                model_->incrementSubTreeStopped();
1411              }
1412              double newCutoff = model2->getMinimizationObjValue() - model2->getDblParam(CbcModel::CbcCutoffIncrement);
1413              if (model2->getSolutionCount()) {
1414                if (!model2->status())
1415                  assert(model2->isProvenOptimal());
1416                memcpy(model_->bestSolution(), subModel->bestSolution(),
1417                  numberColumns * sizeof(double));
1418                model_->setCutoff(newCutoff);
1419              }
1420              delete model2;
1421            } else {
1422              // infeasible - could just be - due to cutoff
1423            }
1424          } else {
1425            // too dangerous at present
1426            assert(subModel->typePresolve() != 2);
1427          }
1428          if (model_->getCutoff() < bestCutoff_) {
1429            // Save solution
1430            if (!bestSolution_)
1431              bestSolution_ = new double[numberColumns];
1432            memcpy(bestSolution_, model_->bestSolution(), numberColumns * sizeof(double));
1433            bestCutoff_ = model_->getCutoff();
1434          }
1435          delete subModel;
1436        }
1437        // we have done search to make sure best general solution
1438        searchType_ = 1;
1439        // Reverse cut weakly
1440        reverseCut(3, rhs_);
1441      } else {
1442        searchType_ = 1;
1443        // delete last cut
1444        deleteCut(cut_);
1445      }
1446    } else {
1447      searchType_ = 1;
1448    }
1449    // save best solution in this subtree
1450    memcpy(savedSolution_, model_->bestSolution(), numberColumns * sizeof(double));
1451    nextStrong_ = false;
1452    rhs_ = range_;
1453    break;
1454  case 4:
1455    // solution not found and subtree not exhausted
1456    if (maxDiversification_) {
1457      if (nextStrong_) {
1458        // Reverse cut weakly
1459        reverseCut(4, rhs_);
1460        model_->setCutoff(1.0e50);
1461        diversification_++;
1462        searchType_ = 0;
1463      } else {
1464        // delete last cut
1465        deleteCut(cut_);
1466        searchType_ = 1;
1467      }
1468      nextStrong_ = true;
1469      rhs_ += range_ / 2;
1470    } else {
1471      // special case when using as heuristic
1472      // Reverse cut weakly if lb -infinity
1473      reverseCut(4, rhs_);
1474      // This will be last try (may hit max time0
1475      lastTry = true;
1476      model_->setCutoff(bestCutoff_);
1477      if (model_->messageHandler()->logLevel() > 1)
1478        printf("Exiting local search with current set of cuts\n");
1479      rhs_ = 1.0e100;
1480      // Can now stop on gap
1481      model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1482      typeCuts_ = -1;
1483    }
1484    break;
1485  }
1486  if (rhs_ < 1.0e30 || lastTry) {
1487    int goodSolution = createCut(savedSolution_, cut_);
1488    if (goodSolution >= 0) {
1489      // Add to global cuts
1490      model_->makeGlobalCut(cut_);
1491      CbcRowCuts *global = model_->globalCuts();
1492      int n = global->sizeRowCuts();
1493      OsiRowCut *rowCut = global->rowCutPtr(n - 1);
1494      if (model_->messageHandler()->logLevel() > 1)
1495        printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
1496          n, rowCut->lb(), rowCut->ub(), rhs_, diversification_);
1497      const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways();
1498      if (debugger) {
1499        if (debugger->invalidCut(*rowCut))
1500          printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1501      }
1502      for (int i = 0; i < n; i++) {
1503        rowCut = global->rowCutPtr(i);
1504        if (model_->messageHandler()->logLevel() > 1)
1505          printf("%d - rhs %g %g\n",
1506            i, rowCut->lb(), rowCut->ub());
1507      }
1508    }
1509    // put back node
1510    startTime_ = static_cast< int >(CoinCpuTime());
1511    startNode_ = model_->getNodeCount();
1512    if (localNode_) {
1513      // save copy of node
1514      CbcNode *localNode2 = new CbcNode(*localNode_);
1515      // But localNode2 now owns cuts so swap
1516      //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
1517      //   localNode_,localNode_->nodeInfo());
1518      nodes_.push_back(localNode_);
1519      localNode_ = localNode2;
1520      std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
1521    }
1522  }
1523  return finished;
1524}
1525// We may have got an intelligent tree so give it one more chance
1526void CbcTreeVariable::endSearch()
1527{
1528  if (typeCuts_ >= 0) {
1529    // copy best solution to model
1530    int numberColumns = model_->getNumCols();
1531    if (bestSolution_ && bestCutoff_ < model_->getCutoff()) {
1532      memcpy(model_->bestSolution(), bestSolution_, numberColumns * sizeof(double));
1533      model_->setCutoff(bestCutoff_);
1534      // recompute objective value
1535      const double *objCoef = model_->getObjCoefficients();
1536      double objOffset = 0.0;
1537      model_->continuousSolver()->getDblParam(OsiObjOffset, objOffset);
1538
1539      // Compute dot product of objCoef and colSol and then adjust by offset
1540      double objValue = -objOffset;
1541      for (int i = 0; i < numberColumns; i++)
1542        objValue += objCoef[i] * bestSolution_[i];
1543      model_->setMinimizationObjValue(objValue);
1544    }
1545    // Can now stop on gap
1546    model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1547  }
1548}
1549// Create cut
1550int CbcTreeVariable::createCut(const double *solution, OsiRowCut &rowCut)
1551{
1552  if (rhs_ > 1.0e20)
1553    return -1;
1554  OsiSolverInterface *solver = model_->solver();
1555  const double *rowLower = solver->getRowLower();
1556  const double *rowUpper = solver->getRowUpper();
1557  //const double * solution = solver->getColSolution();
1558  //const double * objective = solver->getObjCoefficients();
1559  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1560  double primalTolerance;
1561  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1562  // relax
1563  primalTolerance *= 1000.0;
1564
1565  int numberRows = model_->getNumRows();
1566
1567  int numberIntegers = model_->numberIntegers();
1568  const int *integerVariable = model_->integerVariable();
1569  int i;
1570
1571  // Check feasible
1572
1573  double *rowActivity = new double[numberRows];
1574  memset(rowActivity, 0, numberRows * sizeof(double));
1575  solver->getMatrixByCol()->times(solution, rowActivity);
1576  int goodSolution = 0;
1577  // check was feasible
1578  for (i = 0; i < numberRows; i++) {
1579    if (rowActivity[i] < rowLower[i] - primalTolerance) {
1580      goodSolution = -1;
1581    } else if (rowActivity[i] > rowUpper[i] + primalTolerance) {
1582      goodSolution = -1;
1583    }
1584  }
1585  delete[] rowActivity;
1586  for (i = 0; i < numberIntegers; i++) {
1587    int iColumn = integerVariable[i];
1588    double value = solution[iColumn];
1589    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
1590      goodSolution = -1;
1591    }
1592  }
1593  // zap cut
1594  if (goodSolution == 0) {
1595    // Create cut and get total gap
1596    CoinPackedVector cut;
1597    double rhs = rhs_;
1598    double maxValue = 0.0;
1599    for (i = 0; i < numberIntegers; i++) {
1600      int iColumn = integerVariable[i];
1601      double value = floor(solution[iColumn] + 0.5);
1602      if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
1603        continue; // skip as not 0-1
1604      if (originalLower_[i] == originalUpper_[i])
1605        continue;
1606      double mu = 1.0 / (originalUpper_[i] - originalLower_[i]);
1607      if (value == originalLower_[i]) {
1608        rhs += mu * originalLower_[i];
1609        cut.insert(iColumn, 1.0);
1610        maxValue += originalUpper_[i];
1611      } else if (value == originalUpper_[i]) {
1612        rhs -= mu * originalUpper_[i];
1613        cut.insert(iColumn, -1.0);
1614        maxValue += originalLower_[i];
1615      }
1616    }
1617    if (maxValue < rhs - primalTolerance) {
1618      if (model_->messageHandler()->logLevel() > 1)
1619        printf("slack cut\n");
1620      goodSolution = 1;
1621    }
1622    rowCut.setRow(cut);
1623    rowCut.setLb(-COIN_DBL_MAX);
1624    rowCut.setUb(rhs);
1625    rowCut.setGloballyValid();
1626    if (model_->messageHandler()->logLevel() > 1)
1627      printf("Cut size: %i Cut rhs: %g\n", cut.getNumElements(), rhs);
1628#ifdef CBC_DEBUG
1629    if (model_->messageHandler()->logLevel() > 0) {
1630      int k;
1631      for (k = 0; k < cut.getNumElements(); k++) {
1632        printf("%i    %g ", cut.getIndices()[k], cut.getElements()[k]);
1633        if ((k + 1) % 5 == 0)
1634          printf("\n");
1635      }
1636      if (k % 5 != 0)
1637        printf("\n");
1638    }
1639#endif
1640    return goodSolution;
1641  } else {
1642    if (model_->messageHandler()->logLevel() > 1)
1643      printf("Not a good solution\n");
1644    return -1;
1645  }
1646}
1647// Other side of last cut branch
1648void CbcTreeVariable::reverseCut(int state, double bias)
1649{
1650  // find global cut
1651  CbcRowCuts *global = model_->globalCuts();
1652  int n = global->sizeRowCuts();
1653  int i;
1654  OsiRowCut *rowCut = NULL;
1655  for (i = 0; i < n; i++) {
1656    rowCut = global->rowCutPtr(i);
1657    if (cut_ == *rowCut) {
1658      break;
1659    }
1660  }
1661  if (!rowCut) {
1662    // must have got here in odd way e.g. strong branching
1663    return;
1664  }
1665  if (rowCut->lb() > -1.0e10)
1666    return;
1667  // get smallest element
1668  double smallest = COIN_DBL_MAX;
1669  CoinPackedVector row = cut_.row();
1670  for (int k = 0; k < row.getNumElements(); k++)
1671    smallest = CoinMin(smallest, fabs(row.getElements()[k]));
1672  if (!typeCuts_ && !refine_) {
1673    // Reverse cut very very weakly
1674    if (state > 2)
1675      smallest = 0.0;
1676  }
1677  // replace by other way
1678  if (model_->messageHandler()->logLevel() > 1)
1679    printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
1680      i, n, rowCut->lb(), rowCut->ub());
1681  rowCut->setLb(rowCut->ub() + smallest - bias);
1682  rowCut->setUb(COIN_DBL_MAX);
1683  if (model_->messageHandler()->logLevel() > 1)
1684    printf("new rhs %g %g, bias %g smallest %g ",
1685      rowCut->lb(), rowCut->ub(), bias, smallest);
1686  const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways();
1687  if (debugger) {
1688    if (debugger->invalidCut(*rowCut))
1689      printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1690  }
1691}
1692// Delete last cut branch
1693void CbcTreeVariable::deleteCut(OsiRowCut &cut)
1694{
1695  // find global cut
1696  CbcRowCuts *global = model_->globalCuts();
1697  int n = global->sizeRowCuts();
1698  int i;
1699  OsiRowCut *rowCut = NULL;
1700  for (i = 0; i < n; i++) {
1701    rowCut = global->rowCutPtr(i);
1702    if (cut == *rowCut) {
1703      break;
1704    }
1705  }
1706  assert(i < n);
1707  // delete last cut
1708  if (model_->messageHandler()->logLevel() > 1)
1709    printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
1710      i, n, rowCut->lb(), rowCut->ub());
1711  global->eraseRowCut(i);
1712}
1713// Create C++ lines to get to current state
1714void CbcTreeVariable::generateCpp(FILE *fp)
1715{
1716  CbcTreeVariable other;
1717  fprintf(fp, "0#include \"CbcTreeVariable.hpp\"\n");
1718  fprintf(fp, "5  CbcTreeVariable variableTree(cbcModel,NULL);\n");
1719  if (range_ != other.range_)
1720    fprintf(fp, "5  variableTree.setRange(%d);\n", range_);
1721  if (typeCuts_ != other.typeCuts_)
1722    fprintf(fp, "5  variableTree.setTypeCuts(%d);\n", typeCuts_);
1723  if (maxDiversification_ != other.maxDiversification_)
1724    fprintf(fp, "5  variableTree.setMaxDiversification(%d);\n", maxDiversification_);
1725  if (timeLimit_ != other.timeLimit_)
1726    fprintf(fp, "5  variableTree.setTimeLimit(%d);\n", timeLimit_);
1727  if (nodeLimit_ != other.nodeLimit_)
1728    fprintf(fp, "5  variableTree.setNodeLimit(%d);\n", nodeLimit_);
1729  if (refine_ != other.refine_)
1730    fprintf(fp, "5  variableTree.setRefine(%s);\n", refine_ ? "true" : "false");
1731  fprintf(fp, "5  cbcModel->passInTreeHandler(variableTree);\n");
1732}
1733
1734/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1735*/
Note: See TracBrowser for help on using the repository browser.