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

Last change on this file since 2280 was 2280, checked in by forrest, 3 years ago

allow heuristics to see if integers are 'optional'

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