source: stable/2.8/Cbc/src/CbcHeuristicLocal.cpp @ 1820

Last change on this file since 1820 was 1820, checked in by stefan, 7 years ago

sync with trunk rev 1819

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