source: trunk/Cbc/src/CbcHeuristicLocal.cpp @ 2105

Last change on this file since 2105 was 2105, checked in by forrest, 4 years ago

mostly reporting options plus a few tweaks

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