source: trunk/Cbc/src/CbcHeuristicLocal.cpp

Last change on this file was 2631, checked in by unxusr, 3 months ago

memory leak reported by coverity

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 53.8 KB
Line 
1/* $Id: CbcHeuristicLocal.cpp 2631 2019-07-17 10:28:11Z stefan $ */
2// Copyright (C) 2002, 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#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#pragma warning(disable : 4786)
9#endif
10#include <cassert>
11#include <cstdlib>
12#include <cmath>
13#include <cfloat>
14
15#include "OsiSolverInterface.hpp"
16#include "CbcModel.hpp"
17#include "CbcMessage.hpp"
18#include "CbcHeuristicLocal.hpp"
19#include "CbcHeuristicFPump.hpp"
20#include "CbcBranchActual.hpp"
21#include "CbcStrategy.hpp"
22#include "CglPreProcess.hpp"
23
24// Default Constructor
25CbcHeuristicLocal::CbcHeuristicLocal()
26  : CbcHeuristic()
27{
28  numberSolutions_ = 0;
29  swap_ = 0;
30  used_ = NULL;
31  lastRunDeep_ = -1000000;
32  switches_ |= 16; // needs a new solution
33}
34
35// Constructor with model - assumed before cuts
36
37CbcHeuristicLocal::CbcHeuristicLocal(CbcModel &model)
38  : CbcHeuristic(model)
39{
40  numberSolutions_ = 0;
41  swap_ = 0;
42  lastRunDeep_ = -1000000;
43  switches_ |= 16; // needs a new solution
44  // Get a copy of original matrix
45  assert(model.solver());
46  if (model.solver()->getNumRows()) {
47    matrix_ = *model.solver()->getMatrixByCol();
48  }
49  int numberColumns = model.solver()->getNumCols();
50  used_ = new int[numberColumns];
51  memset(used_, 0, numberColumns * sizeof(int));
52}
53
54// Destructor
55CbcHeuristicLocal::~CbcHeuristicLocal()
56{
57  delete[] used_;
58}
59
60// Clone
61CbcHeuristic *
62CbcHeuristicLocal::clone() const
63{
64  return new CbcHeuristicLocal(*this);
65}
66// Create C++ lines to get to current state
67void CbcHeuristicLocal::generateCpp(FILE *fp)
68{
69  CbcHeuristicLocal other;
70  fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n");
71  fprintf(fp, "3  CbcHeuristicLocal heuristicLocal(*cbcModel);\n");
72  CbcHeuristic::generateCpp(fp, "heuristicLocal");
73  if (swap_ != other.swap_)
74    fprintf(fp, "3  heuristicLocal.setSearchType(%d);\n", swap_);
75  else
76    fprintf(fp, "4  heuristicLocal.setSearchType(%d);\n", swap_);
77  fprintf(fp, "3  cbcModel->addHeuristic(&heuristicLocal);\n");
78}
79
80// Copy constructor
81CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal &rhs)
82  : CbcHeuristic(rhs)
83  , matrix_(rhs.matrix_)
84  , numberSolutions_(rhs.numberSolutions_)
85  , swap_(rhs.swap_)
86{
87  if (model_ && rhs.used_) {
88    int numberColumns = model_->solver()->getNumCols();
89    used_ = CoinCopyOfArray(rhs.used_, numberColumns);
90  } else {
91    used_ = NULL;
92  }
93}
94
95// Assignment operator
96CbcHeuristicLocal &
97CbcHeuristicLocal::operator=(const CbcHeuristicLocal &rhs)
98{
99  if (this != &rhs) {
100    CbcHeuristic::operator=(rhs);
101    matrix_ = rhs.matrix_;
102    numberSolutions_ = rhs.numberSolutions_;
103    swap_ = rhs.swap_;
104    delete[] used_;
105    if (model_ && rhs.used_) {
106      int numberColumns = model_->solver()->getNumCols();
107      used_ = CoinCopyOfArray(rhs.used_, numberColumns);
108    } else {
109      used_ = NULL;
110    }
111  }
112  return *this;
113}
114
115// Resets stuff if model changes
116void CbcHeuristicLocal::resetModel(CbcModel * /*model*/)
117{
118  //CbcHeuristic::resetModel(model);
119  delete[] used_;
120  if (model_ && used_) {
121    int numberColumns = model_->solver()->getNumCols();
122    used_ = new int[numberColumns];
123    memset(used_, 0, numberColumns * sizeof(int));
124  } else {
125    used_ = NULL;
126  }
127}
128/*
129  Run a mini-BaB search after fixing all variables not marked as used by
130  solution(). (See comments there for semantics.)
131
132  Return values are:
133    1: smallBranchAndBound found a solution
134    0: everything else
135
136  The degree of overload as return codes from smallBranchAndBound are folded
137  into 0 is such that it's impossible to distinguish return codes that really
138  require attention from a simple `nothing of interest'.
139*/
140// This version fixes stuff and does IP
141int CbcHeuristicLocal::solutionFix(double &objectiveValue,
142  double *newSolution,
143  const int * /*keep*/)
144{
145  /*
146  If when is set to off (0), or set to root (1) and we're not at the root,
147  return. If this heuristic discovered the current solution, don't continue.
148*/
149
150  numCouldRun_++;
151  // See if to do
152  if (!when() || (when() == 1 && model_->phase() != 1))
153    return 0; // switched off
154  // Don't do if it was this heuristic which found solution!
155  if (this == model_->lastHeuristic())
156    return 0;
157  /*
158  Load up a new solver with the solution.
159
160  Why continuousSolver(), as opposed to solver()?
161*/
162  OsiSolverInterface *newSolver = model_->continuousSolver()->clone();
163  const double *colLower = newSolver->getColLower();
164  //const double * colUpper = newSolver->getColUpper();
165
166  int numberIntegers = model_->numberIntegers();
167  const int *integerVariable = model_->integerVariable();
168  /*
169  The net effect here is that anything that hasn't moved from its lower bound
170  will be fixed at lower bound.
171
172  See comments in solution() w.r.t. asymmetric treatment of upper and lower
173  bounds.
174*/
175
176  int i;
177  int nFix = 0;
178  for (i = 0; i < numberIntegers; i++) {
179    int iColumn = integerVariable[i];
180    if (!isHeuristicInteger(newSolver, iColumn))
181      continue;
182    const OsiObject *object = model_->object(i);
183    // get original bounds
184    double originalLower;
185    double originalUpper;
186    getIntegerInformation(object, originalLower, originalUpper);
187    newSolver->setColLower(iColumn, CoinMax(colLower[iColumn], originalLower));
188    if (!used_[iColumn]) {
189      newSolver->setColUpper(iColumn, colLower[iColumn]);
190      nFix++;
191    }
192  }
193  /*
194  Try a `small' branch-and-bound search. The notion here is that we've fixed a
195  lot of variables and reduced the amount of `free' problem to a point where a
196  small BaB search will suffice to fully explore the remaining problem. This
197  routine will execute integer presolve, then call branchAndBound to do the
198  actual search.
199*/
200  int returnCode = 0;
201#ifdef CLP_INVESTIGATE2
202  printf("Fixing %d out of %d (%d continuous)\n",
203    nFix, numberIntegers, newSolver->getNumCols() - numberIntegers);
204#endif
205  if (nFix * 10 <= numberIntegers) {
206    // see if we can fix more
207    int *which = new int[2 * (numberIntegers - nFix)];
208    int *sort = which + (numberIntegers - nFix);
209    int n = 0;
210    for (i = 0; i < numberIntegers; i++) {
211      int iColumn = integerVariable[i];
212      if (!isHeuristicInteger(newSolver, iColumn))
213        continue;
214      if (used_[iColumn]) {
215        which[n] = iColumn;
216        sort[n++] = used_[iColumn];
217      }
218    }
219    CoinSort_2(sort, sort + n, which);
220    // only half fixed in total
221    n = CoinMin(n, numberIntegers / 2 - nFix);
222    int allow = CoinMax(numberSolutions_ - 2, sort[0]);
223    int nFix2 = 0;
224    for (i = 0; i < n; i++) {
225      int iColumn = integerVariable[i];
226      if (!isHeuristicInteger(newSolver, iColumn))
227        continue;
228      if (used_[iColumn] <= allow) {
229        newSolver->setColUpper(iColumn, colLower[iColumn]);
230        nFix2++;
231      } else {
232        break;
233      }
234    }
235    delete[] which;
236    nFix += nFix2;
237#ifdef CLP_INVESTIGATE2
238    printf("Number fixed increased from %d to %d\n",
239      nFix - nFix2, nFix);
240#endif
241  }
242  if (nFix * 10 > numberIntegers) {
243    returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, objectiveValue,
244      objectiveValue, "CbcHeuristicLocal");
245    /*
246  -2 is return due to user event, and -1 is overloaded with what look to be
247  two contradictory meanings.
248*/
249    if (returnCode < 0) {
250      returnCode = 0; // returned on size
251      int numberColumns = newSolver->getNumCols();
252      int numberContinuous = numberColumns - numberIntegers;
253      if (numberContinuous > 2 * numberIntegers && nFix * 10 < numberColumns) {
254#define LOCAL_FIX_CONTINUOUS
255#ifdef LOCAL_FIX_CONTINUOUS
256        //const double * colUpper = newSolver->getColUpper();
257        const double *colLower = newSolver->getColLower();
258        int nAtLb = 0;
259        //double sumDj=0.0;
260        const double *dj = newSolver->getReducedCost();
261        double direction = newSolver->getObjSense();
262        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
263          if (!isHeuristicInteger(newSolver, iColumn)) {
264            if (!used_[iColumn]) {
265              //double djValue = dj[iColumn]*direction;
266              nAtLb++;
267              //sumDj += djValue;
268            }
269          }
270        }
271        if (nAtLb) {
272          // fix some continuous
273          double *sort = new double[nAtLb];
274          int *which = new int[nAtLb];
275          //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
276          int nFix2 = 0;
277          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
278            if (!isHeuristicInteger(newSolver, iColumn)) {
279              if (!used_[iColumn]) {
280                double djValue = dj[iColumn] * direction;
281                if (djValue > 1.0e-6) {
282                  sort[nFix2] = -djValue;
283                  which[nFix2++] = iColumn;
284                }
285              }
286            }
287          }
288          CoinSort_2(sort, sort + nFix2, which);
289          int divisor = 2;
290          nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
291          for (int i = 0; i < nFix2; i++) {
292            int iColumn = which[i];
293            newSolver->setColUpper(iColumn, colLower[iColumn]);
294          }
295          delete[] sort;
296          delete[] which;
297#ifdef CLP_INVESTIGATE2
298          printf("%d integers have zero value, and %d continuous fixed at lb\n",
299            nFix, nFix2);
300#endif
301          returnCode = smallBranchAndBound(newSolver,
302            numberNodes_, newSolution,
303            objectiveValue,
304            objectiveValue, "CbcHeuristicLocal");
305          if (returnCode < 0)
306            returnCode = 0; // returned on size
307        }
308#endif
309      }
310    }
311  }
312  /*
313  If the result is complete exploration with a solution (3) or proven
314  infeasibility (2), we could generate a cut (the AI folks would call it a
315  nogood) to prevent us from going down this route in the future.
316*/
317  if ((returnCode & 2) != 0) {
318    // could add cut
319    returnCode &= ~2;
320  }
321
322  delete newSolver;
323  return returnCode;
324}
325/*
326  First tries setting a variable to better value.  If feasible then
327  tries setting others.  If not feasible then tries swaps
328  Returns 1 if solution, 0 if not
329  The main body of this routine implements an O((q^2)/2) brute force search
330  around the current solution, for q = number of integer variables. Call this
331  the inc/dec heuristic.  For each integer variable x<i>, first decrement the
332  value. Then, for integer variables x<i+1>, ..., x<q-1>, try increment and
333  decrement. If one of these permutations produces a better solution,
334  remember it.  Then repeat, with x<i> incremented. If we find a better
335  solution, update our notion of current solution and continue.
336
337  The net effect is a greedy walk: As each improving pair is found, the
338  current solution is updated and the search continues from this updated
339  solution.
340
341  Way down at the end, we call solutionFix, which will create a drastically
342  restricted problem based on variables marked as used, then do mini-BaC on
343  the restricted problem. This can occur even if we don't try the inc/dec
344  heuristic. This would be more obvious if the inc/dec heuristic were broken
345  out as a separate routine and solutionFix had a name that reflected where
346  it was headed.
347
348  The return code of 0 is grossly overloaded, because it maps to a return
349  code of 0 from solutionFix, which is itself grossly overloaded. See
350  comments in solutionFix and in CbcHeuristic::smallBranchAndBound.
351  */
352int CbcHeuristicLocal::solution(double &solutionValue,
353  double *betterSolution)
354{
355  /*
356  Execute only if a new solution has been discovered since the last time we
357  were called.
358*/
359
360  numCouldRun_++;
361  // See if frequency kills off idea
362  int swap = swap_ % 100;
363  int skip = swap_ / 100;
364  int nodeCount = model_->getNodeCount();
365  if (nodeCount < lastRunDeep_ + skip && nodeCount != lastRunDeep_ + 1)
366    return 0;
367  if (numberSolutions_ == model_->getSolutionCount() && (numberSolutions_ == howOftenShallow_ || nodeCount < lastRunDeep_ + 2 * skip))
368    return 0;
369  howOftenShallow_ = numberSolutions_;
370  numberSolutions_ = model_->getSolutionCount();
371  if (nodeCount < lastRunDeep_ + skip)
372    return 0;
373#ifdef HEURISTIC_INFORM
374  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
375    heuristicName(), numRuns_, numCouldRun_, when_);
376#endif
377  lastRunDeep_ = nodeCount;
378  howOftenShallow_ = numberSolutions_;
379
380  if ((swap % 10) == 2) {
381    // try merge
382    return solutionFix(solutionValue, betterSolution, NULL);
383  }
384  /*
385  Exclude long (column), thin (row) systems.
386
387  Given the n^2 nature of the search, more than 100,000 columns could get
388  expensive. But I don't yet see the rationale for the second part of the
389  condition (cols > 10*rows). And cost is proportional to number of integer
390  variables --- shouldn't we use that?
391
392  Why wait until we have more than one solution?
393*/
394  if ((model_->getNumCols() > 100000 && model_->getNumCols() > 10 * model_->getNumRows()) || numberSolutions_ <= 1)
395    return 0; // probably not worth it
396  // worth trying
397
398  OsiSolverInterface *solver = model_->solver();
399  const double *rowLower = solver->getRowLower();
400  const double *rowUpper = solver->getRowUpper();
401  const double *solution = model_->bestSolution();
402  /*
403  Shouldn't this test be redundant if we've already checked that
404  numberSolutions_ > 1? Stronger: shouldn't this be an assertion?
405*/
406  if (!solution)
407    return 0; // No solution found yet
408  const double *objective = solver->getObjCoefficients();
409  double primalTolerance;
410  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
411
412  int numberRows = matrix_.getNumRows();
413
414  int numberIntegers = model_->numberIntegers();
415  const int *integerVariable = model_->integerVariable();
416
417  int i;
418  double direction = solver->getObjSense();
419  double newSolutionValue = model_->getObjValue() * direction;
420  int returnCode = 0;
421  numRuns_++;
422  // Column copy
423  const double *element = matrix_.getElements();
424  const int *row = matrix_.getIndices();
425  const CoinBigIndex *columnStart = matrix_.getVectorStarts();
426  const int *columnLength = matrix_.getVectorLengths();
427
428  // Get solution array for heuristic solution
429  int numberColumns = solver->getNumCols();
430  double *newSolution = new double[numberColumns];
431  memcpy(newSolution, solution, numberColumns * sizeof(double));
432#ifdef LOCAL_FIX_CONTINUOUS
433  // mark continuous used
434  const double *columnLower = solver->getColLower();
435  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
436    if (!isHeuristicInteger(solver, iColumn)) {
437      if (solution[iColumn] > columnLower[iColumn] + 1.0e-8)
438        used_[iColumn] = numberSolutions_;
439    }
440  }
441#endif
442
443  // way is 1 if down possible, 2 if up possible, 3 if both possible
444  char *way = new char[numberIntegers];
445  // corrected costs
446  double *cost = new double[numberIntegers];
447  // for array to mark infeasible rows after iColumn branch
448  char *mark = new char[numberRows];
449  memset(mark, 0, numberRows);
450  // space to save values so we don't introduce rounding errors
451  double *save = new double[numberRows];
452  /*
453  Force variables within their original bounds, then to the nearest integer.
454  Overall, we seem to be prepared to cope with noninteger bounds. Is this
455  necessary? Seems like we'd be better off to force the bounds to integrality
456  as part of preprocessing.  More generally, why do we need to do this? This
457  solution should have been cleaned and checked when it was accepted as a
458  solution!
459
460  Once the value is set, decide whether we can move up or down.
461
462  The only place that used_ is used is in solutionFix; if a variable is not
463  flagged as used, it will be fixed (at lower bound). Why the asymmetric
464  treatment? This makes some sense for binary variables (for which there are
465  only two options). But for general integer variables, why not make a similar
466  test against the original upper bound?
467*/
468
469  // clean solution
470  for (i = 0; i < numberIntegers; i++) {
471    int iColumn = integerVariable[i];
472    if (!isHeuristicInteger(solver, iColumn))
473      continue;
474    const OsiObject *object = model_->object(i);
475    // get original bounds
476    double originalLower;
477    double originalUpper;
478    getIntegerInformation(object, originalLower, originalUpper);
479    double value = newSolution[iColumn];
480    if (value < originalLower) {
481      value = originalLower;
482      newSolution[iColumn] = value;
483    } else if (value > originalUpper) {
484      value = originalUpper;
485      newSolution[iColumn] = value;
486    }
487    double nearest = floor(value + 0.5);
488    //assert(fabs(value-nearest)<10.0*primalTolerance);
489    value = nearest;
490    newSolution[iColumn] = nearest;
491    // if away from lower bound mark that fact
492    if (nearest > originalLower) {
493      used_[iColumn] = numberSolutions_;
494    }
495    cost[i] = direction * objective[iColumn];
496    /*
497  Given previous computation we're checking that value is at least 1 away
498  from the original bounds.
499*/
500    int iway = 0;
501
502    if (value > originalLower + 0.5)
503      iway = 1;
504    if (value < originalUpper - 0.5)
505      iway |= 2;
506    way[i] = static_cast< char >(iway);
507  }
508  /*
509  Calculate lhs of each constraint for groomed solution.
510*/
511  // get row activities
512  double *rowActivity = new double[numberRows];
513  memset(rowActivity, 0, numberRows * sizeof(double));
514
515  for (i = 0; i < numberColumns; i++) {
516    CoinBigIndex j;
517    double value = newSolution[i];
518    if (value) {
519      for (j = columnStart[i];
520           j < columnStart[i] + columnLength[i]; j++) {
521        int iRow = row[j];
522        rowActivity[iRow] += value * element[j];
523      }
524    }
525  }
526  /*
527  Check that constraints are satisfied. For small infeasibility, force the
528  activity within bound. Again, why is this necessary if the current solution
529  was accepted as a valid solution?
530
531  Why are we scanning past the first unacceptable constraint?
532*/
533  // check was feasible - if not adjust (cleaning may move)
534  // if very infeasible then give up
535  bool tryHeuristic = true;
536  for (i = 0; i < numberRows; i++) {
537    if (rowActivity[i] < rowLower[i]) {
538      if (rowActivity[i] < rowLower[i] - 10.0 * primalTolerance)
539        tryHeuristic = false;
540      rowActivity[i] = rowLower[i];
541    } else if (rowActivity[i] > rowUpper[i]) {
542      if (rowActivity[i] < rowUpper[i] + 10.0 * primalTolerance)
543        tryHeuristic = false;
544      rowActivity[i] = rowUpper[i];
545    }
546  }
547  /*
548  This bit of code is not quite totally redundant: it'll bail at 10,000
549  instead of 100,000. Potentially we can do a lot of work to get here, only
550  to abandon it.
551*/
552  // Switch off if may take too long
553  if (model_->getNumCols() > 10000 && model_->getNumCols() > 10 * model_->getNumRows() && swap < 10)
554    tryHeuristic = false;
555  /*
556  Try the inc/dec heuristic?
557*/
558  if (tryHeuristic) {
559
560    // total change in objective
561    double totalChange = 0.0;
562    // local best change in objective
563    double bestChange = 0.0;
564    // maybe just do 1000
565    int maxIntegers = numberIntegers;
566    // stop if too many goes
567    int maxTries = COIN_INT_MAX;
568    // integerVariable may be randomized copy!
569    int *integerVariable = CoinCopyOfArray(model_->integerVariable(), numberIntegers);
570#ifdef COIN_HAS_CLP
571    OsiClpSolverInterface *clpSolver
572      = dynamic_cast< OsiClpSolverInterface * >(model_->solver());
573    if (clpSolver) {
574      // take out some integers
575      int nn = numberIntegers;
576      numberIntegers = 0;
577      for (int i = 0; i < nn; i++) {
578        int iColumn = integerVariable[i];
579        if (clpSolver->isHeuristicInteger(iColumn))
580          integerVariable[numberIntegers++] = iColumn;
581      }
582    }
583#endif
584    if (swap > 9 && numberIntegers > 500) {
585      int type = swap / 10;
586      if (type == 1) {
587        // reduce
588        maxIntegers = CoinMin(1000, numberIntegers);
589      } else if (type == 2) {
590        // reduce even more
591        maxTries = 100000;
592        maxIntegers = CoinMin(500, numberIntegers);
593      } else if (type > 2) {
594        assert(type < 10);
595        int totals[7] = { 1000, 500, 100, 50, 50, 50, 50 };
596        maxIntegers = CoinMin(totals[type - 3], numberIntegers);
597        double *weight = new double[numberIntegers];
598        for (int i = 0; i < numberIntegers; i++) {
599          weight[i] = model_->randomNumberGenerator()->randomDouble();
600        }
601        CoinSort_2(weight, weight + numberIntegers, integerVariable);
602        delete[] weight;
603      }
604    }
605    /*
606  Outer loop to walk integer variables. Call the current variable x<i>. At the
607  end of this loop, bestChange will contain the best (negative) change in the
608  objective for any single pair.
609
610  The trouble is, we're limited to monotonically increasing improvement.
611  Suppose we discover an improvement of 10 for some pair. If, later in the
612  search, we discover an improvement of 9 for some other pair, we will not use
613  it. That seems wasteful.
614*/
615
616    for (i = 0; i < numberIntegers; i++) {
617      int iColumn = integerVariable[i];
618      bestChange = 0.0;
619      int endInner = CoinMin(numberIntegers, i + maxIntegers);
620
621      double objectiveCoefficient = cost[i];
622      int k;
623      CoinBigIndex j;
624      int goodK = -1;
625      int wayK = -1, wayI = -1;
626      /*
627  Try decrementing x<i>.
628*/
629      if ((way[i] & 1) != 0) {
630        int numberInfeasible = 0;
631        /*
632  Adjust row activities where x<i> has a nonzero coefficient. Save the old
633  values for restoration. Mark any rows that become infeasible as a result
634  of the decrement.
635*/
636        // save row activities and adjust
637        for (j = columnStart[iColumn];
638             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
639          int iRow = row[j];
640          save[iRow] = rowActivity[iRow];
641          rowActivity[iRow] -= element[j];
642          if (rowActivity[iRow] < rowLower[iRow] - primalTolerance || rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
643            // mark row
644            mark[iRow] = 1;
645            numberInfeasible++;
646          }
647        }
648        /*
649  Run through the remaining integer variables. Try increment and decrement on
650  each one. If the potential objective change is better than anything we've
651  seen so far, do a full evaluation of x<k> in that direction.  If we can
652  repair all infeasibilities introduced by pushing x<i> down, we have a
653  winner. Remember the best variable, and the direction for x<i> and x<k>.
654*/
655        // try down
656        for (k = i + 1; k < endInner; k++) {
657          if (!maxTries)
658            break;
659          maxTries--;
660          if ((way[k] & 1) != 0) {
661            // try down
662            if (-objectiveCoefficient - cost[k] < bestChange) {
663              // see if feasible down
664              bool good = true;
665              int numberMarked = 0;
666              int kColumn = integerVariable[k];
667              for (j = columnStart[kColumn];
668                   j < columnStart[kColumn] + columnLength[kColumn]; j++) {
669                int iRow = row[j];
670                double newValue = rowActivity[iRow] - element[j];
671                if (newValue < rowLower[iRow] - primalTolerance || newValue > rowUpper[iRow] + primalTolerance) {
672                  good = false;
673                  break;
674                } else if (mark[iRow]) {
675                  // made feasible
676                  numberMarked++;
677                }
678              }
679              if (good && numberMarked == numberInfeasible) {
680                // better solution
681                goodK = k;
682                wayK = -1;
683                wayI = -1;
684                bestChange = -objectiveCoefficient - cost[k];
685              }
686            }
687          }
688          if ((way[k] & 2) != 0) {
689            // try up
690            if (-objectiveCoefficient + cost[k] < bestChange) {
691              // see if feasible up
692              bool good = true;
693              int numberMarked = 0;
694              int kColumn = integerVariable[k];
695              for (j = columnStart[kColumn];
696                   j < columnStart[kColumn] + columnLength[kColumn]; j++) {
697                int iRow = row[j];
698                double newValue = rowActivity[iRow] + element[j];
699                if (newValue < rowLower[iRow] - primalTolerance || newValue > rowUpper[iRow] + primalTolerance) {
700                  good = false;
701                  break;
702                } else if (mark[iRow]) {
703                  // made feasible
704                  numberMarked++;
705                }
706              }
707              if (good && numberMarked == numberInfeasible) {
708                // better solution
709                goodK = k;
710                wayK = 1;
711                wayI = -1;
712                bestChange = -objectiveCoefficient + cost[k];
713              }
714            }
715          }
716        }
717        /*
718  Remove effect of decrementing x<i> by restoring original lhs values.
719*/
720        // restore row activities
721        for (j = columnStart[iColumn];
722             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
723          int iRow = row[j];
724          rowActivity[iRow] = save[iRow];
725          mark[iRow] = 0;
726        }
727      }
728      /*
729  Try to increment x<i>. Actions as for decrement.
730*/
731      if ((way[i] & 2) != 0) {
732        int numberInfeasible = 0;
733        // save row activities and adjust
734        for (j = columnStart[iColumn];
735             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
736          int iRow = row[j];
737          save[iRow] = rowActivity[iRow];
738          rowActivity[iRow] += element[j];
739          if (rowActivity[iRow] < rowLower[iRow] - primalTolerance || rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
740            // mark row
741            mark[iRow] = 1;
742            numberInfeasible++;
743          }
744        }
745        // try up
746        for (k = i + 1; k < endInner; k++) {
747          if (!maxTries)
748            break;
749          if ((way[k] & 1) != 0) {
750            // try down
751            if (objectiveCoefficient - cost[k] < bestChange) {
752              // see if feasible down
753              bool good = true;
754              int numberMarked = 0;
755              int kColumn = integerVariable[k];
756              for (j = columnStart[kColumn];
757                   j < columnStart[kColumn] + columnLength[kColumn]; j++) {
758                int iRow = row[j];
759                double newValue = rowActivity[iRow] - element[j];
760                if (newValue < rowLower[iRow] - primalTolerance || newValue > rowUpper[iRow] + primalTolerance) {
761                  good = false;
762                  break;
763                } else if (mark[iRow]) {
764                  // made feasible
765                  numberMarked++;
766                }
767              }
768              if (good && numberMarked == numberInfeasible) {
769                // better solution
770                goodK = k;
771                wayK = -1;
772                wayI = 1;
773                bestChange = objectiveCoefficient - cost[k];
774              }
775            }
776          }
777          if ((way[k] & 2) != 0) {
778            // try up
779            if (objectiveCoefficient + cost[k] < bestChange) {
780              // see if feasible up
781              bool good = true;
782              int numberMarked = 0;
783              int kColumn = integerVariable[k];
784              for (j = columnStart[kColumn];
785                   j < columnStart[kColumn] + columnLength[kColumn]; j++) {
786                int iRow = row[j];
787                double newValue = rowActivity[iRow] + element[j];
788                if (newValue < rowLower[iRow] - primalTolerance || newValue > rowUpper[iRow] + primalTolerance) {
789                  good = false;
790                  break;
791                } else if (mark[iRow]) {
792                  // made feasible
793                  numberMarked++;
794                }
795              }
796              if (good && numberMarked == numberInfeasible) {
797                // better solution
798                goodK = k;
799                wayK = 1;
800                wayI = 1;
801                bestChange = objectiveCoefficient + cost[k];
802              }
803            }
804          }
805        }
806        // restore row activities
807        for (j = columnStart[iColumn];
808             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
809          int iRow = row[j];
810          rowActivity[iRow] = save[iRow];
811          mark[iRow] = 0;
812        }
813      }
814      /*
815  We've found a pair x<i> and x<k> which produce a better solution. Update our
816  notion of current solution to match.
817
818  Why does this not update newSolutionValue?
819*/
820      if (goodK >= 0) {
821        // we found something - update solution
822        for (j = columnStart[iColumn];
823             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
824          int iRow = row[j];
825          rowActivity[iRow] += wayI * element[j];
826        }
827        newSolution[iColumn] += wayI;
828        int kColumn = integerVariable[goodK];
829        for (j = columnStart[kColumn];
830             j < columnStart[kColumn] + columnLength[kColumn]; j++) {
831          int iRow = row[j];
832          rowActivity[iRow] += wayK * element[j];
833        }
834        newSolution[kColumn] += wayK;
835        /*
836  Adjust motion range for x<k>. We may have banged up against a bound with that
837  last move.
838*/
839        // See if k can go further ?
840        const OsiObject *object = model_->object(goodK);
841        // get original bounds
842        double originalLower;
843        double originalUpper;
844        getIntegerInformation(object, originalLower, originalUpper);
845
846        double value = newSolution[kColumn];
847        int iway = 0;
848
849        if (value > originalLower + 0.5)
850          iway = 1;
851        if (value < originalUpper - 0.5)
852          iway |= 2;
853        way[goodK] = static_cast< char >(iway);
854        totalChange += bestChange;
855      }
856    }
857    /*
858  End of loop to try increment/decrement of integer variables.
859
860  newSolutionValue does not necessarily match the current newSolution, and
861  bestChange simply reflects the best single change. Still, that's sufficient
862  to indicate that there's been at least one change. Check that we really do
863  have a valid solution.
864*/
865    if (totalChange + newSolutionValue < solutionValue) {
866      // paranoid check
867      memset(rowActivity, 0, numberRows * sizeof(double));
868
869      for (i = 0; i < numberColumns; i++) {
870        CoinBigIndex j;
871        double value = newSolution[i];
872        if (value) {
873          for (j = columnStart[i];
874               j < columnStart[i] + columnLength[i]; j++) {
875            int iRow = row[j];
876            rowActivity[iRow] += value * element[j];
877          }
878        }
879      }
880      int numberBad = 0;
881      double sumBad = 0.0;
882      // check was approximately feasible
883      for (i = 0; i < numberRows; i++) {
884        if (rowActivity[i] < rowLower[i]) {
885          sumBad += rowLower[i] - rowActivity[i];
886          if (rowActivity[i] < rowLower[i] - 10.0 * primalTolerance)
887            numberBad++;
888        } else if (rowActivity[i] > rowUpper[i]) {
889          sumBad += rowUpper[i] - rowActivity[i];
890          if (rowActivity[i] > rowUpper[i] + 10.0 * primalTolerance)
891            numberBad++;
892        }
893      }
894      if (!numberBad) {
895        for (i = 0; i < numberIntegers; i++) {
896          int iColumn = integerVariable[i];
897          const OsiObject *object = model_->object(i);
898          // get original bounds
899          double originalLower;
900          double originalUpper;
901          getIntegerInformation(object, originalLower, originalUpper);
902
903          double value = newSolution[iColumn];
904          // if away from lower bound mark that fact
905          if (value > originalLower) {
906            used_[iColumn] = numberSolutions_;
907          }
908        }
909        /*
910  Copy the solution to the array returned to the client. Grab a basis from
911  the solver (which, if it exists, is almost certainly infeasible, but it
912  should be ok for a dual start). The value returned as solutionValue is
913  conservative because of handling of newSolutionValue and bestChange, as
914  described above.
915*/
916        // new solution
917        memcpy(betterSolution, newSolution, numberColumns * sizeof(double));
918        CoinWarmStartBasis *basis = dynamic_cast< CoinWarmStartBasis * >(solver->getWarmStart());
919        if (basis) {
920          model_->setBestSolutionBasis(*basis);
921          delete basis;
922        }
923        returnCode = 1;
924        solutionValue = newSolutionValue + bestChange;
925      } else {
926        // bad solution - should not happen so debug if see message
927        COIN_DETAIL_PRINT(printf("Local search got bad solution with %d infeasibilities summing to %g\n",
928          numberBad, sumBad));
929      }
930    }
931    // This is just a copy!
932    delete[] integerVariable;
933  }
934  /*
935  We're done. Clean up.
936*/
937  delete[] newSolution;
938  delete[] rowActivity;
939  delete[] way;
940  delete[] cost;
941  delete[] save;
942  delete[] mark;
943  /*
944  Do we want to try swapping values between solutions?
945  swap_ is set elsewhere; it's not adjusted during heuristic execution.
946
947  Again, redundant test. We shouldn't be here if numberSolutions_ = 1.
948*/
949  if (numberSolutions_ > 1 && (swap % 10) == 1) {
950    // try merge
951    int returnCode2 = solutionFix(solutionValue, betterSolution, NULL);
952    if (returnCode2)
953      returnCode = 1;
954  }
955  return returnCode;
956}
957// update model
958void CbcHeuristicLocal::setModel(CbcModel *model)
959{
960  model_ = model;
961  // Get a copy of original matrix
962  assert(model_->solver());
963  if (model_->solver()->getNumRows()) {
964    matrix_ = *model_->solver()->getMatrixByCol();
965  }
966  delete[] used_;
967  int numberColumns = model->solver()->getNumCols();
968  used_ = new int[numberColumns];
969  memset(used_, 0, numberColumns * sizeof(int));
970}
971
972// Default Constructor
973CbcHeuristicProximity::CbcHeuristicProximity()
974  : CbcHeuristic()
975{
976  increment_ = 0.01;
977  feasibilityPump_ = NULL;
978  numberSolutions_ = 0;
979  used_ = NULL;
980  lastRunDeep_ = -1000000;
981  switches_ |= 16; // needs a new solution
982}
983
984// Constructor with model - assumed before cuts
985
986CbcHeuristicProximity::CbcHeuristicProximity(CbcModel &model)
987  : CbcHeuristic(model)
988{
989  increment_ = 0.01;
990  feasibilityPump_ = NULL;
991  numberSolutions_ = 0;
992  lastRunDeep_ = -1000000;
993  switches_ |= 16; // needs a new solution
994  int numberColumns = model.solver()->getNumCols();
995  used_ = new int[numberColumns];
996  memset(used_, 0, numberColumns * sizeof(int));
997}
998
999// Destructor
1000CbcHeuristicProximity::~CbcHeuristicProximity()
1001{
1002  delete feasibilityPump_;
1003  delete[] used_;
1004}
1005
1006// Clone
1007CbcHeuristic *
1008CbcHeuristicProximity::clone() const
1009{
1010  return new CbcHeuristicProximity(*this);
1011}
1012// Create C++ lines to get to current state
1013void CbcHeuristicProximity::generateCpp(FILE *fp)
1014{
1015  CbcHeuristicProximity other;
1016  fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n");
1017  fprintf(fp, "3  CbcHeuristicProximity heuristicProximity(*cbcModel);\n");
1018  CbcHeuristic::generateCpp(fp, "heuristicProximity");
1019  fprintf(fp, "3  cbcModel->addHeuristic(&heuristicProximity);\n");
1020}
1021
1022// Copy constructor
1023CbcHeuristicProximity::CbcHeuristicProximity(const CbcHeuristicProximity &rhs)
1024  : CbcHeuristic(rhs)
1025  , numberSolutions_(rhs.numberSolutions_)
1026{
1027  increment_ = rhs.increment_;
1028  feasibilityPump_ = NULL;
1029  if (model_ && rhs.used_) {
1030    int numberColumns = model_->solver()->getNumCols();
1031    used_ = CoinCopyOfArray(rhs.used_, numberColumns);
1032    if (rhs.feasibilityPump_)
1033      feasibilityPump_ = new CbcHeuristicFPump(*rhs.feasibilityPump_);
1034  } else {
1035    used_ = NULL;
1036  }
1037}
1038
1039// Assignment operator
1040CbcHeuristicProximity &
1041CbcHeuristicProximity::operator=(const CbcHeuristicProximity &rhs)
1042{
1043  if (this != &rhs) {
1044    CbcHeuristic::operator=(rhs);
1045    increment_ = rhs.increment_;
1046    numberSolutions_ = rhs.numberSolutions_;
1047    delete[] used_;
1048    delete feasibilityPump_;
1049    feasibilityPump_ = NULL;
1050    if (model_ && rhs.used_) {
1051      int numberColumns = model_->solver()->getNumCols();
1052      used_ = CoinCopyOfArray(rhs.used_, numberColumns);
1053      if (rhs.feasibilityPump_)
1054        feasibilityPump_ = new CbcHeuristicFPump(*rhs.feasibilityPump_);
1055    } else {
1056      used_ = NULL;
1057    }
1058  }
1059  return *this;
1060}
1061
1062// Resets stuff if model changes
1063void CbcHeuristicProximity::resetModel(CbcModel * /*model*/)
1064{
1065  //CbcHeuristic::resetModel(model);
1066  delete[] used_;
1067  if (model_ && used_) {
1068    int numberColumns = model_->solver()->getNumCols();
1069    used_ = new int[numberColumns];
1070    memset(used_, 0, numberColumns * sizeof(int));
1071  } else {
1072    used_ = NULL;
1073  }
1074}
1075/*
1076  Run a mini-BaB search after changing objective
1077
1078  Return values are:
1079    1: smallBranchAndBound found a solution
1080    0: everything else
1081
1082  The degree of overload as return codes from smallBranchAndBound are folded
1083  into 0 is such that it's impossible to distinguish return codes that really
1084  require attention from a simple `nothing of interest'.
1085*/
1086int CbcHeuristicProximity::solution(double &solutionValue,
1087  double *betterSolution)
1088{
1089  if (feasibilityPumpOptions_ == -3 && numCouldRun_ == 0 && !feasibilityPump_) {
1090    // clone feasibility pump
1091    for (int i = 0; i < model_->numberHeuristics(); i++) {
1092      const CbcHeuristicFPump *pump = dynamic_cast< const CbcHeuristicFPump * >(model_->heuristic(i));
1093      if (pump) {
1094        feasibilityPump_ = new CbcHeuristicFPump(*pump);
1095        break;
1096      }
1097    }
1098  }
1099  /*
1100  Execute only if a new solution has been discovered since the last time we
1101  were called.
1102*/
1103
1104  numCouldRun_++;
1105  int nodeCount = model_->getNodeCount();
1106  if (numberSolutions_ == model_->getSolutionCount())
1107    return 0;
1108  if (!model_->bestSolution())
1109    return 0; // odd - because in parallel mode
1110  numberSolutions_ = model_->getSolutionCount();
1111  lastRunDeep_ = nodeCount;
1112  numRuns_++;
1113  //howOftenShallow_ = numberSolutions_;
1114
1115  /*
1116  Load up a new solver with the solution.
1117
1118  Why continuousSolver(), as opposed to solver()?
1119*/
1120  OsiSolverInterface *newSolver = model_->continuousSolver()->clone();
1121  int numberColumns = newSolver->getNumCols();
1122  double *obj = CoinCopyOfArray(newSolver->getObjCoefficients(), numberColumns);
1123  int *indices = new int[numberColumns];
1124  int n = 0;
1125  for (int i = 0; i < numberColumns; i++) {
1126    if (obj[i]) {
1127      indices[n] = i;
1128      obj[n++] = obj[i];
1129    }
1130  }
1131  double cutoff = model_->getCutoff();
1132  assert(cutoff < 1.0e20);
1133  if (model_->getCutoffIncrement() < 1.0e-4) {
1134    cutoff -= increment_;
1135  }
1136  double offset;
1137  newSolver->getDblParam(OsiObjOffset, offset);
1138  newSolver->setDblParam(OsiObjOffset, 0.0);
1139  newSolver->addRow(n, indices, obj, -COIN_DBL_MAX, cutoff + offset);
1140  delete[] indices;
1141  memset(obj, 0, numberColumns * sizeof(double));
1142  newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e20);
1143  int numberIntegers = model_->numberIntegers();
1144  const int *integerVariable = model_->integerVariable();
1145  const double *solutionIn = model_->bestSolution();
1146  for (int i = 0; i < numberIntegers; i++) {
1147    int iColumn = integerVariable[i];
1148    if (!isHeuristicInteger(newSolver, iColumn))
1149      continue;
1150    if (fabs(solutionIn[iColumn]) < 1.0e-5)
1151      obj[iColumn] = 1.0;
1152    else if (fabs(solutionIn[iColumn] - 1.0) < 1.0e-5)
1153      obj[iColumn] = -1.0;
1154  }
1155  newSolver->setObjective(obj);
1156  delete[] obj;
1157  //newSolver->writeMps("xxxx");
1158  int maxSolutions = model_->getMaximumSolutions();
1159  model_->setMaximumSolutions(1);
1160  bool pumpAdded = false;
1161  if (feasibilityPumpOptions_ == -3 && feasibilityPump_) {
1162    // add back feasibility pump
1163    pumpAdded = true;
1164    for (int i = 0; i < model_->numberHeuristics(); i++) {
1165      const CbcHeuristicFPump *pump = dynamic_cast< const CbcHeuristicFPump * >(model_->heuristic(i));
1166      if (pump) {
1167        pumpAdded = false;
1168        break;
1169      }
1170    }
1171    if (pumpAdded)
1172      model_->addHeuristic(feasibilityPump_);
1173  }
1174  int returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
1175    1.0e20, "CbcHeuristicProximity");
1176  if (pumpAdded) {
1177    // take off feasibility pump
1178    int lastHeuristic = model_->numberHeuristics() - 1;
1179    model_->setNumberHeuristics(lastHeuristic);
1180    delete model_->heuristic(lastHeuristic);
1181  }
1182  model_->setMaximumSolutions(maxSolutions);
1183  /*
1184  -2 is return due to user event, and -1 is overloaded with what look to be
1185  two contradictory meanings.
1186*/
1187  if (returnCode < 0) {
1188    returnCode = 0;
1189  }
1190  /*
1191  If the result is complete exploration with a solution (3) or proven
1192  infeasibility (2), we could generate a cut (the AI folks would call it a
1193  nogood) to prevent us from going down this route in the future.
1194*/
1195  if ((returnCode & 2) != 0) {
1196    // could add cut
1197    returnCode &= ~2;
1198  }
1199  char proxPrint[200];
1200  if ((returnCode & 1) != 0) {
1201    // redo objective
1202    OsiSolverInterface *solver = model_->continuousSolver();
1203    const double *obj = solver->getObjCoefficients();
1204    solutionValue = -offset;
1205    int sumIncrease = 0.0;
1206    int sumDecrease = 0.0;
1207    int numberIncrease = 0;
1208    int numberDecrease = 0;
1209    for (int i = 0; i < numberColumns; i++) {
1210      solutionValue += obj[i] * betterSolution[i];
1211      if (isHeuristicInteger(solver, i)) {
1212        int change = static_cast< int >(floor(solutionIn[i] - betterSolution[i] + 0.5));
1213        if (change > 0) {
1214          numberIncrease++;
1215          sumIncrease += change;
1216        } else if (change < 0) {
1217          numberDecrease++;
1218          sumDecrease -= change;
1219        }
1220      }
1221    }
1222    sprintf(proxPrint, "Proximity search ran %d nodes (out of %d) - in new solution %d increased (%d), %d decreased (%d)",
1223      numberNodesDone_, numberNodes_,
1224      numberIncrease, sumIncrease, numberDecrease, sumDecrease);
1225    if (!numberIncrease && !numberDecrease) {
1226      // somehow tolerances are such that we can slip through
1227      // change for next time
1228      increment_ += CoinMax(increment_, fabs(solutionValue + offset) * 1.0e-10);
1229    }
1230  } else {
1231    sprintf(proxPrint, "Proximity search ran %d nodes - no new solution",
1232      numberNodesDone_);
1233  }
1234  model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1235    << proxPrint
1236    << CoinMessageEol;
1237
1238  delete newSolver;
1239  return returnCode;
1240}
1241// update model
1242void CbcHeuristicProximity::setModel(CbcModel *model)
1243{
1244  model_ = model;
1245  // Get a copy of original matrix
1246  assert(model_->solver());
1247  delete[] used_;
1248  int numberColumns = model->solver()->getNumCols();
1249  used_ = new int[numberColumns];
1250  memset(used_, 0, numberColumns * sizeof(int));
1251}
1252
1253// Default Constructor
1254CbcHeuristicNaive::CbcHeuristicNaive()
1255  : CbcHeuristic()
1256{
1257  large_ = 1.0e6;
1258}
1259
1260// Constructor with model - assumed before cuts
1261
1262CbcHeuristicNaive::CbcHeuristicNaive(CbcModel &model)
1263  : CbcHeuristic(model)
1264{
1265  large_ = 1.0e6;
1266}
1267
1268// Destructor
1269CbcHeuristicNaive::~CbcHeuristicNaive()
1270{
1271}
1272
1273// Clone
1274CbcHeuristic *
1275CbcHeuristicNaive::clone() const
1276{
1277  return new CbcHeuristicNaive(*this);
1278}
1279// Create C++ lines to get to current state
1280void CbcHeuristicNaive::generateCpp(FILE *fp)
1281{
1282  CbcHeuristicNaive other;
1283  fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n");
1284  fprintf(fp, "3  CbcHeuristicNaive naive(*cbcModel);\n");
1285  CbcHeuristic::generateCpp(fp, "naive");
1286  if (large_ != other.large_)
1287    fprintf(fp, "3  naive.setLarge(%g);\n", large_);
1288  else
1289    fprintf(fp, "4  naive.setLarge(%g);\n", large_);
1290  fprintf(fp, "3  cbcModel->addHeuristic(&naive);\n");
1291}
1292
1293// Copy constructor
1294CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive &rhs)
1295  : CbcHeuristic(rhs)
1296  , large_(rhs.large_)
1297{
1298}
1299
1300// Assignment operator
1301CbcHeuristicNaive &
1302CbcHeuristicNaive::operator=(const CbcHeuristicNaive &rhs)
1303{
1304  if (this != &rhs) {
1305    CbcHeuristic::operator=(rhs);
1306    large_ = rhs.large_;
1307  }
1308  return *this;
1309}
1310
1311// Resets stuff if model changes
1312void CbcHeuristicNaive::resetModel(CbcModel *model)
1313{
1314  CbcHeuristic::resetModel(model);
1315}
1316int CbcHeuristicNaive::solution(double &solutionValue,
1317  double *betterSolution)
1318{
1319  numCouldRun_++;
1320  // See if to do
1321  bool atRoot = model_->getNodeCount() == 0;
1322  int passNumber = model_->getCurrentPassNumber();
1323  if (!when() || (when() == 1 && model_->phase() != 1) || !atRoot || passNumber > 1)
1324    return 0; // switched off
1325  // Don't do if it was this heuristic which found solution!
1326  if (this == model_->lastHeuristic())
1327    return 0;
1328  numRuns_++;
1329  double cutoff;
1330  model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
1331  double direction = model_->solver()->getObjSense();
1332  cutoff *= direction;
1333  cutoff = CoinMin(cutoff, solutionValue);
1334  OsiSolverInterface *solver = model_->continuousSolver();
1335  if (!solver)
1336    solver = model_->solver();
1337  const double *colLower = solver->getColLower();
1338  const double *colUpper = solver->getColUpper();
1339  const double *objective = solver->getObjCoefficients();
1340
1341  int numberColumns = model_->getNumCols();
1342  int numberIntegers = model_->numberIntegers();
1343  const int *integerVariable = model_->integerVariable();
1344
1345  int i;
1346  bool solutionFound = false;
1347  CoinWarmStartBasis saveBasis;
1348  CoinWarmStartBasis *basis = dynamic_cast< CoinWarmStartBasis * >(solver->getWarmStart());
1349  if (basis) {
1350    saveBasis = *basis;
1351    delete basis;
1352  }
1353  // First just fix all integers as close to zero as possible
1354  OsiSolverInterface *newSolver = cloneBut(7); // wassolver->clone();
1355  for (i = 0; i < numberIntegers; i++) {
1356    int iColumn = integerVariable[i];
1357    if (!isHeuristicInteger(newSolver, iColumn))
1358      continue;
1359    double lower = colLower[iColumn];
1360    double upper = colUpper[iColumn];
1361    double value;
1362    if (lower > 0.0)
1363      value = lower;
1364    else if (upper < 0.0)
1365      value = upper;
1366    else
1367      value = 0.0;
1368    newSolver->setColLower(iColumn, value);
1369    newSolver->setColUpper(iColumn, value);
1370  }
1371  newSolver->initialSolve();
1372  if (newSolver->isProvenOptimal()) {
1373    double solValue = newSolver->getObjValue() * direction;
1374    if (solValue < cutoff) {
1375      // we have a solution
1376      solutionFound = true;
1377      solutionValue = solValue;
1378      memcpy(betterSolution, newSolver->getColSolution(),
1379        numberColumns * sizeof(double));
1380      COIN_DETAIL_PRINT(printf("Naive fixing close to zero gave solution of %g\n", solutionValue));
1381      cutoff = solValue - model_->getCutoffIncrement();
1382    }
1383  }
1384  // Now fix all integers as close to zero if not zero or large cost
1385  int nFix = 0;
1386  for (i = 0; i < numberIntegers; i++) {
1387    int iColumn = integerVariable[i];
1388    if (!isHeuristicInteger(newSolver, iColumn))
1389      continue;
1390    double lower = colLower[iColumn];
1391    double upper = colUpper[iColumn];
1392    double value;
1393    if (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) {
1394      nFix++;
1395      if (lower > 0.0)
1396        value = lower;
1397      else if (upper < 0.0)
1398        value = upper;
1399      else
1400        value = 0.0;
1401      newSolver->setColLower(iColumn, value);
1402      newSolver->setColUpper(iColumn, value);
1403    } else {
1404      // set back to original
1405      newSolver->setColLower(iColumn, lower);
1406      newSolver->setColUpper(iColumn, upper);
1407    }
1408  }
1409  const double *solution = solver->getColSolution();
1410  if (nFix) {
1411    newSolver->setWarmStart(&saveBasis);
1412    newSolver->setColSolution(solution);
1413    newSolver->initialSolve();
1414    if (newSolver->isProvenOptimal()) {
1415      double solValue = newSolver->getObjValue() * direction;
1416      if (solValue < cutoff) {
1417        // try branch and bound
1418        double *newSolution = new double[numberColumns];
1419        COIN_DETAIL_PRINT(printf("%d fixed after fixing costs\n", nFix));
1420        int returnCode = smallBranchAndBound(newSolver,
1421          numberNodes_, newSolution,
1422          solutionValue,
1423          solutionValue, "CbcHeuristicNaive1");
1424        if (returnCode < 0)
1425          returnCode = 0; // returned on size
1426        if ((returnCode & 2) != 0) {
1427          // could add cut
1428          returnCode &= ~2;
1429        }
1430        if (returnCode == 1) {
1431          // solution
1432          solutionFound = true;
1433          memcpy(betterSolution, newSolution,
1434            numberColumns * sizeof(double));
1435          COIN_DETAIL_PRINT(printf("Naive fixing zeros gave solution of %g\n", solutionValue));
1436          cutoff = solutionValue - model_->getCutoffIncrement();
1437        }
1438        delete[] newSolution;
1439      }
1440    }
1441  }
1442#if 1
1443  newSolver->setObjSense(-direction); // maximize
1444  newSolver->setWarmStart(&saveBasis);
1445  newSolver->setColSolution(solution);
1446  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1447    double value = solution[iColumn];
1448    double lower = colLower[iColumn];
1449    double upper = colUpper[iColumn];
1450    double newLower;
1451    double newUpper;
1452    if (isHeuristicInteger(newSolver, iColumn)) {
1453      newLower = CoinMax(lower, floor(value) - 2.0);
1454      newUpper = CoinMin(upper, ceil(value) + 2.0);
1455    } else {
1456      newLower = CoinMax(lower, value - 1.0e5);
1457      newUpper = CoinMin(upper, value + 1.0e-5);
1458    }
1459    newSolver->setColLower(iColumn, newLower);
1460    newSolver->setColUpper(iColumn, newUpper);
1461  }
1462  newSolver->initialSolve();
1463  if (newSolver->isProvenOptimal()) {
1464    double solValue = newSolver->getObjValue() * direction;
1465    if (solValue < cutoff) {
1466      nFix = 0;
1467      newSolver->setObjSense(direction); // correct direction
1468      //const double * thisSolution = newSolver->getColSolution();
1469      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1470        double value = solution[iColumn];
1471        double lower = colLower[iColumn];
1472        double upper = colUpper[iColumn];
1473        double newLower = lower;
1474        double newUpper = upper;
1475        if (isHeuristicInteger(newSolver, iColumn)) {
1476          if (value < lower + 1.0e-6) {
1477            nFix++;
1478            newUpper = lower;
1479          } else if (value > upper - 1.0e-6) {
1480            nFix++;
1481            newLower = upper;
1482          } else {
1483            newLower = CoinMax(lower, floor(value) - 2.0);
1484            newUpper = CoinMin(upper, ceil(value) + 2.0);
1485          }
1486        }
1487        newSolver->setColLower(iColumn, newLower);
1488        newSolver->setColUpper(iColumn, newUpper);
1489      }
1490      // try branch and bound
1491      double *newSolution = new double[numberColumns];
1492      COIN_DETAIL_PRINT(printf("%d fixed after maximizing\n", nFix));
1493      int returnCode = smallBranchAndBound(newSolver,
1494        numberNodes_, newSolution,
1495        solutionValue,
1496        solutionValue, "CbcHeuristicNaive1");
1497      if (returnCode < 0)
1498        returnCode = 0; // returned on size
1499      if ((returnCode & 2) != 0) {
1500        // could add cut
1501        returnCode &= ~2;
1502      }
1503      if (returnCode == 1) {
1504        // solution
1505        solutionFound = true;
1506        memcpy(betterSolution, newSolution,
1507          numberColumns * sizeof(double));
1508        COIN_DETAIL_PRINT(printf("Naive maximizing gave solution of %g\n", solutionValue));
1509        cutoff = solutionValue - model_->getCutoffIncrement();
1510      }
1511      delete[] newSolution;
1512    }
1513  }
1514#endif
1515  delete newSolver;
1516  return solutionFound ? 1 : 0;
1517}
1518// update model
1519void CbcHeuristicNaive::setModel(CbcModel *model)
1520{
1521  model_ = model;
1522}
1523// Default Constructor
1524CbcHeuristicCrossover::CbcHeuristicCrossover()
1525  : CbcHeuristic()
1526  , numberSolutions_(0)
1527  , useNumber_(3)
1528{
1529  setWhen(1);
1530}
1531
1532// Constructor with model - assumed before cuts
1533
1534CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel &model)
1535  : CbcHeuristic(model)
1536  , numberSolutions_(0)
1537  , useNumber_(3)
1538{
1539  setWhen(1);
1540  for (int i = 0; i < 10; i++)
1541    random_[i] = model.randomNumberGenerator()->randomDouble();
1542}
1543
1544// Destructor
1545CbcHeuristicCrossover::~CbcHeuristicCrossover()
1546{
1547}
1548
1549// Clone
1550CbcHeuristic *
1551CbcHeuristicCrossover::clone() const
1552{
1553  return new CbcHeuristicCrossover(*this);
1554}
1555// Create C++ lines to get to current state
1556void CbcHeuristicCrossover::generateCpp(FILE *fp)
1557{
1558  CbcHeuristicCrossover other;
1559  fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n");
1560  fprintf(fp, "3  CbcHeuristicCrossover crossover(*cbcModel);\n");
1561  CbcHeuristic::generateCpp(fp, "crossover");
1562  if (useNumber_ != other.useNumber_)
1563    fprintf(fp, "3  crossover.setNumberSolutions(%d);\n", useNumber_);
1564  else
1565    fprintf(fp, "4  crossover.setNumberSolutions(%d);\n", useNumber_);
1566  fprintf(fp, "3  cbcModel->addHeuristic(&crossover);\n");
1567}
1568
1569// Copy constructor
1570CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover &rhs)
1571  : CbcHeuristic(rhs)
1572  , attempts_(rhs.attempts_)
1573  , numberSolutions_(rhs.numberSolutions_)
1574  , useNumber_(rhs.useNumber_)
1575{
1576  memcpy(random_, rhs.random_, 10 * sizeof(double));
1577}
1578
1579// Assignment operator
1580CbcHeuristicCrossover &
1581CbcHeuristicCrossover::operator=(const CbcHeuristicCrossover &rhs)
1582{
1583  if (this != &rhs) {
1584    CbcHeuristic::operator=(rhs);
1585    useNumber_ = rhs.useNumber_;
1586    attempts_ = rhs.attempts_;
1587    numberSolutions_ = rhs.numberSolutions_;
1588    memcpy(random_, rhs.random_, 10 * sizeof(double));
1589  }
1590  return *this;
1591}
1592
1593// Resets stuff if model changes
1594void CbcHeuristicCrossover::resetModel(CbcModel *model)
1595{
1596  CbcHeuristic::resetModel(model);
1597}
1598int CbcHeuristicCrossover::solution(double &solutionValue,
1599  double *betterSolution)
1600{
1601  if (when_ == 0)
1602    return 0;
1603  numCouldRun_++;
1604  bool useBest = (numberSolutions_ != model_->getSolutionCount());
1605  if (!useBest && (when_ % 10) == 1)
1606    return 0;
1607  numberSolutions_ = model_->getSolutionCount();
1608  OsiSolverInterface *continuousSolver = model_->continuousSolver();
1609  int useNumber = CoinMin(model_->numberSavedSolutions(), useNumber_);
1610  if (useNumber < 2 || !continuousSolver)
1611    return 0;
1612  // Fix later
1613  if (!useBest)
1614    abort();
1615  numRuns_++;
1616  double cutoff;
1617  model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
1618  double direction = model_->solver()->getObjSense();
1619  cutoff *= direction;
1620  cutoff = CoinMin(cutoff, solutionValue);
1621  OsiSolverInterface *solver = cloneBut(2);
1622  // But reset bounds
1623  solver->setColLower(continuousSolver->getColLower());
1624  solver->setColUpper(continuousSolver->getColUpper());
1625  int numberColumns = solver->getNumCols();
1626  // Fixed
1627  double *fixed = new double[numberColumns];
1628  for (int i = 0; i < numberColumns; i++)
1629    fixed[i] = -COIN_DBL_MAX;
1630  int whichSolution[10];
1631  for (int i = 0; i < useNumber; i++)
1632    whichSolution[i] = i;
1633  for (int i = 0; i < useNumber; i++) {
1634    int k = whichSolution[i];
1635    const double *solution = model_->savedSolution(k);
1636    for (int j = 0; j < numberColumns; j++) {
1637      if (isHeuristicInteger(solver, j)) {
1638        if (fixed[j] == -COIN_DBL_MAX)
1639          fixed[j] = floor(solution[j] + 0.5);
1640        else if (fabs(fixed[j] - solution[j]) > 1.0e-7)
1641          fixed[j] = COIN_DBL_MAX;
1642      }
1643    }
1644  }
1645  const double *colLower = solver->getColLower();
1646  for (int i = 0; i < numberColumns; i++) {
1647    if (isHeuristicInteger(solver, i)) {
1648      double value = fixed[i];
1649      if (value != COIN_DBL_MAX) {
1650        if (when_ < 10) {
1651          solver->setColLower(i, value);
1652          solver->setColUpper(i, value);
1653        } else if (value == colLower[i]) {
1654          solver->setColUpper(i, value);
1655        }
1656      }
1657    }
1658  }
1659  int returnCode = smallBranchAndBound(solver, numberNodes_, betterSolution,
1660    solutionValue,
1661    solutionValue, "CbcHeuristicCrossover");
1662  if (returnCode < 0)
1663    returnCode = 0; // returned on size
1664  if ((returnCode & 2) != 0) {
1665    // could add cut
1666    returnCode &= ~2;
1667  }
1668
1669  delete[] fixed;
1670  delete solver;
1671  return returnCode;
1672}
1673// update model
1674void CbcHeuristicCrossover::setModel(CbcModel *model)
1675{
1676  model_ = model;
1677  if (model) {
1678    for (int i = 0; i < 10; i++)
1679      random_[i] = model->randomNumberGenerator()->randomDouble();
1680  }
1681}
1682
1683/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1684*/
Note: See TracBrowser for help on using the repository browser.