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

Last change on this file since 2094 was 2094, checked in by forrest, 5 years ago

for memory leaks and heuristics and some experimental stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 58.3 KB
Line 
1/* $Id: CbcHeuristicLocal.cpp 2094 2014-11-18 11:15:36Z 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())
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        if (((swap/10) &1) != 0) {
569          maxIntegers = CoinMin(1000,numberIntegers);
570        }
571/*
572  Outer loop to walk integer variables. Call the current variable x<i>. At the
573  end of this loop, bestChange will contain the best (negative) change in the
574  objective for any single pair.
575
576  The trouble is, we're limited to monotonically increasing improvement.
577  Suppose we discover an improvement of 10 for some pair. If, later in the
578  search, we discover an improvement of 9 for some other pair, we will not use
579  it. That seems wasteful.
580*/
581
582        for (i = 0; i < numberIntegers; i++) {
583            int iColumn = integerVariable[i];
584            bestChange = 0.0;
585            int endInner = CoinMin(numberIntegers,i+maxIntegers);
586
587            double objectiveCoefficient = cost[i];
588            int k;
589            int j;
590            int goodK = -1;
591            int wayK = -1, wayI = -1;
592/*
593  Try decrementing x<i>.
594*/
595            if ((way[i]&1) != 0) {
596                int numberInfeasible = 0;
597/*
598  Adjust row activities where x<i> has a nonzero coefficient. Save the old
599  values for restoration. Mark any rows that become infeasible as a result
600  of the decrement.
601*/
602                // save row activities and adjust
603                for (j = columnStart[iColumn];
604                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
605                    int iRow = row[j];
606                    save[iRow] = rowActivity[iRow];
607                    rowActivity[iRow] -= element[j];
608                    if (rowActivity[iRow] < rowLower[iRow] - primalTolerance ||
609                            rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
610                        // mark row
611                        mark[iRow] = 1;
612                        numberInfeasible++;
613                    }
614                }
615  /*
616  Run through the remaining integer variables. Try increment and decrement on
617  each one. If the potential objective change is better than anything we've
618  seen so far, do a full evaluation of x<k> in that direction.  If we can
619  repair all infeasibilities introduced by pushing x<i> down, we have a
620  winner. Remember the best variable, and the direction for x<i> and x<k>.
621*/
622              // try down
623                for (k = i + 1; k < endInner; k++) {
624                    if ((way[k]&1) != 0) {
625                        // try down
626                        if (-objectiveCoefficient - cost[k] < bestChange) {
627                            // see if feasible down
628                            bool good = true;
629                            int numberMarked = 0;
630                            int kColumn = integerVariable[k];
631                            for (j = columnStart[kColumn];
632                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
633                                int iRow = row[j];
634                                double newValue = rowActivity[iRow] - element[j];
635                                if (newValue < rowLower[iRow] - primalTolerance ||
636                                        newValue > rowUpper[iRow] + primalTolerance) {
637                                    good = false;
638                                    break;
639                                } else if (mark[iRow]) {
640                                    // made feasible
641                                    numberMarked++;
642                                }
643                            }
644                            if (good && numberMarked == numberInfeasible) {
645                                // better solution
646                                goodK = k;
647                                wayK = -1;
648                                wayI = -1;
649                                bestChange = -objectiveCoefficient - cost[k];
650                            }
651                        }
652                    }
653                    if ((way[k]&2) != 0) {
654                        // try up
655                        if (-objectiveCoefficient + cost[k] < bestChange) {
656                            // see if feasible up
657                            bool good = true;
658                            int numberMarked = 0;
659                            int kColumn = integerVariable[k];
660                            for (j = columnStart[kColumn];
661                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
662                                int iRow = row[j];
663                                double newValue = rowActivity[iRow] + element[j];
664                                if (newValue < rowLower[iRow] - primalTolerance ||
665                                        newValue > rowUpper[iRow] + primalTolerance) {
666                                    good = false;
667                                    break;
668                                } else if (mark[iRow]) {
669                                    // made feasible
670                                    numberMarked++;
671                                }
672                            }
673                            if (good && numberMarked == numberInfeasible) {
674                                // better solution
675                                goodK = k;
676                                wayK = 1;
677                                wayI = -1;
678                                bestChange = -objectiveCoefficient + cost[k];
679                            }
680                        }
681                    }
682                }
683/*
684  Remove effect of decrementing x<i> by restoring original lhs values.
685*/
686                // restore row activities
687                for (j = columnStart[iColumn];
688                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
689                    int iRow = row[j];
690                    rowActivity[iRow] = save[iRow];
691                    mark[iRow] = 0;
692                }
693            }
694/*
695  Try to increment x<i>. Actions as for decrement.
696*/
697            if ((way[i]&2) != 0) {
698                int numberInfeasible = 0;
699                // save row activities and adjust
700                for (j = columnStart[iColumn];
701                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
702                    int iRow = row[j];
703                    save[iRow] = rowActivity[iRow];
704                    rowActivity[iRow] += element[j];
705                    if (rowActivity[iRow] < rowLower[iRow] - primalTolerance ||
706                            rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
707                        // mark row
708                        mark[iRow] = 1;
709                        numberInfeasible++;
710                    }
711                }
712                // try up
713                for (k = i + 1; k < endInner; k++) {
714                    if ((way[k]&1) != 0) {
715                        // try down
716                        if (objectiveCoefficient - cost[k] < bestChange) {
717                            // see if feasible down
718                            bool good = true;
719                            int numberMarked = 0;
720                            int kColumn = integerVariable[k];
721                            for (j = columnStart[kColumn];
722                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
723                                int iRow = row[j];
724                                double newValue = rowActivity[iRow] - element[j];
725                                if (newValue < rowLower[iRow] - primalTolerance ||
726                                        newValue > rowUpper[iRow] + primalTolerance) {
727                                    good = false;
728                                    break;
729                                } else if (mark[iRow]) {
730                                    // made feasible
731                                    numberMarked++;
732                                }
733                            }
734                            if (good && numberMarked == numberInfeasible) {
735                                // better solution
736                                goodK = k;
737                                wayK = -1;
738                                wayI = 1;
739                                bestChange = objectiveCoefficient - cost[k];
740                            }
741                        }
742                    }
743                    if ((way[k]&2) != 0) {
744                        // try up
745                        if (objectiveCoefficient + cost[k] < bestChange) {
746                            // see if feasible up
747                            bool good = true;
748                            int numberMarked = 0;
749                            int kColumn = integerVariable[k];
750                            for (j = columnStart[kColumn];
751                                    j < columnStart[kColumn] + columnLength[kColumn]; j++) {
752                                int iRow = row[j];
753                                double newValue = rowActivity[iRow] + element[j];
754                                if (newValue < rowLower[iRow] - primalTolerance ||
755                                        newValue > rowUpper[iRow] + primalTolerance) {
756                                    good = false;
757                                    break;
758                                } else if (mark[iRow]) {
759                                    // made feasible
760                                    numberMarked++;
761                                }
762                            }
763                            if (good && numberMarked == numberInfeasible) {
764                                // better solution
765                                goodK = k;
766                                wayK = 1;
767                                wayI = 1;
768                                bestChange = objectiveCoefficient + cost[k];
769                            }
770                        }
771                    }
772                }
773                // restore row activities
774                for (j = columnStart[iColumn];
775                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
776                    int iRow = row[j];
777                    rowActivity[iRow] = save[iRow];
778                    mark[iRow] = 0;
779                }
780            }
781/*
782  We've found a pair x<i> and x<k> which produce a better solution. Update our
783  notion of current solution to match.
784
785  Why does this not update newSolutionValue?
786*/
787            if (goodK >= 0) {
788                // we found something - update solution
789                for (j = columnStart[iColumn];
790                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
791                    int iRow = row[j];
792                    rowActivity[iRow]  += wayI * element[j];
793                }
794                newSolution[iColumn] += wayI;
795                int kColumn = integerVariable[goodK];
796                for (j = columnStart[kColumn];
797                        j < columnStart[kColumn] + columnLength[kColumn]; j++) {
798                    int iRow = row[j];
799                    rowActivity[iRow]  += wayK * element[j];
800                }
801                newSolution[kColumn] += wayK;
802/*
803  Adjust motion range for x<k>. We may have banged up against a bound with that
804  last move.
805*/
806               // See if k can go further ?
807                const OsiObject * object = model_->object(goodK);
808                // get original bounds
809                double originalLower;
810                double originalUpper;
811                getIntegerInformation( object, originalLower, originalUpper);
812
813                double value = newSolution[kColumn];
814                int iway = 0;
815
816                if (value > originalLower + 0.5)
817                    iway = 1;
818                if (value < originalUpper - 0.5)
819                    iway |= 2;
820                way[goodK] = static_cast<char>(iway);
821                totalChange += bestChange;
822            }
823        }
824/*
825  End of loop to try increment/decrement of integer variables.
826
827  newSolutionValue does not necessarily match the current newSolution, and
828  bestChange simply reflects the best single change. Still, that's sufficient
829  to indicate that there's been at least one change. Check that we really do
830  have a valid solution.
831*/
832        if (totalChange + newSolutionValue < solutionValue) {
833            // paranoid check
834            memset(rowActivity, 0, numberRows*sizeof(double));
835
836            for (i = 0; i < numberColumns; i++) {
837                int j;
838                double value = newSolution[i];
839                if (value) {
840                    for (j = columnStart[i];
841                            j < columnStart[i] + columnLength[i]; j++) {
842                        int iRow = row[j];
843                        rowActivity[iRow] += value * element[j];
844                    }
845                }
846            }
847            int numberBad = 0;
848            double sumBad = 0.0;
849            // check was approximately feasible
850            for (i = 0; i < numberRows; i++) {
851                if (rowActivity[i] < rowLower[i]) {
852                    sumBad += rowLower[i] - rowActivity[i];
853                    if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance)
854                        numberBad++;
855                } else if (rowActivity[i] > rowUpper[i]) {
856                    sumBad += rowUpper[i] - rowActivity[i];
857                    if (rowActivity[i] > rowUpper[i] + 10.0*primalTolerance)
858                        numberBad++;
859                }
860            }
861            if (!numberBad) {
862                for (i = 0; i < numberIntegers; i++) {
863                    int iColumn = integerVariable[i];
864                    const OsiObject * object = model_->object(i);
865                    // get original bounds
866                    double originalLower;
867                    double originalUpper;
868                    getIntegerInformation( object, originalLower, originalUpper);
869
870                    double value = newSolution[iColumn];
871                    // if away from lower bound mark that fact
872                    if (value > originalLower) {
873                        used_[iColumn] = numberSolutions_;
874                    }
875                }
876/*
877  Copy the solution to the array returned to the client. Grab a basis from
878  the solver (which, if it exists, is almost certainly infeasible, but it
879  should be ok for a dual start). The value returned as solutionValue is
880  conservative because of handling of newSolutionValue and bestChange, as
881  described above.
882*/
883                // new solution
884                memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
885                CoinWarmStartBasis * basis =
886                    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
887                if (basis) {
888                    model_->setBestSolutionBasis(* basis);
889                    delete basis;
890                }
891                returnCode = 1;
892                solutionValue = newSolutionValue + bestChange;
893            } else {
894                // bad solution - should not happen so debug if see message
895                COIN_DETAIL_PRINT(printf("Local search got bad solution with %d infeasibilities summing to %g\n",
896                                         numberBad, sumBad));
897            }
898        }
899    }
900/*
901  We're done. Clean up.
902*/
903    delete [] newSolution;
904    delete [] rowActivity;
905    delete [] way;
906    delete [] cost;
907    delete [] save;
908    delete [] mark;
909/*
910  Do we want to try swapping values between solutions?
911  swap_ is set elsewhere; it's not adjusted during heuristic execution.
912
913  Again, redundant test. We shouldn't be here if numberSolutions_ = 1.
914*/
915    if (numberSolutions_ > 1 && (swap%10) == 1) {
916        // try merge
917        int returnCode2 = solutionFix( solutionValue, betterSolution, NULL);
918        if (returnCode2)
919            returnCode = 1;
920    }
921    return returnCode;
922}
923// update model
924void CbcHeuristicLocal::setModel(CbcModel * model)
925{
926    model_ = model;
927    // Get a copy of original matrix
928    assert(model_->solver());
929    if (model_->solver()->getNumRows()) {
930        matrix_ = *model_->solver()->getMatrixByCol();
931    }
932    delete [] used_;
933    int numberColumns = model->solver()->getNumCols();
934    used_ = new int[numberColumns];
935    memset(used_, 0, numberColumns*sizeof(int));
936}
937
938// Default Constructor
939CbcHeuristicProximity::CbcHeuristicProximity()
940        : CbcHeuristic()
941{
942    increment_ = 0.01;
943    feasibilityPump_ = NULL;
944    numberSolutions_ = 0;
945    used_ = NULL;
946    lastRunDeep_ = -1000000;
947    switches_ |= 16; // needs a new solution
948}
949
950// Constructor with model - assumed before cuts
951
952CbcHeuristicProximity::CbcHeuristicProximity(CbcModel & model)
953        : CbcHeuristic(model)
954{
955    increment_ = 0.01;
956    feasibilityPump_ = NULL;
957    numberSolutions_ = 0;
958    lastRunDeep_ = -1000000;
959    switches_ |= 16; // needs a new solution
960    int numberColumns = model.solver()->getNumCols();
961    used_ = new int[numberColumns];
962    memset(used_, 0, numberColumns*sizeof(int));
963}
964
965// Destructor
966CbcHeuristicProximity::~CbcHeuristicProximity ()
967{
968    delete feasibilityPump_;
969    delete [] used_;
970}
971
972// Clone
973CbcHeuristic *
974CbcHeuristicProximity::clone() const
975{
976    return new CbcHeuristicProximity(*this);
977}
978// Create C++ lines to get to current state
979void
980CbcHeuristicProximity::generateCpp( FILE * fp)
981{
982    CbcHeuristicProximity other;
983    fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n");
984    fprintf(fp, "3  CbcHeuristicProximity heuristicProximity(*cbcModel);\n");
985    CbcHeuristic::generateCpp(fp, "heuristicProximity");
986    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicProximity);\n");
987}
988
989// Copy constructor
990CbcHeuristicProximity::CbcHeuristicProximity(const CbcHeuristicProximity & rhs)
991  :
992  CbcHeuristic(rhs),
993  numberSolutions_(rhs.numberSolutions_)
994{
995    increment_ = rhs.increment_;
996    feasibilityPump_ = NULL;
997    if (model_ && rhs.used_) {
998        int numberColumns = model_->solver()->getNumCols();
999        used_ = CoinCopyOfArray(rhs.used_, numberColumns);
1000        if (rhs.feasibilityPump_)
1001          feasibilityPump_ = new CbcHeuristicFPump(*rhs.feasibilityPump_);
1002    } else {
1003        used_ = NULL;
1004    }
1005}
1006
1007// Assignment operator
1008CbcHeuristicProximity &
1009CbcHeuristicProximity::operator=( const CbcHeuristicProximity & rhs)
1010{
1011    if (this != &rhs) {
1012        CbcHeuristic::operator=(rhs);
1013        increment_ = rhs.increment_;
1014        numberSolutions_ = rhs.numberSolutions_;
1015        delete [] used_;
1016        delete feasibilityPump_;
1017        feasibilityPump_ = NULL;
1018        if (model_ && rhs.used_) {
1019            int numberColumns = model_->solver()->getNumCols();
1020            used_ = CoinCopyOfArray(rhs.used_, numberColumns);
1021            if (rhs.feasibilityPump_)
1022              feasibilityPump_ = new CbcHeuristicFPump(*rhs.feasibilityPump_);
1023        } else {
1024            used_ = NULL;
1025        }
1026    }
1027    return *this;
1028}
1029
1030// Resets stuff if model changes
1031void
1032CbcHeuristicProximity::resetModel(CbcModel * /*model*/)
1033{
1034    //CbcHeuristic::resetModel(model);
1035    delete [] used_;
1036    if (model_ && used_) {
1037        int numberColumns = model_->solver()->getNumCols();
1038        used_ = new int[numberColumns];
1039        memset(used_, 0, numberColumns*sizeof(int));
1040    } else {
1041        used_ = NULL;
1042    }
1043}
1044/*
1045  Run a mini-BaB search after changing objective
1046
1047  Return values are:
1048    1: smallBranchAndBound found a solution
1049    0: everything else
1050
1051  The degree of overload as return codes from smallBranchAndBound are folded
1052  into 0 is such that it's impossible to distinguish return codes that really
1053  require attention from a simple `nothing of interest'.
1054*/
1055int
1056CbcHeuristicProximity::solution(double & solutionValue,
1057                            double * betterSolution)
1058{
1059  if (feasibilityPumpOptions_ == -3 && numCouldRun_==0 &&
1060      !feasibilityPump_ ) {
1061    // clone feasibility pump
1062    for (int i = 0; i < model_->numberHeuristics(); i++) {
1063      const CbcHeuristicFPump* pump =
1064        dynamic_cast<const CbcHeuristicFPump*>(model_->heuristic(i));
1065      if (pump) {
1066        feasibilityPump_ = new CbcHeuristicFPump(*pump);
1067        break;
1068      }
1069    }
1070  }
1071/*
1072  Execute only if a new solution has been discovered since the last time we
1073  were called.
1074*/
1075
1076  numCouldRun_++;
1077  int nodeCount = model_->getNodeCount();
1078  if (numberSolutions_ == model_->getSolutionCount())
1079    return 0;
1080  if (!model_->bestSolution())
1081    return 0; // odd - because in parallel mode
1082  numberSolutions_ = model_->getSolutionCount();
1083  lastRunDeep_ = nodeCount;
1084  numRuns_++;
1085  //howOftenShallow_ = numberSolutions_;
1086 
1087/*
1088  Load up a new solver with the solution.
1089
1090  Why continuousSolver(), as opposed to solver()?
1091*/
1092  OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1093  int numberColumns=newSolver->getNumCols();
1094  double * obj = CoinCopyOfArray(newSolver->getObjCoefficients(),numberColumns);
1095  int * indices = new int [numberColumns];
1096  int n=0;
1097  for (int i=0;i<numberColumns;i++) {
1098    if (obj[i]) {
1099      indices[n]=i;
1100      obj[n++]=obj[i];
1101    }
1102  }
1103  double cutoff=model_->getCutoff();
1104  assert (cutoff<1.0e20);
1105  if (model_->getCutoffIncrement()<1.0e-4) {
1106    cutoff -= increment_;
1107  }
1108  double offset;
1109  newSolver->getDblParam(OsiObjOffset, offset);
1110  newSolver->setDblParam(OsiObjOffset, 0.0);
1111  newSolver->addRow(n,indices,obj,-COIN_DBL_MAX,cutoff+offset);
1112  delete [] indices;
1113  memset(obj,0,numberColumns*sizeof(double));
1114  newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e20);
1115  int numberIntegers = model_->numberIntegers();
1116  const int * integerVariable = model_->integerVariable();
1117  const double * solutionIn = model_->bestSolution();
1118  for (int i = 0; i < numberIntegers; i++) {
1119    int iColumn = integerVariable[i];
1120    if (fabs(solutionIn[iColumn])<1.0e-5) 
1121      obj[iColumn]=1.0;
1122    else if (fabs(solutionIn[iColumn]-1.0)<1.0e-5) 
1123      obj[iColumn]=-1.0;
1124  }
1125  newSolver->setObjective(obj);
1126  delete [] obj;
1127  //newSolver->writeMps("xxxx");
1128  int maxSolutions = model_->getMaximumSolutions();
1129  model_->setMaximumSolutions(1); 
1130  bool pumpAdded = false;
1131  if (feasibilityPumpOptions_ == -3 && feasibilityPump_) {
1132    // add back feasibility pump
1133    pumpAdded = true;
1134    for (int i = 0; i < model_->numberHeuristics(); i++) {
1135      const CbcHeuristicFPump* pump =
1136        dynamic_cast<const CbcHeuristicFPump*>(model_->heuristic(i));
1137      if (pump) {
1138        pumpAdded = false;
1139        break;
1140      }
1141    }
1142    if (pumpAdded) 
1143      model_->addHeuristic(feasibilityPump_);
1144  }
1145  int returnCode = 
1146    smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
1147                        1.0e20, "CbcHeuristicProximity");
1148  if (pumpAdded) {
1149    // take off feasibility pump
1150    int lastHeuristic = model_->numberHeuristics()-1;
1151    model_->setNumberHeuristics(lastHeuristic);
1152    delete model_->heuristic(lastHeuristic);
1153  }
1154  model_->setMaximumSolutions(maxSolutions); 
1155 /*
1156  -2 is return due to user event, and -1 is overloaded with what look to be
1157  two contradictory meanings.
1158*/
1159  if (returnCode < 0) {
1160    returnCode = 0; 
1161  }
1162/*
1163  If the result is complete exploration with a solution (3) or proven
1164  infeasibility (2), we could generate a cut (the AI folks would call it a
1165  nogood) to prevent us from going down this route in the future.
1166*/
1167  if ((returnCode&2) != 0) {
1168    // could add cut
1169    returnCode &= ~2;
1170  }
1171  char proxPrint[200];
1172  if ((returnCode&1) != 0) {
1173    // redo objective
1174    const double * obj = model_->continuousSolver()->getObjCoefficients();
1175    solutionValue = - offset;
1176    int sumIncrease=0.0;
1177    int sumDecrease=0.0;
1178    int numberIncrease=0;
1179    int numberDecrease=0;
1180    for (int i=0;i<numberColumns;i++) {
1181      solutionValue += obj[i]*betterSolution[i];
1182      if (model_->isInteger(i)) {
1183        int change=static_cast<int>(floor(solutionIn[i]-betterSolution[i]+0.5));
1184        if (change>0) {
1185          numberIncrease++;
1186          sumIncrease+=change;
1187        } else if (change<0) {
1188          numberDecrease++;
1189          sumDecrease-=change;
1190        }
1191      }
1192    }
1193    sprintf(proxPrint,"Proximity search ran %d nodes (out of %d) - in new solution %d increased (%d), %d decreased (%d)",
1194            numberNodesDone_,numberNodes_,
1195            numberIncrease,sumIncrease,numberDecrease,sumDecrease);
1196    if (!numberIncrease&&!numberDecrease) {
1197      // somehow tolerances are such that we can slip through
1198      // change for next time
1199      increment_ += CoinMax(increment_,fabs(solutionValue+offset)*1.0e-10);
1200    }
1201  } else {
1202    sprintf(proxPrint,"Proximity search ran %d nodes - no new solution",
1203            numberNodesDone_);
1204  }
1205  model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1206    << proxPrint
1207    << CoinMessageEol;
1208 
1209  delete newSolver;
1210  return returnCode;
1211}
1212// update model
1213void CbcHeuristicProximity::setModel(CbcModel * model)
1214{
1215    model_ = model;
1216    // Get a copy of original matrix
1217    assert(model_->solver());
1218    delete [] used_;
1219    int numberColumns = model->solver()->getNumCols();
1220    used_ = new int[numberColumns];
1221    memset(used_, 0, numberColumns*sizeof(int));
1222}
1223
1224// Default Constructor
1225CbcHeuristicNaive::CbcHeuristicNaive()
1226        : CbcHeuristic()
1227{
1228    large_ = 1.0e6;
1229}
1230
1231// Constructor with model - assumed before cuts
1232
1233CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model)
1234        : CbcHeuristic(model)
1235{
1236    large_ = 1.0e6;
1237}
1238
1239// Destructor
1240CbcHeuristicNaive::~CbcHeuristicNaive ()
1241{
1242}
1243
1244// Clone
1245CbcHeuristic *
1246CbcHeuristicNaive::clone() const
1247{
1248    return new CbcHeuristicNaive(*this);
1249}
1250// Create C++ lines to get to current state
1251void
1252CbcHeuristicNaive::generateCpp( FILE * fp)
1253{
1254    CbcHeuristicNaive other;
1255    fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n");
1256    fprintf(fp, "3  CbcHeuristicNaive naive(*cbcModel);\n");
1257    CbcHeuristic::generateCpp(fp, "naive");
1258    if (large_ != other.large_)
1259        fprintf(fp, "3  naive.setLarge(%g);\n", large_);
1260    else
1261        fprintf(fp, "4  naive.setLarge(%g);\n", large_);
1262    fprintf(fp, "3  cbcModel->addHeuristic(&naive);\n");
1263}
1264
1265// Copy constructor
1266CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs)
1267        :
1268        CbcHeuristic(rhs),
1269        large_(rhs.large_)
1270{
1271}
1272
1273// Assignment operator
1274CbcHeuristicNaive &
1275CbcHeuristicNaive::operator=( const CbcHeuristicNaive & rhs)
1276{
1277    if (this != &rhs) {
1278        CbcHeuristic::operator=(rhs);
1279        large_ = rhs.large_;
1280    }
1281    return *this;
1282}
1283
1284// Resets stuff if model changes
1285void
1286CbcHeuristicNaive::resetModel(CbcModel * model)
1287{
1288    CbcHeuristic::resetModel(model);
1289}
1290int
1291CbcHeuristicNaive::solution(double & solutionValue,
1292                            double * betterSolution)
1293{
1294    numCouldRun_++;
1295    // See if to do
1296    bool atRoot = model_->getNodeCount() == 0;
1297    int passNumber = model_->getCurrentPassNumber();
1298    if (!when() || (when() == 1 && model_->phase() != 1) || !atRoot || passNumber > 1)
1299        return 0; // switched off
1300    // Don't do if it was this heuristic which found solution!
1301    if (this == model_->lastHeuristic())
1302        return 0;
1303    numRuns_++;
1304    double cutoff;
1305    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
1306    double direction = model_->solver()->getObjSense();
1307    cutoff *= direction;
1308    cutoff = CoinMin(cutoff, solutionValue);
1309    OsiSolverInterface * solver = model_->continuousSolver();
1310    if (!solver)
1311        solver = model_->solver();
1312    const double * colLower = solver->getColLower();
1313    const double * colUpper = solver->getColUpper();
1314    const double * objective = solver->getObjCoefficients();
1315
1316    int numberColumns = model_->getNumCols();
1317    int numberIntegers = model_->numberIntegers();
1318    const int * integerVariable = model_->integerVariable();
1319
1320    int i;
1321    bool solutionFound = false;
1322    CoinWarmStartBasis saveBasis;
1323    CoinWarmStartBasis * basis =
1324        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1325    if (basis) {
1326        saveBasis = * basis;
1327        delete basis;
1328    }
1329    // First just fix all integers as close to zero as possible
1330    OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone();
1331    for (i = 0; i < numberIntegers; i++) {
1332        int iColumn = integerVariable[i];
1333        double lower = colLower[iColumn];
1334        double upper = colUpper[iColumn];
1335        double value;
1336        if (lower > 0.0)
1337            value = lower;
1338        else if (upper < 0.0)
1339            value = upper;
1340        else
1341            value = 0.0;
1342        newSolver->setColLower(iColumn, value);
1343        newSolver->setColUpper(iColumn, value);
1344    }
1345    newSolver->initialSolve();
1346    if (newSolver->isProvenOptimal()) {
1347        double solValue = newSolver->getObjValue() * direction ;
1348        if (solValue < cutoff) {
1349            // we have a solution
1350            solutionFound = true;
1351            solutionValue = solValue;
1352            memcpy(betterSolution, newSolver->getColSolution(),
1353                   numberColumns*sizeof(double));
1354            COIN_DETAIL_PRINT(printf("Naive fixing close to zero gave solution of %g\n", solutionValue));
1355            cutoff = solValue - model_->getCutoffIncrement();
1356        }
1357    }
1358    // Now fix all integers as close to zero if not zero or large cost
1359    int nFix = 0;
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 (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) {
1366            nFix++;
1367            if (lower > 0.0)
1368                value = lower;
1369            else if (upper < 0.0)
1370                value = upper;
1371            else
1372                value = 0.0;
1373            newSolver->setColLower(iColumn, value);
1374            newSolver->setColUpper(iColumn, value);
1375        } else {
1376            // set back to original
1377            newSolver->setColLower(iColumn, lower);
1378            newSolver->setColUpper(iColumn, upper);
1379        }
1380    }
1381    const double * solution = solver->getColSolution();
1382    if (nFix) {
1383        newSolver->setWarmStart(&saveBasis);
1384        newSolver->setColSolution(solution);
1385        newSolver->initialSolve();
1386        if (newSolver->isProvenOptimal()) {
1387            double solValue = newSolver->getObjValue() * direction ;
1388            if (solValue < cutoff) {
1389                // try branch and bound
1390                double * newSolution = new double [numberColumns];
1391                COIN_DETAIL_PRINT(printf("%d fixed after fixing costs\n", nFix));
1392                int returnCode = smallBranchAndBound(newSolver,
1393                                                     numberNodes_, newSolution,
1394                                                     solutionValue,
1395                                                     solutionValue, "CbcHeuristicNaive1");
1396                if (returnCode < 0)
1397                    returnCode = 0; // returned on size
1398                if ((returnCode&2) != 0) {
1399                    // could add cut
1400                    returnCode &= ~2;
1401                }
1402                if (returnCode == 1) {
1403                    // solution
1404                    solutionFound = true;
1405                    memcpy(betterSolution, newSolution,
1406                           numberColumns*sizeof(double));
1407                    COIN_DETAIL_PRINT(printf("Naive fixing zeros gave solution of %g\n", solutionValue));
1408                    cutoff = solutionValue - model_->getCutoffIncrement();
1409                }
1410                delete [] newSolution;
1411            }
1412        }
1413    }
1414#if 1
1415    newSolver->setObjSense(-direction); // maximize
1416    newSolver->setWarmStart(&saveBasis);
1417    newSolver->setColSolution(solution);
1418    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1419        double value = solution[iColumn];
1420        double lower = colLower[iColumn];
1421        double upper = colUpper[iColumn];
1422        double newLower;
1423        double newUpper;
1424        if (newSolver->isInteger(iColumn)) {
1425            newLower = CoinMax(lower, floor(value) - 2.0);
1426            newUpper = CoinMin(upper, ceil(value) + 2.0);
1427        } else {
1428            newLower = CoinMax(lower, value - 1.0e5);
1429            newUpper = CoinMin(upper, value + 1.0e-5);
1430        }
1431        newSolver->setColLower(iColumn, newLower);
1432        newSolver->setColUpper(iColumn, newUpper);
1433    }
1434    newSolver->initialSolve();
1435    if (newSolver->isProvenOptimal()) {
1436        double solValue = newSolver->getObjValue() * direction ;
1437        if (solValue < cutoff) {
1438            nFix = 0;
1439            newSolver->setObjSense(direction); // correct direction
1440            //const double * thisSolution = newSolver->getColSolution();
1441            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1442                double value = solution[iColumn];
1443                double lower = colLower[iColumn];
1444                double upper = colUpper[iColumn];
1445                double newLower = lower;
1446                double newUpper = upper;
1447                if (newSolver->isInteger(iColumn)) {
1448                    if (value < lower + 1.0e-6) {
1449                        nFix++;
1450                        newUpper = lower;
1451                    } else if (value > upper - 1.0e-6) {
1452                        nFix++;
1453                        newLower = upper;
1454                    } else {
1455                        newLower = CoinMax(lower, floor(value) - 2.0);
1456                        newUpper = CoinMin(upper, ceil(value) + 2.0);
1457                    }
1458                }
1459                newSolver->setColLower(iColumn, newLower);
1460                newSolver->setColUpper(iColumn, newUpper);
1461            }
1462            // try branch and bound
1463            double * newSolution = new double [numberColumns];
1464            COIN_DETAIL_PRINT(printf("%d fixed after maximizing\n", nFix));
1465            int returnCode = smallBranchAndBound(newSolver,
1466                                                 numberNodes_, newSolution,
1467                                                 solutionValue,
1468                                                 solutionValue, "CbcHeuristicNaive1");
1469            if (returnCode < 0)
1470                returnCode = 0; // returned on size
1471            if ((returnCode&2) != 0) {
1472                // could add cut
1473                returnCode &= ~2;
1474            }
1475            if (returnCode == 1) {
1476                // solution
1477                solutionFound = true;
1478                memcpy(betterSolution, newSolution,
1479                       numberColumns*sizeof(double));
1480                COIN_DETAIL_PRINT(printf("Naive maximizing gave solution of %g\n", solutionValue));
1481                cutoff = solutionValue - model_->getCutoffIncrement();
1482            }
1483            delete [] newSolution;
1484        }
1485    }
1486#endif
1487    delete newSolver;
1488    return solutionFound ? 1 : 0;
1489}
1490// update model
1491void CbcHeuristicNaive::setModel(CbcModel * model)
1492{
1493    model_ = model;
1494}
1495// Default Constructor
1496CbcHeuristicCrossover::CbcHeuristicCrossover()
1497        : CbcHeuristic(),
1498        numberSolutions_(0),
1499        useNumber_(3)
1500{
1501    setWhen(1);
1502}
1503
1504// Constructor with model - assumed before cuts
1505
1506CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel & model)
1507        : CbcHeuristic(model),
1508        numberSolutions_(0),
1509        useNumber_(3)
1510{
1511    setWhen(1);
1512    for (int i = 0; i < 10; i++)
1513        random_[i] = model.randomNumberGenerator()->randomDouble();
1514}
1515
1516// Destructor
1517CbcHeuristicCrossover::~CbcHeuristicCrossover ()
1518{
1519}
1520
1521// Clone
1522CbcHeuristic *
1523CbcHeuristicCrossover::clone() const
1524{
1525    return new CbcHeuristicCrossover(*this);
1526}
1527// Create C++ lines to get to current state
1528void
1529CbcHeuristicCrossover::generateCpp( FILE * fp)
1530{
1531    CbcHeuristicCrossover other;
1532    fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n");
1533    fprintf(fp, "3  CbcHeuristicCrossover crossover(*cbcModel);\n");
1534    CbcHeuristic::generateCpp(fp, "crossover");
1535    if (useNumber_ != other.useNumber_)
1536        fprintf(fp, "3  crossover.setNumberSolutions(%d);\n", useNumber_);
1537    else
1538        fprintf(fp, "4  crossover.setNumberSolutions(%d);\n", useNumber_);
1539    fprintf(fp, "3  cbcModel->addHeuristic(&crossover);\n");
1540}
1541
1542// Copy constructor
1543CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover & rhs)
1544        :
1545        CbcHeuristic(rhs),
1546        attempts_(rhs.attempts_),
1547        numberSolutions_(rhs.numberSolutions_),
1548        useNumber_(rhs.useNumber_)
1549{
1550    memcpy(random_, rhs.random_, 10*sizeof(double));
1551}
1552
1553// Assignment operator
1554CbcHeuristicCrossover &
1555CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover & rhs)
1556{
1557    if (this != &rhs) {
1558        CbcHeuristic::operator=(rhs);
1559        useNumber_ = rhs.useNumber_;
1560        attempts_ = rhs.attempts_;
1561        numberSolutions_ = rhs.numberSolutions_;
1562        memcpy(random_, rhs.random_, 10*sizeof(double));
1563    }
1564    return *this;
1565}
1566
1567// Resets stuff if model changes
1568void
1569CbcHeuristicCrossover::resetModel(CbcModel * model)
1570{
1571    CbcHeuristic::resetModel(model);
1572}
1573int
1574CbcHeuristicCrossover::solution(double & solutionValue,
1575                                double * betterSolution)
1576{
1577    if (when_ == 0)
1578        return 0;
1579    numCouldRun_++;
1580    bool useBest = (numberSolutions_ != model_->getSolutionCount());
1581    if (!useBest && (when_ % 10) == 1)
1582        return 0;
1583    numberSolutions_ = model_->getSolutionCount();
1584    OsiSolverInterface * continuousSolver = model_->continuousSolver();
1585    int useNumber = CoinMin(model_->numberSavedSolutions(), useNumber_);
1586    if (useNumber < 2 || !continuousSolver)
1587        return 0;
1588    // Fix later
1589    if (!useBest)
1590        abort();
1591    numRuns_++;
1592    double cutoff;
1593    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
1594    double direction = model_->solver()->getObjSense();
1595    cutoff *= direction;
1596    cutoff = CoinMin(cutoff, solutionValue);
1597    OsiSolverInterface * solver = cloneBut(2);
1598    // But reset bounds
1599    solver->setColLower(continuousSolver->getColLower());
1600    solver->setColUpper(continuousSolver->getColUpper());
1601    int numberColumns = solver->getNumCols();
1602    // Fixed
1603    double * fixed = new double [numberColumns];
1604    for (int i = 0; i < numberColumns; i++)
1605        fixed[i] = -COIN_DBL_MAX;
1606    int whichSolution[10];
1607    for (int i = 0; i < useNumber; i++)
1608        whichSolution[i] = i;
1609    for (int i = 0; i < useNumber; i++) {
1610        int k = whichSolution[i];
1611        const double * solution = model_->savedSolution(k);
1612        for (int j = 0; j < numberColumns; j++) {
1613            if (solver->isInteger(j)) {
1614                if (fixed[j] == -COIN_DBL_MAX)
1615                    fixed[j] = floor(solution[j] + 0.5);
1616                else if (fabs(fixed[j] - solution[j]) > 1.0e-7)
1617                    fixed[j] = COIN_DBL_MAX;
1618            }
1619        }
1620    }
1621    const double * colLower = solver->getColLower();
1622    for (int i = 0; i < numberColumns; i++) {
1623        if (solver->isInteger(i)) {
1624            double value = fixed[i];
1625            if (value != COIN_DBL_MAX) {
1626                if (when_ < 10) {
1627                    solver->setColLower(i, value);
1628                    solver->setColUpper(i, value);
1629                } else if (value == colLower[i]) {
1630                    solver->setColUpper(i, value);
1631                }
1632            }
1633        }
1634    }
1635    int returnCode = smallBranchAndBound(solver, numberNodes_, betterSolution,
1636                                         solutionValue,
1637                                         solutionValue, "CbcHeuristicCrossover");
1638    if (returnCode < 0)
1639        returnCode = 0; // returned on size
1640    if ((returnCode&2) != 0) {
1641        // could add cut
1642        returnCode &= ~2;
1643    }
1644
1645    delete solver;
1646    return returnCode;
1647}
1648// update model
1649void CbcHeuristicCrossover::setModel(CbcModel * model)
1650{
1651    model_ = model;
1652    if (model) {
1653        for (int i = 0; i < 10; i++)
1654            random_[i] = model->randomNumberGenerator()->randomDouble();
1655    }
1656}
1657
1658
Note: See TracBrowser for help on using the repository browser.