source: trunk/Cbc/src/CbcHeuristicFPump.cpp @ 1802

Last change on this file since 1802 was 1802, checked in by forrest, 7 years ago

add Proximity heuristic (Fischetti and Monaci) - shouldn't break anything

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 131.5 KB
Line 
1/* $Id: CbcHeuristicFPump.cpp 1802 2012-11-21 09:38:56Z forrest $ */
2// Copyright (C) 2004, 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#ifdef COIN_HAS_CLP
18#include "OsiClpSolverInterface.hpp"
19#endif
20#include "CbcMessage.hpp"
21#include "CbcHeuristicFPump.hpp"
22#include "CbcBranchActual.hpp"
23#include "CbcBranchDynamic.hpp"
24#include "CoinHelperFunctions.hpp"
25#include "CoinWarmStartBasis.hpp"
26#include "CoinTime.hpp"
27#include "CbcEventHandler.hpp"
28
29
30// Default Constructor
31CbcHeuristicFPump::CbcHeuristicFPump()
32        : CbcHeuristic(),
33        startTime_(0.0),
34        maximumTime_(0.0),
35        fakeCutoff_(COIN_DBL_MAX),
36        absoluteIncrement_(0.0),
37        relativeIncrement_(0.0),
38        defaultRounding_(0.49999),
39        initialWeight_(0.0),
40        weightFactor_(0.1),
41        artificialCost_(COIN_DBL_MAX),
42        iterationRatio_(0.0),
43        reducedCostMultiplier_(1.0),
44        maximumPasses_(100),
45        maximumRetries_(1),
46        accumulate_(0),
47        fixOnReducedCosts_(1),
48        roundExpensive_(false)
49{
50    setWhen(1);
51}
52
53// Constructor from model
54CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model,
55                                     double downValue, bool roundExpensive)
56        : CbcHeuristic(model),
57        startTime_(0.0),
58        maximumTime_(0.0),
59        fakeCutoff_(COIN_DBL_MAX),
60        absoluteIncrement_(0.0),
61        relativeIncrement_(0.0),
62        defaultRounding_(downValue),
63        initialWeight_(0.0),
64        weightFactor_(0.1),
65        artificialCost_(COIN_DBL_MAX),
66        iterationRatio_(0.0),
67        reducedCostMultiplier_(1.0),
68        maximumPasses_(100),
69        maximumRetries_(1),
70        accumulate_(0),
71        fixOnReducedCosts_(1),
72        roundExpensive_(roundExpensive)
73{
74    setWhen(1);
75}
76
77// Destructor
78CbcHeuristicFPump::~CbcHeuristicFPump ()
79{
80}
81
82// Clone
83CbcHeuristic *
84CbcHeuristicFPump::clone() const
85{
86    return new CbcHeuristicFPump(*this);
87}
88// Create C++ lines to get to current state
89void
90CbcHeuristicFPump::generateCpp( FILE * fp)
91{
92    CbcHeuristicFPump other;
93    fprintf(fp, "0#include \"CbcHeuristicFPump.hpp\"\n");
94    fprintf(fp, "3  CbcHeuristicFPump heuristicFPump(*cbcModel);\n");
95    CbcHeuristic::generateCpp(fp, "heuristicFPump");
96    if (maximumPasses_ != other.maximumPasses_)
97        fprintf(fp, "3  heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_);
98    else
99        fprintf(fp, "4  heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_);
100    if (maximumRetries_ != other.maximumRetries_)
101        fprintf(fp, "3  heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_);
102    else
103        fprintf(fp, "4  heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_);
104    if (accumulate_ != other.accumulate_)
105        fprintf(fp, "3  heuristicFPump.setAccumulate(%d);\n", accumulate_);
106    else
107        fprintf(fp, "4  heuristicFPump.setAccumulate(%d);\n", accumulate_);
108    if (fixOnReducedCosts_ != other.fixOnReducedCosts_)
109        fprintf(fp, "3  heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_);
110    else
111        fprintf(fp, "4  heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_);
112    if (maximumTime_ != other.maximumTime_)
113        fprintf(fp, "3  heuristicFPump.setMaximumTime(%g);\n", maximumTime_);
114    else
115        fprintf(fp, "4  heuristicFPump.setMaximumTime(%g);\n", maximumTime_);
116    if (fakeCutoff_ != other.fakeCutoff_)
117        fprintf(fp, "3  heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_);
118    else
119        fprintf(fp, "4  heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_);
120    if (absoluteIncrement_ != other.absoluteIncrement_)
121        fprintf(fp, "3  heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_);
122    else
123        fprintf(fp, "4  heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_);
124    if (relativeIncrement_ != other.relativeIncrement_)
125        fprintf(fp, "3  heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_);
126    else
127        fprintf(fp, "4  heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_);
128    if (defaultRounding_ != other.defaultRounding_)
129        fprintf(fp, "3  heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_);
130    else
131        fprintf(fp, "4  heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_);
132    if (initialWeight_ != other.initialWeight_)
133        fprintf(fp, "3  heuristicFPump.setInitialWeight(%g);\n", initialWeight_);
134    else
135        fprintf(fp, "4  heuristicFPump.setInitialWeight(%g);\n", initialWeight_);
136    if (weightFactor_ != other.weightFactor_)
137        fprintf(fp, "3  heuristicFPump.setWeightFactor(%g);\n", weightFactor_);
138    else
139        fprintf(fp, "4  heuristicFPump.setWeightFactor(%g);\n", weightFactor_);
140    if (artificialCost_ != other.artificialCost_)
141        fprintf(fp, "3  heuristicFPump.setArtificialCost(%g);\n", artificialCost_);
142    else
143        fprintf(fp, "4  heuristicFPump.setArtificialCost(%g);\n", artificialCost_);
144    if (iterationRatio_ != other.iterationRatio_)
145        fprintf(fp, "3  heuristicFPump.setIterationRatio(%g);\n", iterationRatio_);
146    else
147        fprintf(fp, "4  heuristicFPump.setIterationRatio(%g);\n", iterationRatio_);
148    if (reducedCostMultiplier_ != other.reducedCostMultiplier_)
149        fprintf(fp, "3  heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_);
150    else
151        fprintf(fp, "4  heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_);
152    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicFPump);\n");
153}
154
155// Copy constructor
156CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs)
157        :
158        CbcHeuristic(rhs),
159        startTime_(rhs.startTime_),
160        maximumTime_(rhs.maximumTime_),
161        fakeCutoff_(rhs.fakeCutoff_),
162        absoluteIncrement_(rhs.absoluteIncrement_),
163        relativeIncrement_(rhs.relativeIncrement_),
164        defaultRounding_(rhs.defaultRounding_),
165        initialWeight_(rhs.initialWeight_),
166        weightFactor_(rhs.weightFactor_),
167        artificialCost_(rhs.artificialCost_),
168        iterationRatio_(rhs.iterationRatio_),
169        reducedCostMultiplier_(rhs.reducedCostMultiplier_),
170        maximumPasses_(rhs.maximumPasses_),
171        maximumRetries_(rhs.maximumRetries_),
172        accumulate_(rhs.accumulate_),
173        fixOnReducedCosts_(rhs.fixOnReducedCosts_),
174        roundExpensive_(rhs.roundExpensive_)
175{
176}
177
178// Assignment operator
179CbcHeuristicFPump &
180CbcHeuristicFPump::operator=( const CbcHeuristicFPump & rhs)
181{
182    if (this != &rhs) {
183        CbcHeuristic::operator=(rhs);
184        startTime_ = rhs.startTime_;
185        maximumTime_ = rhs.maximumTime_;
186        fakeCutoff_ = rhs.fakeCutoff_;
187        absoluteIncrement_ = rhs.absoluteIncrement_;
188        relativeIncrement_ = rhs.relativeIncrement_;
189        defaultRounding_ = rhs.defaultRounding_;
190        initialWeight_ = rhs.initialWeight_;
191        weightFactor_ = rhs.weightFactor_;
192        artificialCost_ = rhs.artificialCost_;
193        iterationRatio_ = rhs.iterationRatio_;
194        reducedCostMultiplier_ = rhs.reducedCostMultiplier_;
195        maximumPasses_ = rhs.maximumPasses_;
196        maximumRetries_ = rhs.maximumRetries_;
197        accumulate_ = rhs.accumulate_;
198        fixOnReducedCosts_ = rhs.fixOnReducedCosts_;
199        roundExpensive_ = rhs.roundExpensive_;
200    }
201    return *this;
202}
203
204// Resets stuff if model changes
205void
206CbcHeuristicFPump::resetModel(CbcModel * )
207{
208}
209
210/**************************BEGIN MAIN PROCEDURE ***********************************/
211
212// See if feasibility pump will give better solution
213// Sets value of solution
214// Returns 1 if solution, 0 if not
215int
216CbcHeuristicFPump::solution(double & solutionValue,
217                            double * betterSolution)
218{
219    startTime_ = CoinCpuTime();
220    numCouldRun_++;
221    double incomingObjective = solutionValue;
222#define LEN_PRINT 200
223    char pumpPrint[LEN_PRINT];
224    pumpPrint[0] = '\0';
225/*
226  Decide if we want to run. Standard values for when are described in
227  CbcHeuristic.hpp. If we're off, or running only at root and this isn't the
228  root, bail out.
229
230  The double test (against phase, then atRoot and passNumber) has a fair bit
231  of redundancy, but the results will differ depending on whether we're
232  actually at the root of the main search tree or at the root of a small tree
233  (recursive call to branchAndBound).
234
235  FPump also supports some exotic values (11 -- 15) for when, described in
236  CbcHeuristicFPump.hpp.
237*/
238    if (!when() || (when() == 1 && model_->phase() != 1))
239        return 0; // switched off
240    // See if at root node
241    bool atRoot = model_->getNodeCount() == 0;
242    int passNumber = model_->getCurrentPassNumber();
243    // just do once
244    if (!atRoot)
245        return 0;
246    int options = feasibilityPumpOptions_;
247    if ((options % 1000000) > 0) {
248        int kOption = options / 1000000;
249        options = options % 1000000;
250        /*
251          Add 10 to do even if solution
252          1 - do after cuts
253          2 - do after cuts (not before)
254          3 - not used do after every cut round (and after cuts)
255          k not used do after every (k-2)th round
256        */
257        if (kOption < 10 && model_->getSolutionCount())
258            return 0;
259        if (model_->getSolutionCount())
260            kOption = kOption % 10;
261        bool good;
262        if (kOption == 1) {
263            good = (passNumber == 999999);
264        } else if (kOption == 2) {
265            good = (passNumber == 999999);
266            passNumber = 2; // so won't run before
267            //} else if (kOption==3) {
268            //good = true;
269        } else {
270            //good = (((passNumber-1)%(kOption-2))==0);
271            good = false;
272        }
273        if (passNumber != 1 && !good)
274            return 0;
275    } else {
276        if (passNumber != 1)
277            return 0;
278    }
279    // loop round doing repeated pumps
280    double cutoff;
281    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
282    double direction = model_->solver()->getObjSense();
283    cutoff *= direction;
284    int numberBandBsolutions = 0;
285    double firstCutoff = fabs(cutoff);
286    cutoff = CoinMin(cutoff, solutionValue);
287    // check plausible and space for rounded solution
288    int numberColumns = model_->getNumCols();
289    int numberIntegers = model_->numberIntegers();
290    const int * integerVariableOrig = model_->integerVariable();
291    double iterationLimit = -1.0;
292    //iterationRatio_=1.0;
293    if (iterationRatio_ > 0.0)
294        iterationLimit = (2 * model_->solver()->getNumRows() + 2 * numberColumns) *
295                         iterationRatio_;
296    int totalNumberIterations = 0;
297    int averageIterationsPerTry = -1;
298    int numberIterationsLastPass = 0;
299    // 1. initially check 0-1
300/*
301  I'm skeptical of the above comment, but it's likely accurate as the default.
302  Bit 4 or bit 8 needs to be set in order to consider working with general
303  integers.
304*/
305    int i, j;
306    int general = 0;
307    int * integerVariable = new int[numberIntegers];
308    const double * lower = model_->solver()->getColLower();
309    const double * upper = model_->solver()->getColUpper();
310    bool doGeneral = (accumulate_ & 4) != 0;
311    j = 0;
312/*
313  Scan the objects, recording the columns and counting general integers.
314
315  Seems like the NDEBUG tests could be made into an applicability test. If
316  a scan of the objects reveals complex objects, just clean up and return
317  failure.
318*/
319    for (i = 0; i < numberIntegers; i++) {
320        int iColumn = integerVariableOrig[i];
321#ifndef NDEBUG
322        const OsiObject * object = model_->object(i);
323        const CbcSimpleInteger * integerObject =
324            dynamic_cast<const  CbcSimpleInteger *> (object);
325        const OsiSimpleInteger * integerObject2 =
326            dynamic_cast<const  OsiSimpleInteger *> (object);
327        assert(integerObject || integerObject2);
328#endif
329        if (upper[iColumn] - lower[iColumn] > 1.000001) {
330            general++;
331            if (doGeneral)
332                integerVariable[j++] = iColumn;
333        } else {
334            integerVariable[j++] = iColumn;
335        }
336    }
337/*
338  If 2/3 of integers are general integers, and we're not going to work with
339  them, might as well go home.
340
341  The else case is unclear to me. We reach it if general integers are less than
342  2/3 of the total, or if either of bit 4 or 8 is set. But only bit 8 is used
343  in the decision. (Let manyGen = 1 if more than 2/3 of integers are general
344  integers. Then a k-map on manyGen, bit4, and bit8 shows it clearly.)
345
346  So there's something odd here. In the case where bit4 = 1 and bit8 = 0,
347  we've included general integers in integerVariable, but we're not going to
348  process them.
349*/
350    if (general*3 > 2*numberIntegers && !doGeneral) {
351        delete [] integerVariable;
352        return 0;
353    } else if ((accumulate_&4) == 0) {
354        doGeneral = false;
355        j = 0;
356        for (i = 0; i < numberIntegers; i++) {
357            int iColumn = integerVariableOrig[i];
358            if (upper[iColumn] - lower[iColumn] < 1.000001)
359                integerVariable[j++] = iColumn;
360        }
361    }
362    if (!general)
363        doGeneral = false;
364#ifdef CLP_INVESTIGATE
365    if (doGeneral)
366        printf("DOing general with %d out of %d\n", general, numberIntegers);
367#endif
368/*
369  This `closest solution' will satisfy integrality, but violate some other
370  constraints?
371*/
372    // For solution closest to feasible if none found
373    int * closestSolution = general ? NULL : new int[numberIntegers];
374    double closestObjectiveValue = COIN_DBL_MAX;
375
376    int numberIntegersOrig = numberIntegers;
377    numberIntegers = j;
378    double * newSolution = new double [numberColumns];
379    double newSolutionValue = COIN_DBL_MAX;
380    int maxSolutions = model_->getMaximumSolutions();
381    int numberSolutions=0;
382    bool solutionFound = false;
383    int * usedColumn = NULL;
384    double * lastSolution = NULL;
385    int fixContinuous = 0;
386    bool fixInternal = false;
387    bool allSlack = false;
388    if (when_ >= 21 && when_ <= 25) {
389        when_ -= 10;
390        allSlack = true;
391    }
392    double time1 = CoinCpuTime();
393/*
394  Obtain a relaxed lp solution.
395*/
396    model_->solver()->resolve();
397    if (!model_->solver()->isProvenOptimal()) {
398        // presumably max time or some such
399        return 0;
400    }
401    numRuns_++;
402    if (cutoff < 1.0e50 && false) {
403        // Fix on djs
404        double direction = model_->solver()->getObjSense() ;
405        double gap = cutoff - model_->solver()->getObjValue() * direction ;
406        double tolerance;
407        model_->solver()->getDblParam(OsiDualTolerance, tolerance) ;
408        if (gap > 0.0) {
409            gap += 100.0 * tolerance;
410            int nFix = model_->solver()->reducedCostFix(gap);
411            printf("dj fixing fixed %d variables\n", nFix);
412        }
413    }
414/*
415  I have no idea why we're doing this, except perhaps that saveBasis will be
416  automagically deleted on exit from the routine.
417*/
418    CoinWarmStartBasis saveBasis;
419    CoinWarmStartBasis * basis =
420        dynamic_cast<CoinWarmStartBasis *>(model_->solver()->getWarmStart()) ;
421    if (basis) {
422        saveBasis = * basis;
423        delete basis;
424    }
425    double continuousObjectiveValue = model_->solver()->getObjValue() * model_->solver()->getObjSense();
426    double * firstPerturbedObjective = NULL;
427    double * firstPerturbedSolution = NULL;
428    double firstPerturbedValue = COIN_DBL_MAX;
429    if (when_ >= 11 && when_ <= 15) {
430        fixInternal = when_ > 11 && when_ < 15;
431        if (when_ < 13)
432            fixContinuous = 0;
433        else if (when_ != 14)
434            fixContinuous = 1;
435        else
436            fixContinuous = 2;
437        when_ = 1;
438        if ((accumulate_&1) != 0) {
439            usedColumn = new int [numberColumns];
440            for (int i = 0; i < numberColumns; i++)
441                usedColumn[i] = -1;
442        }
443        lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(), numberColumns);
444    }
445    int finalReturnCode = 0;
446    int totalNumberPasses = 0;
447    int numberTries = 0;
448    CoinWarmStartBasis bestBasis;
449    bool exitAll = false;
450    //double saveBestObjective = model_->getMinimizationObjValue();
451    OsiSolverInterface * solver = NULL;
452    double artificialFactor = 0.00001;
453    // also try rounding!
454    double * roundingSolution = new double[numberColumns];
455    double roundingObjective = solutionValue;
456    CbcRounding roundingHeuristic(*model_);
457    int dualPass = 0;
458    int secondPassOpt = 0;
459#define RAND_RAND
460#ifdef RAND_RAND
461    int offRandom = 0;
462#endif
463    int maximumAllowed = -1;
464    bool moreIterations = false;
465    if (options > 0) {
466        if (options >= 1000)
467            maximumAllowed = options / 1000;
468        int options2 = (options % 1000) / 100;
469#ifdef RAND_RAND
470        offRandom = options2 & 1;
471#endif
472        moreIterations = (options2 & 2) != 0;
473        secondPassOpt = (options / 10) % 10;
474        /* 1 to 7 - re-use solution
475           8 use dual and current solution(ish)
476           9 use dual and allslack
477           1 - primal and mod obj
478           2 - dual and mod obj
479           3 - primal and no mod obj
480           add 3 to redo current solution
481        */
482        if (secondPassOpt >= 8) {
483            dualPass = secondPassOpt - 7;
484            secondPassOpt = 0;
485        }
486    }
487    // Number of passes to do
488    int maximumPasses = maximumPasses_;
489#ifdef COIN_HAS_CLP
490    if (maximumPasses == 30) {
491        OsiClpSolverInterface * clpSolver
492        = dynamic_cast<OsiClpSolverInterface *> (model_->solver());
493        if (clpSolver && clpSolver->fakeObjective())
494            maximumPasses = 100; // feasibility problem?
495    }
496#endif
497#ifdef RAND_RAND
498    double * randomFactor = new double [numberColumns];
499    for (int i = 0; i < numberColumns; i++) {
500        double value = floor(1.0e3 * randomNumberGenerator_.randomDouble());
501        randomFactor[i] = 1.0 + value * 1.0e-4;
502    }
503#endif
504    // guess exact multiple of objective
505    double exactMultiple = model_->getCutoffIncrement();
506    exactMultiple *= 2520;
507    if (fabs(exactMultiple / 0.999 - floor(exactMultiple / 0.999 + 0.5)) < 1.0e-9)
508        exactMultiple /= 2520.0 * 0.999;
509    else if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-9)
510        exactMultiple /= 2520.0;
511    else
512        exactMultiple = 0.0;
513    // check for rounding errors (only for integral case)
514    if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-8)
515        exactMultiple = floor(exactMultiple + 0.5);
516    //printf("exact multiple %g\n",exactMultiple);
517    // Clone solver for rounding
518    OsiSolverInterface * clonedSolver = cloneBut(2); // wasmodel_->solver()->clone();
519    while (!exitAll) {
520        // Cutoff rhs
521        double useRhs = COIN_DBL_MAX;
522        double useOffset = 0.0;
523        int numberPasses = 0;
524        artificialFactor *= 10.0;
525        int lastMove = (!numberTries) ? -10 : 1000000;
526        double lastSumInfeas = COIN_DBL_MAX;
527        numberTries++;
528        // Clone solver - otherwise annoys root node computations
529        solver = cloneBut(2); // was model_->solver()->clone();
530#ifdef COIN_HAS_CLP
531        {
532            OsiClpSolverInterface * clpSolver
533            = dynamic_cast<OsiClpSolverInterface *> (solver);
534            if (clpSolver) {
535                // better to clean up using primal?
536                ClpSimplex * lp = clpSolver->getModelPtr();
537                int options = lp->specialOptions();
538                lp->setSpecialOptions(options | 8192);
539                //lp->setSpecialOptions(options|0x01000000);
540#ifdef CLP_INVESTIGATE
541                clpSolver->setHintParam(OsiDoReducePrint, false, OsiHintTry);
542                lp->setLogLevel(CoinMax(1, lp->logLevel()));
543#endif
544            }
545        }
546#endif
547        if (CoinMin(fakeCutoff_, cutoff) < 1.0e50) {
548            // Fix on djs
549            double direction = solver->getObjSense() ;
550            double gap = CoinMin(fakeCutoff_, cutoff) - solver->getObjValue() * direction ;
551            double tolerance;
552            solver->getDblParam(OsiDualTolerance, tolerance) ;
553            if (gap > 0.0 && (fixOnReducedCosts_ == 1 || (numberTries == 1 && fixOnReducedCosts_ == 2))) {
554                gap += 100.0 * tolerance;
555                gap *= reducedCostMultiplier_;
556                int nFix = solver->reducedCostFix(gap);
557                if (nFix) {
558                    sprintf(pumpPrint, "Reduced cost fixing fixed %d variables on major pass %d", nFix, numberTries);
559                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
560                    << pumpPrint
561                    << CoinMessageEol;
562                    //pumpPrint[0]='\0';
563                }
564            }
565        }
566        // if cutoff exists then add constraint
567        bool useCutoff = (fabs(cutoff) < 1.0e20 && (fakeCutoff_ != COIN_DBL_MAX || numberTries > 1));
568        // but there may be a close one
569        if (firstCutoff < 2.0*solutionValue && numberTries == 1 && CoinMin(cutoff, fakeCutoff_) < 1.0e20)
570            useCutoff = true;
571        if (useCutoff) {
572            double rhs = CoinMin(cutoff, fakeCutoff_);
573            const double * objective = solver->getObjCoefficients();
574            int numberColumns = solver->getNumCols();
575            int * which = new int[numberColumns];
576            double * els = new double[numberColumns];
577            int nel = 0;
578            for (int i = 0; i < numberColumns; i++) {
579                double value = objective[i];
580                if (value) {
581                    which[nel] = i;
582                    els[nel++] = direction * value;
583                }
584            }
585            solver->getDblParam(OsiObjOffset, useOffset);
586#ifdef COIN_DEVELOP
587            if (useOffset)
588                printf("CbcHeuristicFPump obj offset %g\n", useOffset);
589#endif
590            useOffset *= direction;
591            // Tweak rhs and save
592            useRhs = rhs;
593#ifdef JJF_ZERO
594            double tempValue = 60.0 * useRhs;
595            if (fabs(tempValue - floor(tempValue + 0.5)) < 1.0e-7 && rhs != fakeCutoff_) {
596                // add a little
597                useRhs += 1.0e-5;
598            }
599#endif
600            solver->addRow(nel, which, els, -COIN_DBL_MAX, useRhs + useOffset);
601            delete [] which;
602            delete [] els;
603            bool takeHint;
604            OsiHintStrength strength;
605            solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
606            solver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
607            solver->resolve();
608            solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
609            if (!solver->isProvenOptimal()) {
610                // presumably max time or some such
611                exitAll = true;
612                break;
613            }
614        }
615        solver->setDblParam(OsiDualObjectiveLimit, 1.0e50);
616        solver->resolve();
617        // Solver may not be feasible
618        if (!solver->isProvenOptimal()) {
619            exitAll = true;
620            break;
621        }
622        const double * lower = solver->getColLower();
623        const double * upper = solver->getColUpper();
624        const double * solution = solver->getColSolution();
625        if (lastSolution)
626            memcpy(lastSolution, solution, numberColumns*sizeof(double));
627        double primalTolerance;
628        solver->getDblParam(OsiPrimalTolerance, primalTolerance);
629
630        // 2 space for last rounded solutions
631#define NUMBER_OLD 4
632        double ** oldSolution = new double * [NUMBER_OLD];
633        for (j = 0; j < NUMBER_OLD; j++) {
634            oldSolution[j] = new double[numberColumns];
635            for (i = 0; i < numberColumns; i++) oldSolution[j][i] = -COIN_DBL_MAX;
636        }
637
638        // 3. Replace objective with an initial 0-valued objective
639        double * saveObjective = new double [numberColumns];
640        memcpy(saveObjective, solver->getObjCoefficients(), numberColumns*sizeof(double));
641        for (i = 0; i < numberColumns; i++) {
642            solver->setObjCoeff(i, 0.0);
643        }
644        bool finished = false;
645        double direction = solver->getObjSense();
646        int returnCode = 0;
647        bool takeHint;
648        OsiHintStrength strength;
649        solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
650        solver->setHintParam(OsiDoDualInResolve, false);
651        //solver->messageHandler()->setLogLevel(0);
652
653        // 4. Save objective offset so we can see progress
654        double saveOffset;
655        solver->getDblParam(OsiObjOffset, saveOffset);
656        // Get amount for original objective
657        double scaleFactor = 0.0;
658#ifdef COIN_DEVELOP
659        double largestCost = 0.0;
660        int nArtificial = 0;
661#endif
662        for (i = 0; i < numberColumns; i++) {
663            double value = saveObjective[i];
664            scaleFactor += value * value;
665#ifdef COIN_DEVELOP
666            largestCost = CoinMax(largestCost, fabs(value));
667            if (value*direction >= artificialCost_)
668                nArtificial++;
669#endif
670        }
671        if (scaleFactor)
672            scaleFactor = (initialWeight_ * sqrt(static_cast<double> (numberIntegers))) / sqrt(scaleFactor);
673#ifdef CLP_INVESTIGATE
674#ifdef COIN_DEVELOP
675        if (scaleFactor || nArtificial)
676            printf("Using %g fraction of original objective (decay %g) - largest %g - %d artificials\n", scaleFactor, weightFactor_,
677                   largestCost, nArtificial);
678#else
679        if (scaleFactor)
680            printf("Using %g fraction of original objective (decay %g)\n",
681                   scaleFactor, weightFactor_);
682#endif
683#endif
684        // This is an array of sums of infeasibilities so can see if "bobbling"
685#define SIZE_BOBBLE 20
686        double saveSumInf[SIZE_BOBBLE];
687        CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
688        // 0 before doing anything
689        int bobbleMode = 0;
690        // 5. MAIN WHILE LOOP
691        //bool newLineNeeded=false;
692/*
693  finished occurs exactly twice in this routine: immediately above, where it's
694  set to false, and here in the loop condition.
695*/
696        while (!finished) {
697            double newTrueSolutionValue = 0.0;
698            double newSumInfeas = 0.0;
699            int newNumberInfeas = 0;
700            returnCode = 0;
701            if (model_->maximumSecondsReached()) {
702                exitAll = true;
703                break;
704            }
705            // see what changed
706            if (usedColumn) {
707                for (i = 0; i < numberColumns; i++) {
708                    if (fabs(solution[i] - lastSolution[i]) > 1.0e-8)
709                        usedColumn[i] = numberPasses;
710                    lastSolution[i] = solution[i];
711                }
712            }
713            if (averageIterationsPerTry >= 0) {
714                int n = totalNumberIterations - numberIterationsLastPass;
715                double perPass = totalNumberIterations /
716                                 (totalNumberPasses + numberPasses + 1.0e-5);
717                perPass /= (solver->getNumRows() + numberColumns);
718                double test = moreIterations ? 0.3 : 0.05;
719                if (n > CoinMax(20000, 3*averageIterationsPerTry)
720                        && (switches_&2) == 0 && maximumPasses<200 && perPass>test) {
721                    exitAll = true;
722                }
723            }
724            // Exit on exact total number if maximumPasses large
725            if ((maximumPasses >= 200 || (switches_&2) != 0)
726                    && numberPasses + totalNumberPasses >=
727                    maximumPasses)
728                exitAll = true;
729            bool exitThis = false;
730            if (iterationLimit < 0.0) {
731                if (numberPasses >= maximumPasses) {
732                    // If going well then keep going if maximumPasses small
733                    if (lastMove < numberPasses - 4 || lastMove == 1000000)
734                        exitThis = true;
735                    if (maximumPasses > 20 || numberPasses >= 40)
736                        exitThis = true;
737                }
738            }
739            if (iterationLimit > 0.0 && totalNumberIterations > iterationLimit
740                    && numberPasses > 15) {
741                // exiting on iteration count
742                exitAll = true;
743            } else if (maximumPasses<30 && numberPasses>100) {
744                // too many passes anyway
745                exitAll = true;
746            }
747            if (maximumTime_ > 0.0 && CoinCpuTime() >= startTime_ + maximumTime_) {
748                exitAll = true;
749                // force exit
750                switches_ |= 2048;
751            }
752            if (exitAll || exitThis)
753                break;
754            memcpy(newSolution, solution, numberColumns*sizeof(double));
755            int flip;
756            if (numberPasses == 0 && false) {
757                // always use same seed
758                randomNumberGenerator_.setSeed(987654321);
759            }
760            returnCode = rounds(solver, newSolution,/*saveObjective,*/
761                                numberIntegers, integerVariable,
762                                /*pumpPrint,*/numberPasses,
763                                /*roundExpensive_,*/defaultRounding_, &flip);
764            if (numberPasses == 0 && false) {
765                // Make sure random will be different
766                for (i = 1; i < numberTries; i++)
767                    randomNumberGenerator_.randomDouble();
768            }
769            numberPasses++;
770            if (returnCode) {
771                // SOLUTION IS INTEGER
772                // Put back correct objective
773                for (i = 0; i < numberColumns; i++)
774                    solver->setObjCoeff(i, saveObjective[i]);
775
776                // solution - but may not be better
777                // Compute using dot product
778                solver->setDblParam(OsiObjOffset, saveOffset);
779                newSolutionValue = -saveOffset;
780                for (  i = 0 ; i < numberColumns ; i++ )
781                    newSolutionValue += saveObjective[i] * newSolution[i];
782                newSolutionValue *= direction;
783                sprintf(pumpPrint, "Solution found of %g", newSolutionValue);
784                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
785                << pumpPrint
786                << CoinMessageEol;
787                //newLineNeeded=false;
788                if (newSolutionValue < solutionValue) {
789                    double saveValue = solutionValue;
790                    if (!doGeneral) {
791                        int numberLeft = 0;
792                        for (i = 0; i < numberIntegersOrig; i++) {
793                            int iColumn = integerVariableOrig[i];
794                            double value = floor(newSolution[iColumn] + 0.5);
795                            if (solver->isBinary(iColumn)) {
796                                solver->setColLower(iColumn, value);
797                                solver->setColUpper(iColumn, value);
798                            } else {
799                                if (fabs(value - newSolution[iColumn]) > 1.0e-7)
800                                    numberLeft++;
801                            }
802                        }
803                        if (numberLeft) {
804                            sprintf(pumpPrint, "Branch and bound needed to clear up %d general integers", numberLeft);
805                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
806                            << pumpPrint
807                            << CoinMessageEol;
808                            returnCode = smallBranchAndBound(solver, numberNodes_, newSolution, newSolutionValue,
809                                                             solutionValue, "CbcHeuristicFpump");
810                            if (returnCode < 0) {
811                                if (returnCode == -2)
812                                    exitAll = true;
813                                returnCode = 0; // returned on size or event
814                            }
815                            if ((returnCode&2) != 0) {
816                                // could add cut
817                                returnCode &= ~2;
818                            }
819                            if (returnCode != 1)
820                                newSolutionValue = saveValue;
821                            if (returnCode && newSolutionValue < saveValue)
822                                numberBandBsolutions++;
823                        }
824                    }
825                    if (returnCode && newSolutionValue < saveValue) {
826                        memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
827                        solutionFound = true;
828                        if (exitNow(newSolutionValue))
829                            exitAll = true;
830                        CoinWarmStartBasis * basis =
831                            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
832                        if (basis) {
833                            bestBasis = * basis;
834                            delete basis;
835                            int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
836                            if (action == 0) {
837                                double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
838                                double saveObjectiveValue = model_->getMinimizationObjValue();
839                                model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
840                                if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
841                                    model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
842                                delete [] saveOldSolution;
843                            }
844                            if (action == 0 || model_->maximumSecondsReached()) {
845                                exitAll = true; // exit
846                                break;
847                            }
848                        }
849                        if ((accumulate_&1) != 0) {
850                            model_->incrementUsed(betterSolution); // for local search
851                        }
852                        solutionValue = newSolutionValue;
853                        solutionFound = true;
854                        numberSolutions++;
855                        if (numberSolutions>=maxSolutions)
856                          exitAll = true;
857                        if (general && saveValue != newSolutionValue) {
858                            sprintf(pumpPrint, "Cleaned solution of %g", solutionValue);
859                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
860                            << pumpPrint
861                            << CoinMessageEol;
862                        }
863                        if (exitNow(newSolutionValue))
864                            exitAll = true;
865                    } else {
866                        sprintf(pumpPrint, "Mini branch and bound could not fix general integers");
867                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
868                        << pumpPrint
869                        << CoinMessageEol;
870                    }
871                } else {
872                    sprintf(pumpPrint, "After further testing solution no better than previous of %g", solutionValue);
873                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
874                    << pumpPrint
875                    << CoinMessageEol;
876                    //newLineNeeded=false;
877                    returnCode = 0;
878                }
879                break;
880            } else {
881                // SOLUTION IS not INTEGER
882                // 1. check for loop
883                bool matched;
884                for (int k = NUMBER_OLD - 1; k > 0; k--) {
885                    double * b = oldSolution[k];
886                    matched = true;
887                    for (i = 0; i < numberIntegers; i++) {
888                        int iColumn = integerVariable[i];
889                        if (newSolution[iColumn] != b[iColumn]) {
890                            matched = false;
891                            break;
892                        }
893                    }
894                    if (matched) break;
895                }
896                int numberPerturbed = 0;
897                if (matched || numberPasses % 100 == 0) {
898                    // perturbation
899                    //sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
900                    //newLineNeeded=true;
901                    double factorX[10] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
902                    double factor = 1.0;
903                    double target = -1.0;
904                    double * randomX = new double [numberIntegers];
905                    for (i = 0; i < numberIntegers; i++)
906                        randomX[i] = CoinMax(0.0, randomNumberGenerator_.randomDouble() - 0.3);
907                    for (int k = 0; k < 10; k++) {
908#ifdef COIN_DEVELOP_x
909                        printf("kpass %d\n", k);
910#endif
911                        int numberX[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
912                        for (i = 0; i < numberIntegers; i++) {
913                            int iColumn = integerVariable[i];
914                            double value = randomX[i];
915                            double difference = fabs(solution[iColumn] - newSolution[iColumn]);
916                            for (int j = 0; j < 10; j++) {
917                                if (difference + value*factorX[j] > 0.5)
918                                    numberX[j]++;
919                            }
920                        }
921                        if (target < 0.0) {
922                            if (numberX[9] <= 200)
923                                break; // not very many changes
924                            target = CoinMax(200.0, CoinMin(0.05 * numberX[9], 1000.0));
925                        }
926                        int iX = -1;
927                        int iBand = -1;
928                        for (i = 0; i < 10; i++) {
929#ifdef COIN_DEVELOP_x
930                            printf("** %d changed at %g\n", numberX[i], factorX[i]);
931#endif
932                            if (numberX[i] >= target && numberX[i] < 2.0*target && iX < 0)
933                                iX = i;
934                            if (iBand<0 && numberX[i]>target) {
935                                iBand = i;
936                                factor = factorX[i];
937                            }
938                        }
939                        if (iX >= 0) {
940                            factor = factorX[iX];
941                            break;
942                        } else {
943                            assert (iBand >= 0);
944                            double hi = factor;
945                            double lo = (iBand > 0) ? factorX[iBand-1] : 0.0;
946                            double diff = (hi - lo) / 9.0;
947                            for (i = 0; i < 10; i++) {
948                                factorX[i] = lo;
949                                lo += diff;
950                            }
951                        }
952                    }
953                    for (i = 0; i < numberIntegers; i++) {
954                        int iColumn = integerVariable[i];
955                        double value = randomX[i];
956                        double difference = fabs(solution[iColumn] - newSolution[iColumn]);
957                        if (difference + value*factor > 0.5) {
958                            numberPerturbed++;
959                            if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
960                                newSolution[iColumn] += 1.0;
961                            } else if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
962                                newSolution[iColumn] -= 1.0;
963                            } else {
964                                // general integer
965                                if (difference + value > 0.75)
966                                    newSolution[iColumn] += 1.0;
967                                else
968                                    newSolution[iColumn] -= 1.0;
969                            }
970                        }
971                    }
972                    delete [] randomX;
973                } else {
974                    for (j = NUMBER_OLD - 1; j > 0; j--) {
975                        for (i = 0; i < numberColumns; i++) oldSolution[j][i] = oldSolution[j-1][i];
976                    }
977                    for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
978                }
979
980                // 2. update the objective function based on the new rounded solution
981                double offset = 0.0;
982                double costValue = (1.0 - scaleFactor) * solver->getObjSense();
983                int numberChanged = 0;
984                const double * oldObjective = solver->getObjCoefficients();
985                for (i = 0; i < numberColumns; i++) {
986                    // below so we can keep original code and allow for objective
987                    int iColumn = i;
988                    // Special code for "artificials"
989                    if (direction*saveObjective[iColumn] >= artificialCost_) {
990                        //solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]);
991                        solver->setObjCoeff(iColumn, (artificialFactor*saveObjective[iColumn]) / artificialCost_);
992                    }
993                    if (!solver->isBinary(iColumn) && !doGeneral)
994                        continue;
995                    // deal with fixed variables (i.e., upper=lower)
996                    if (fabs(lower[iColumn] - upper[iColumn]) < primalTolerance || !solver->isInteger(iColumn)) {
997                        //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
998                        //else                                       solver->setObjCoeff(iColumn,costValue);
999                        continue;
1000                    }
1001                    double newValue = 0.0;
1002                    if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
1003                        newValue = costValue + scaleFactor * saveObjective[iColumn];
1004                    } else {
1005                        if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
1006                            newValue = -costValue + scaleFactor * saveObjective[iColumn];
1007                        }
1008                    }
1009#ifdef RAND_RAND
1010                    if (!offRandom)
1011                        newValue *= randomFactor[iColumn];
1012#endif
1013                    if (newValue != oldObjective[iColumn]) {
1014                        numberChanged++;
1015                    }
1016                    solver->setObjCoeff(iColumn, newValue);
1017                    offset += costValue * newSolution[iColumn];
1018                }
1019                solver->setDblParam(OsiObjOffset, -offset);
1020                if (!general && false) {
1021                    // Solve in two goes - first keep satisfied ones fixed
1022                    double * saveLower = new double [numberIntegers];
1023                    double * saveUpper = new double [numberIntegers];
1024                    for (i = 0; i < numberIntegers; i++) {
1025                        int iColumn = integerVariable[i];
1026                        saveLower[i] = COIN_DBL_MAX;
1027                        saveUpper[i] = -COIN_DBL_MAX;
1028                        if (solution[iColumn] < lower[iColumn] + primalTolerance) {
1029                            saveUpper[i] = upper[iColumn];
1030                            solver->setColUpper(iColumn, lower[iColumn]);
1031                        } else if (solution[iColumn] > upper[iColumn] - primalTolerance) {
1032                            saveLower[i] = lower[iColumn];
1033                            solver->setColLower(iColumn, upper[iColumn]);
1034                        }
1035                    }
1036                    solver->resolve();
1037                    if (!solver->isProvenOptimal()) {
1038                        // presumably max time or some such
1039                        exitAll = true;
1040                        break;
1041                    }
1042                    for (i = 0; i < numberIntegers; i++) {
1043                        int iColumn = integerVariable[i];
1044                        if (saveLower[i] != COIN_DBL_MAX)
1045                            solver->setColLower(iColumn, saveLower[i]);
1046                        if (saveUpper[i] != -COIN_DBL_MAX)
1047                            solver->setColUpper(iColumn, saveUpper[i]);
1048                        saveUpper[i] = -COIN_DBL_MAX;
1049                    }
1050                    memcpy(newSolution, solution, numberColumns*sizeof(double));
1051                    int flip;
1052                    returnCode = rounds(solver, newSolution,/*saveObjective,*/
1053                                        numberIntegers, integerVariable,
1054                                        /*pumpPrint,*/numberPasses,
1055                                        /*roundExpensive_,*/defaultRounding_, &flip);
1056                    numberPasses++;
1057                    if (returnCode) {
1058                        // solution - but may not be better
1059                        // Compute using dot product
1060                        double newSolutionValue = -saveOffset;
1061                        for (  i = 0 ; i < numberColumns ; i++ )
1062                            newSolutionValue += saveObjective[i] * newSolution[i];
1063                        newSolutionValue *= direction;
1064                        sprintf(pumpPrint, "Intermediate solution found of %g", newSolutionValue);
1065                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1066                        << pumpPrint
1067                        << CoinMessageEol;
1068                        if (newSolutionValue < solutionValue) {
1069                            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1070                            CoinWarmStartBasis * basis =
1071                                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1072                            solutionFound = true;
1073                            numberSolutions++;
1074                            if (numberSolutions>=maxSolutions)
1075                              exitAll = true;
1076                            if (exitNow(newSolutionValue))
1077                                exitAll = true;
1078                            if (basis) {
1079                                bestBasis = * basis;
1080                                delete basis;
1081                                int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
1082                                if (!action) {
1083                                    double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
1084                                    double saveObjectiveValue = model_->getMinimizationObjValue();
1085                                    model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
1086                                    if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
1087                                        model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
1088                                    delete [] saveOldSolution;
1089                                }
1090                                if (!action || model_->maximumSecondsReached()) {
1091                                    exitAll = true; // exit
1092                                    break;
1093                                }
1094                            }
1095                            if ((accumulate_&1) != 0) {
1096                                model_->incrementUsed(betterSolution); // for local search
1097                            }
1098                            solutionValue = newSolutionValue;
1099                            solutionFound = true;
1100                            numberSolutions++;
1101                            if (numberSolutions>=maxSolutions)
1102                              exitAll = true;
1103                            if (exitNow(newSolutionValue))
1104                                exitAll = true;
1105                        } else {
1106                            returnCode = 0;
1107                        }
1108                    }
1109                }
1110                int numberIterations = 0;
1111                if (!doGeneral) {
1112                    // faster to do from all slack!!!!
1113                    if (allSlack) {
1114                        CoinWarmStartBasis dummy;
1115                        solver->setWarmStart(&dummy);
1116                    }
1117#ifdef COIN_DEVELOP
1118                    printf("%d perturbed out of %d columns (%d changed)\n", numberPerturbed, numberColumns, numberChanged);
1119#endif
1120                    bool takeHint;
1121                    OsiHintStrength strength;
1122                    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
1123                    if (dualPass && numberChanged > 2) {
1124                        solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
1125                        if (dualPass == 1 && 2*numberChanged < numberColumns &&
1126                                (numberChanged < 5000 || 6*numberChanged < numberColumns)) {
1127                            // but we need to make infeasible
1128                            CoinWarmStartBasis * basis =
1129                                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1130                            if (basis) {
1131                                // modify
1132                                const double * lower = solver->getColLower();
1133                                const double * upper = solver->getColUpper();
1134                                double * solution = CoinCopyOfArray(solver->getColSolution(),
1135                                                                    numberColumns);
1136                                const double * objective = solver->getObjCoefficients();
1137                                int nChanged = 0;
1138                                for (i = 0; i < numberIntegersOrig; i++) {
1139                                    int iColumn = integerVariableOrig[i];
1140#ifdef RAND_RAND
1141                                    if (nChanged > numberChanged)
1142                                        break;
1143#endif
1144                                    if (objective[iColumn] > 0.0) {
1145                                        if (basis->getStructStatus(iColumn) ==
1146                                                CoinWarmStartBasis::atUpperBound) {
1147                                            solution[iColumn] = lower[iColumn];
1148                                            basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
1149                                            nChanged++;
1150                                        }
1151                                    } else if (objective[iColumn] < 0.0) {
1152                                        if (basis->getStructStatus(iColumn) ==
1153                                                CoinWarmStartBasis::atLowerBound) {
1154                                            solution[iColumn] = upper[iColumn];
1155                                            basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
1156                                            nChanged++;
1157                                        }
1158                                    }
1159                                }
1160                                if (!nChanged) {
1161                                    for (i = 0; i < numberIntegersOrig; i++) {
1162                                        int iColumn = integerVariableOrig[i];
1163                                        if (objective[iColumn] > 0.0) {
1164                                            if (basis->getStructStatus(iColumn) ==
1165                                                    CoinWarmStartBasis::basic) {
1166                                                solution[iColumn] = lower[iColumn];
1167                                                basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
1168                                                break;
1169                                            }
1170                                        } else if (objective[iColumn] < 0.0) {
1171                                            if (basis->getStructStatus(iColumn) ==
1172                                                    CoinWarmStartBasis::basic) {
1173                                                solution[iColumn] = upper[iColumn];
1174                                                basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
1175                                                break;
1176                                            }
1177                                        }
1178                                    }
1179                                }
1180                                solver->setColSolution(solution);
1181                                delete [] solution;
1182                                solver->setWarmStart(basis);
1183                                delete basis;
1184                            }
1185                        } else {
1186                            // faster to do from all slack!!!! ???
1187                            CoinWarmStartBasis dummy;
1188                            solver->setWarmStart(&dummy);
1189                        }
1190                    }
1191                    if (numberTries > 1 && numberPasses == 1 && firstPerturbedObjective) {
1192                        // Modify to use convex combination
1193                        // use basis from first time
1194                        solver->setWarmStart(&saveBasis);
1195                        // and objective
1196                        if (secondPassOpt < 3 || (secondPassOpt >= 4 && secondPassOpt < 6))
1197                            solver->setObjective(firstPerturbedObjective);
1198                        // and solution
1199                        solver->setColSolution(firstPerturbedSolution);
1200                        //if (secondPassOpt==2||secondPassOpt==5||
1201                        if (firstPerturbedValue > cutoff)
1202                            solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
1203                    }
1204                    solver->resolve();
1205                    if (!solver->isProvenOptimal()) {
1206                        // presumably max time or some such
1207                        exitAll = true;
1208                        break;
1209                    }
1210                    solver->setHintParam(OsiDoDualInResolve, takeHint);
1211                    newTrueSolutionValue = -saveOffset;
1212                    newSumInfeas = 0.0;
1213                    newNumberInfeas = 0;
1214                    {
1215                        const double * newSolution = solver->getColSolution();
1216                        for (  i = 0 ; i < numberColumns ; i++ ) {
1217                            if (solver->isInteger(i)) {
1218                                double value = newSolution[i];
1219                                double nearest = floor(value + 0.5);
1220                                newSumInfeas += fabs(value - nearest);
1221                                if (fabs(value - nearest) > 1.0e-6)
1222                                    newNumberInfeas++;
1223                            }
1224                            newTrueSolutionValue += saveObjective[i] * newSolution[i];
1225                        }
1226                        newTrueSolutionValue *= direction;
1227                        if (numberPasses == 1 && secondPassOpt) {
1228                            if (numberTries == 1 || secondPassOpt > 3) {
1229                                // save basis
1230                                CoinWarmStartBasis * basis =
1231                                    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1232                                if (basis) {
1233                                    saveBasis = * basis;
1234                                    delete basis;
1235                                }
1236                                delete [] firstPerturbedObjective;
1237                                delete [] firstPerturbedSolution;
1238                                firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
1239                                firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(), numberColumns);
1240                                firstPerturbedValue = newTrueSolutionValue;
1241                            }
1242                        }
1243                        if (newNumberInfeas && newNumberInfeas < 15) {
1244#ifdef JJF_ZERO
1245                            roundingObjective = solutionValue;
1246                            OsiSolverInterface * saveSolver = model_->swapSolver(solver);
1247                            double * currentObjective =
1248                                CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
1249                            solver->setObjective(saveObjective);
1250                            double saveOffset2;
1251                            solver->getDblParam(OsiObjOffset, saveOffset2);
1252                            solver->setDblParam(OsiObjOffset, saveOffset);
1253                            int ifSol = roundingHeuristic.solution(roundingObjective, roundingSolution);
1254                            solver->setObjective(currentObjective);
1255                            solver->setDblParam(OsiObjOffset, saveOffset2);
1256                            delete [] currentObjective;
1257                            model_->swapSolver(saveSolver);
1258                            if (ifSol > 0)
1259                                abort();
1260#endif
1261                            int numberRows = solver->getNumRows();
1262                            double * rowActivity = new double[numberRows];
1263                            memset(rowActivity, 0, numberRows*sizeof(double));
1264                            int * which = new int[newNumberInfeas];
1265                            int * stack = new int[newNumberInfeas+1];
1266                            double * baseValue = new double[newNumberInfeas];
1267                            int * whichRow = new int[numberRows];
1268                            double * rowValue = new double[numberRows];
1269                            memset(rowValue, 0, numberRows*sizeof(double));
1270                            int nRow = 0;
1271                            // Column copy
1272                            const double * element = solver->getMatrixByCol()->getElements();
1273                            const int * row = solver->getMatrixByCol()->getIndices();
1274                            const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
1275                            const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
1276                            int n = 0;
1277                            double contrib = 0.0;
1278                            for (  i = 0 ; i < numberColumns ; i++ ) {
1279                                double value = newSolution[i];
1280                                if (solver->isInteger(i)) {
1281                                    double nearest = floor(value + 0.5);
1282                                    if (fabs(value - nearest) > 1.0e-6) {
1283                                        //printf("Column %d value %g\n",i,value);
1284                                        for (CoinBigIndex j = columnStart[i];
1285                                                j < columnStart[i] + columnLength[i]; j++) {
1286                                            int iRow = row[j];
1287                                            //printf("row %d element %g\n",iRow,element[j]);
1288                                            if (!rowValue[iRow]) {
1289                                                rowValue[iRow] = 1.0;
1290                                                whichRow[nRow++] = iRow;
1291                                            }
1292                                        }
1293                                        baseValue[n] = floor(value);
1294                                        contrib += saveObjective[i] * value;
1295                                        value = 0.0;
1296                                        stack[n] = 0;
1297                                        which[n++] = i;
1298                                    }
1299                                }
1300                                for (CoinBigIndex j = columnStart[i];
1301                                        j < columnStart[i] + columnLength[i]; j++) {
1302                                    int iRow = row[j];
1303                                    rowActivity[iRow] += value * element[j];
1304                                }
1305                            }
1306                            if (newNumberInfeas < 15) {
1307                                stack[n] = newNumberInfeas + 100;
1308                                int iStack = n;
1309                                memset(rowValue, 0, numberRows*sizeof(double));
1310                                const double * rowLower = solver->getRowLower();
1311                                const double * rowUpper = solver->getRowUpper();
1312                                while (iStack >= 0) {
1313                                    double contrib2 = 0.0;
1314                                    // Could do faster
1315                                    for (int k = 0 ; k < n ; k++ ) {
1316                                        i = which[k];
1317                                        double value = baseValue[k] + stack[k];
1318                                        contrib2 += saveObjective[i] * value;
1319                                        for (CoinBigIndex j = columnStart[i];
1320                                                j < columnStart[i] + columnLength[i]; j++) {
1321                                            int iRow = row[j];
1322                                            rowValue[iRow] += value * element[j];
1323                                        }
1324                                    }
1325                                    // check if feasible
1326                                    bool feasible = true;
1327                                    for (int k = 0; k < nRow; k++) {
1328                                        i = whichRow[k];
1329                                        double value = rowValue[i] + rowActivity[i];
1330                                        rowValue[i] = 0.0;
1331                                        if (value < rowLower[i] - 1.0e-7 ||
1332                                                value > rowUpper[i] + 1.0e-7)
1333                                            feasible = false;
1334                                    }
1335                                    if (feasible) {
1336                                        double newObj = newTrueSolutionValue * direction;
1337                                        newObj += contrib2 - contrib;
1338                                        newObj *= direction;
1339#ifdef COIN_DEVELOP
1340                                        printf("FFFeasible! - obj %g\n", newObj);
1341#endif
1342                                        if (newObj < roundingObjective - 1.0e-6) {
1343#ifdef COIN_DEVELOP
1344                                            printf("FBetter\n");
1345#endif
1346                                            roundingObjective = newObj;
1347                                            memcpy(roundingSolution, newSolution, numberColumns*sizeof(double));
1348                                            for (int k = 0 ; k < n ; k++ ) {
1349                                                i = which[k];
1350                                                double value = baseValue[k] + stack[k];
1351                                                roundingSolution[i] = value;
1352                                            }
1353                                        }
1354                                    }
1355                                    while (iStack >= 0 && stack[iStack]) {
1356                                        stack[iStack]--;
1357                                        iStack--;
1358                                    }
1359                                    if (iStack >= 0) {
1360                                        stack[iStack] = 1;
1361                                        iStack = n;
1362                                        stack[n] = 1;
1363                                    }
1364                                }
1365                            }
1366                            delete [] rowActivity;
1367                            delete [] which;
1368                            delete [] stack;
1369                            delete [] baseValue;
1370                            delete [] whichRow;
1371                            delete [] rowValue;
1372                        }
1373                    }
1374                    if (true) {
1375                        OsiSolverInterface * saveSolver = model_->swapSolver(clonedSolver);
1376                        clonedSolver->setColSolution(solver->getColSolution());
1377                        CbcRounding heuristic1(*model_);
1378                        heuristic1.setHeuristicName("rounding in feaspump!");
1379                        heuristic1.setWhen(1);
1380                        roundingObjective = CoinMin(roundingObjective, solutionValue);
1381                        double testSolutionValue = newTrueSolutionValue;
1382                        int returnCode = heuristic1.solution(roundingObjective,
1383                                                             roundingSolution,
1384                                                             testSolutionValue) ;
1385                        if (returnCode == 1) {
1386#ifdef COIN_DEVELOP
1387                            printf("rounding obj of %g?\n", roundingObjective);
1388#endif
1389                            //roundingObjective = newSolutionValue;
1390                            //} else {
1391                            //roundingObjective = COIN_DBL_MAX;
1392                        }
1393                        model_->swapSolver(saveSolver);
1394                    }
1395                    if (!solver->isProvenOptimal()) {
1396                        // presumably max time or some such
1397                        exitAll = true;
1398                        break;
1399                    }
1400                    // in case very dubious solver
1401                    lower = solver->getColLower();
1402                    upper = solver->getColUpper();
1403                    solution = solver->getColSolution();
1404                    numberIterations = solver->getIterationCount();
1405                } else {
1406                    int * addStart = new int[2*general+1];
1407                    int * addIndex = new int[4*general];
1408                    double * addElement = new double[4*general];
1409                    double * addLower = new double[2*general];
1410                    double * addUpper = new double[2*general];
1411                    double * obj = new double[general];
1412                    int nAdd = 0;
1413                    for (i = 0; i < numberIntegers; i++) {
1414                        int iColumn = integerVariable[i];
1415                        if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
1416                                newSolution[iColumn] < upper[iColumn] - primalTolerance) {
1417                            assert (upper[iColumn] - lower[iColumn] > 1.00001);
1418                            obj[nAdd] = 1.0;
1419                            addLower[nAdd] = 0.0;
1420                            addUpper[nAdd] = COIN_DBL_MAX;
1421                            nAdd++;
1422                        }
1423                    }
1424                    OsiSolverInterface * solver2 = solver;
1425                    if (nAdd) {
1426                        CoinZeroN(addStart, nAdd + 1);
1427                        solver2 = solver->clone();
1428                        solver2->addCols(nAdd, addStart, NULL, NULL, addLower, addUpper, obj);
1429                        // feasible solution
1430                        double * sol = new double[nAdd+numberColumns];
1431                        memcpy(sol, solution, numberColumns*sizeof(double));
1432                        // now rows
1433                        int nAdd = 0;
1434                        int nEl = 0;
1435                        int nAddRow = 0;
1436                        for (i = 0; i < numberIntegers; i++) {
1437                            int iColumn = integerVariable[i];
1438                            if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
1439                                    newSolution[iColumn] < upper[iColumn] - primalTolerance) {
1440                                addLower[nAddRow] = -newSolution[iColumn];;
1441                                addUpper[nAddRow] = COIN_DBL_MAX;
1442                                addIndex[nEl] = iColumn;
1443                                addElement[nEl++] = -1.0;
1444                                addIndex[nEl] = numberColumns + nAdd;
1445                                addElement[nEl++] = 1.0;
1446                                nAddRow++;
1447                                addStart[nAddRow] = nEl;
1448                                addLower[nAddRow] = newSolution[iColumn];;
1449                                addUpper[nAddRow] = COIN_DBL_MAX;
1450                                addIndex[nEl] = iColumn;
1451                                addElement[nEl++] = 1.0;
1452                                addIndex[nEl] = numberColumns + nAdd;
1453                                addElement[nEl++] = 1.0;
1454                                nAddRow++;
1455                                addStart[nAddRow] = nEl;
1456                                sol[nAdd+numberColumns] = fabs(sol[iColumn] - newSolution[iColumn]);
1457                                nAdd++;
1458                            }
1459                        }
1460                        solver2->setColSolution(sol);
1461                        delete [] sol;
1462                        solver2->addRows(nAddRow, addStart, addIndex, addElement, addLower, addUpper);
1463                    }
1464                    delete [] addStart;
1465                    delete [] addIndex;
1466                    delete [] addElement;
1467                    delete [] addLower;
1468                    delete [] addUpper;
1469                    delete [] obj;
1470                    solver2->resolve();
1471                    if (!solver2->isProvenOptimal()) {
1472                        // presumably max time or some such
1473                        exitAll = true;
1474                        break;
1475                    }
1476                    //assert (solver2->isProvenOptimal());
1477                    if (nAdd) {
1478                        solver->setColSolution(solver2->getColSolution());
1479                        numberIterations = solver2->getIterationCount();
1480                        delete solver2;
1481                    } else {
1482                        numberIterations = solver->getIterationCount();
1483                    }
1484                    lower = solver->getColLower();
1485                    upper = solver->getColUpper();
1486                    solution = solver->getColSolution();
1487                    newTrueSolutionValue = -saveOffset;
1488                    newSumInfeas = 0.0;
1489                    newNumberInfeas = 0;
1490                    {
1491                        const double * newSolution = solver->getColSolution();
1492                        for (  i = 0 ; i < numberColumns ; i++ ) {
1493                            if (solver->isInteger(i)) {
1494                                double value = newSolution[i];
1495                                double nearest = floor(value + 0.5);
1496                                newSumInfeas += fabs(value - nearest);
1497                                if (fabs(value - nearest) > 1.0e-6)
1498                                    newNumberInfeas++;
1499                            }
1500                            newTrueSolutionValue += saveObjective[i] * newSolution[i];
1501                        }
1502                        newTrueSolutionValue *= direction;
1503                    }
1504                }
1505                if (lastMove != 1000000) {
1506                    if (newSumInfeas < lastSumInfeas) {
1507                        lastMove = numberPasses;
1508                        lastSumInfeas = newSumInfeas;
1509                    } else if (newSumInfeas > lastSumInfeas + 1.0e-5) {
1510                        lastMove = 1000000; // going up
1511                    }
1512                }
1513                totalNumberIterations += numberIterations;
1514                if (solver->getNumRows() < 3000)
1515                    sprintf(pumpPrint, "Pass %3d: suminf. %10.5f (%d) obj. %g iterations %d",
1516                            numberPasses + totalNumberPasses,
1517                            newSumInfeas, newNumberInfeas,
1518                            newTrueSolutionValue, numberIterations);
1519                else
1520                    sprintf(pumpPrint, "Pass %3d: (%.2f seconds) suminf. %10.5f (%d) obj. %g iterations %d", numberPasses + totalNumberPasses,
1521                            model_->getCurrentSeconds(), newSumInfeas, newNumberInfeas,
1522                            newTrueSolutionValue, numberIterations);
1523                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1524                << pumpPrint
1525                << CoinMessageEol;
1526                if (closestSolution && solver->getObjValue() < closestObjectiveValue) {
1527                    int i;
1528                    const double * objective = solver->getObjCoefficients();
1529                    for (i = 0; i < numberIntegersOrig; i++) {
1530                        int iColumn = integerVariableOrig[i];
1531                        if (objective[iColumn] > 0.0)
1532                            closestSolution[i] = 0;
1533                        else
1534                            closestSolution[i] = 1;
1535                    }
1536                    closestObjectiveValue = solver->getObjValue();
1537                }
1538                // See if we need to think about changing rhs
1539                if ((switches_&12) != 0 && useRhs < 1.0e50) {
1540                    double oldRhs = useRhs;
1541                    bool trying = false;
1542                    if ((switches_&4) != 0 && numberPasses && (numberPasses % 50) == 0) {
1543                        if (solutionValue > 1.0e20) {
1544                            // only if no genuine solution
1545                            double gap = useRhs - continuousObjectiveValue;
1546                            useRhs += 0.1 * gap;
1547                            if (exactMultiple) {
1548                                useRhs = exactMultiple * ceil(useRhs / exactMultiple);
1549                                useRhs = CoinMax(useRhs, oldRhs + exactMultiple);
1550                            }
1551                            trying = true;
1552                        }
1553                    }
1554                    if ((switches_&8) != 0) {
1555                        // Put in new suminf and check
1556                        double largest = newSumInfeas;
1557                        double smallest = newSumInfeas;
1558                        for (int i = 0; i < SIZE_BOBBLE - 1; i++) {
1559                            double value = saveSumInf[i+1];
1560                            saveSumInf[i] = value;
1561                            largest = CoinMax(largest, value);
1562                            smallest = CoinMin(smallest, value);
1563                        }
1564                        saveSumInf[SIZE_BOBBLE-1] = newSumInfeas;
1565                        if (smallest*1.5 > largest && smallest > 2.0) {
1566                            if (bobbleMode == 0) {
1567                                // go closer
1568                                double gap = oldRhs - continuousObjectiveValue;
1569                                useRhs -= 0.4 * gap;
1570                                if (exactMultiple) {
1571                                    double value = floor(useRhs / exactMultiple);
1572                                    useRhs = CoinMin(value * exactMultiple, oldRhs - exactMultiple);
1573                                }
1574                                if (useRhs < continuousObjectiveValue) {
1575                                    // skip decrease
1576                                    bobbleMode = 1;
1577                                    useRhs = oldRhs;
1578                                }
1579                            }
1580                            if (bobbleMode) {
1581                                trying = true;
1582                                // weaken
1583                                if (solutionValue < 1.0e20) {
1584                                    double gap = solutionValue - oldRhs;
1585                                    useRhs += 0.3 * gap;
1586                                } else {
1587                                    double gap = oldRhs - continuousObjectiveValue;
1588                                    useRhs += 0.05 * gap;
1589                                }
1590                                if (exactMultiple) {
1591                                    double value = ceil(useRhs / exactMultiple);
1592                                    useRhs = CoinMin(value * exactMultiple,
1593                                                     solutionValue - exactMultiple);
1594                                }
1595                            }
1596                            bobbleMode++;
1597                            // reset
1598                            CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
1599                        }
1600                    }
1601                    if (useRhs != oldRhs) {
1602                        // tidy up
1603                        if (exactMultiple) {
1604                            double value = floor(useRhs / exactMultiple);
1605                            double bestPossible = ceil(continuousObjectiveValue / exactMultiple);
1606                            useRhs = CoinMax(value, bestPossible) * exactMultiple;
1607                        } else {
1608                            useRhs = CoinMax(useRhs, continuousObjectiveValue);
1609                        }
1610                        int k = solver->getNumRows() - 1;
1611                        solver->setRowUpper(k, useRhs + useOffset);
1612                        bool takeHint;
1613                        OsiHintStrength strength;
1614                        solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
1615                        if (useRhs < oldRhs) {
1616                            solver->setHintParam(OsiDoDualInResolve, true);
1617                            solver->resolve();
1618                        } else if (useRhs > oldRhs) {
1619                            solver->setHintParam(OsiDoDualInResolve, false);
1620                            solver->resolve();
1621                        }
1622                        solver->setHintParam(OsiDoDualInResolve, takeHint);
1623                        if (!solver->isProvenOptimal()) {
1624                            // presumably max time or some such
1625                            exitAll = true;
1626                            break;
1627                        }
1628                    } else if (trying) {
1629                        // doesn't look good
1630                        break;
1631                    }
1632                }
1633            }
1634            // reduce scale factor
1635            scaleFactor *= weightFactor_;
1636        } // END WHILE
1637        // see if rounding worked!
1638        if (roundingObjective < solutionValue) {
1639            if (roundingObjective < solutionValue - 1.0e-6*fabs(roundingObjective)) {
1640                sprintf(pumpPrint, "Rounding solution of %g is better than previous of %g\n",
1641                        roundingObjective, solutionValue);
1642                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1643                << pumpPrint
1644                << CoinMessageEol;
1645            }
1646            solutionValue = roundingObjective;
1647            newSolutionValue = solutionValue;
1648            memcpy(betterSolution, roundingSolution, numberColumns*sizeof(double));
1649            solutionFound = true;
1650            numberSolutions++;
1651            if (numberSolutions>=maxSolutions)
1652              exitAll = true;
1653            if (exitNow(roundingObjective))
1654                exitAll = true;
1655        }
1656        if (!solutionFound) {
1657            sprintf(pumpPrint, "No solution found this major pass");
1658            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1659            << pumpPrint
1660            << CoinMessageEol;
1661        }
1662        //}
1663        delete solver;
1664        solver = NULL;
1665        for ( j = 0; j < NUMBER_OLD; j++)
1666            delete [] oldSolution[j];
1667        delete [] oldSolution;
1668        delete [] saveObjective;
1669        if (usedColumn && !exitAll) {
1670            OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1671            const double * colLower = newSolver->getColLower();
1672            const double * colUpper = newSolver->getColUpper();
1673            bool stopBAB = false;
1674            int allowedPass = -1;
1675            if (maximumAllowed > 0)
1676                allowedPass = CoinMax(numberPasses - maximumAllowed, -1);
1677            while (!stopBAB) {
1678                stopBAB = true;
1679                int i;
1680                int nFix = 0;
1681                int nFixI = 0;
1682                int nFixC = 0;
1683                int nFixC2 = 0;
1684                for (i = 0; i < numberIntegersOrig; i++) {
1685                    int iColumn = integerVariableOrig[i];
1686                    //const OsiObject * object = model_->object(i);
1687                    //double originalLower;
1688                    //double originalUpper;
1689                    //getIntegerInformation( object,originalLower, originalUpper);
1690                    //assert(colLower[iColumn]==originalLower);
1691                    //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
1692                    newSolver->setColLower(iColumn, colLower[iColumn]);
1693                    //assert(colUpper[iColumn]==originalUpper);
1694                    //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
1695                    newSolver->setColUpper(iColumn, colUpper[iColumn]);
1696                    if (usedColumn[iColumn] <= allowedPass) {
1697                        double value = lastSolution[iColumn];
1698                        double nearest = floor(value + 0.5);
1699                        if (fabs(value - nearest) < 1.0e-7) {
1700                            if (nearest == colLower[iColumn]) {
1701                                newSolver->setColUpper(iColumn, colLower[iColumn]);
1702                                nFix++;
1703                            } else if (nearest == colUpper[iColumn]) {
1704                                newSolver->setColLower(iColumn, colUpper[iColumn]);
1705                                nFix++;
1706                            } else if (fixInternal) {
1707                                newSolver->setColLower(iColumn, nearest);
1708                                newSolver->setColUpper(iColumn, nearest);
1709                                nFix++;
1710                                nFixI++;
1711                            }
1712                        }
1713                    }
1714                }
1715                if (fixContinuous) {
1716                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1717                        if (!newSolver->isInteger(iColumn) && usedColumn[iColumn] <= allowedPass) {
1718                            double value = lastSolution[iColumn];
1719                            if (value < colLower[iColumn] + 1.0e-8) {
1720                                newSolver->setColUpper(iColumn, colLower[iColumn]);
1721                                nFixC++;
1722                            } else if (value > colUpper[iColumn] - 1.0e-8) {
1723                                newSolver->setColLower(iColumn, colUpper[iColumn]);
1724                                nFixC++;
1725                            } else if (fixContinuous == 2) {
1726                                newSolver->setColLower(iColumn, value);
1727                                newSolver->setColUpper(iColumn, value);
1728                                nFixC++;
1729                                nFixC2++;
1730                            }
1731                        }
1732                    }
1733                }
1734                newSolver->initialSolve();
1735                if (!newSolver->isProvenOptimal()) {
1736                    //newSolver->writeMps("bad.mps");
1737                    //assert (newSolver->isProvenOptimal());
1738                    exitAll = true;
1739                    break;
1740                }
1741                sprintf(pumpPrint, "Before mini branch and bound, %d integers at bound fixed and %d continuous",
1742                        nFix, nFixC);
1743                if (nFixC2 + nFixI != 0)
1744                    sprintf(pumpPrint + strlen(pumpPrint), " of which %d were internal integer and %d internal continuous",
1745                            nFixI, nFixC2);
1746                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1747                << pumpPrint
1748                << CoinMessageEol;
1749                double saveValue = newSolutionValue;
1750                if (newSolutionValue - model_->getCutoffIncrement()
1751                        > continuousObjectiveValue - 1.0e-7) {
1752                    double saveFraction = fractionSmall_;
1753                    if (numberTries > 1 && !numberBandBsolutions)
1754                        fractionSmall_ *= 0.5;
1755                    // Give branch and bound a bit more freedom
1756                    double cutoff2 = newSolutionValue +
1757                                     CoinMax(model_->getCutoffIncrement(), 1.0e-3);
1758                    int returnCode2 = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
1759                                                          cutoff2, "CbcHeuristicLocalAfterFPump");
1760                    fractionSmall_ = saveFraction;
1761                    if (returnCode2 < 0) {
1762                        if (returnCode2 == -2) {
1763                            exitAll = true;
1764                            returnCode = 0;
1765                        } else {
1766                            returnCode2 = 0; // returned on size - try changing
1767                            //#define ROUND_AGAIN
1768#ifdef ROUND_AGAIN
1769                            if (numberTries == 1 && numberPasses > 20 && allowedPass < numberPasses - 1) {
1770                                allowedPass = (numberPasses + allowedPass) >> 1;
1771                                sprintf(pumpPrint,
1772                                        "Fixing all variables which were last changed on pass %d and trying again",
1773                                        allowedPass);
1774                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1775                                << pumpPrint
1776                                << CoinMessageEol;
1777                                stopBAB = false;
1778                                continue;
1779                            }
1780#endif
1781                        }
1782                    }
1783                    if ((returnCode2&2) != 0) {
1784                        // could add cut
1785                        returnCode2 &= ~2;
1786                    }
1787                    if (returnCode2)
1788                        numberBandBsolutions++;
1789                } else {
1790                    // no need
1791                    exitAll = true;
1792                    //returnCode=0;
1793                }
1794                // recompute solution value
1795                if (returnCode && true) {
1796                    delete newSolver;
1797                    newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1798                    newSolutionValue = -saveOffset;
1799                    double newSumInfeas = 0.0;
1800                    const double * obj = newSolver->getObjCoefficients();
1801                    for (int i = 0 ; i < numberColumns ; i++ ) {
1802                        if (newSolver->isInteger(i)) {
1803                            double value = newSolution[i];
1804                            double nearest = floor(value + 0.5);
1805                            newSumInfeas += fabs(value - nearest);
1806                        }
1807                        newSolutionValue += obj[i] * newSolution[i];
1808                    }
1809                    newSolutionValue *= direction;
1810                }
1811                bool gotSolution = false;
1812                if (returnCode && newSolutionValue < saveValue) {
1813                    sprintf(pumpPrint, "Mini branch and bound improved solution from %g to %g (%.2f seconds)",
1814                            saveValue, newSolutionValue, model_->getCurrentSeconds());
1815                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1816                    << pumpPrint
1817                    << CoinMessageEol;
1818                    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1819                    gotSolution = true;
1820                    if (fixContinuous && nFixC + nFixC2 > 0) {
1821                        // may be able to do even better
1822                        int nFixed = 0;
1823                        const double * lower = model_->solver()->getColLower();
1824                        const double * upper = model_->solver()->getColUpper();
1825                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1826                            double value = newSolution[iColumn];
1827                            if (newSolver->isInteger(iColumn)) {
1828                                value = floor(newSolution[iColumn] + 0.5);
1829                                newSolver->setColLower(iColumn, value);
1830                                newSolver->setColUpper(iColumn, value);
1831                                nFixed++;
1832                            } else {
1833                                newSolver->setColLower(iColumn, lower[iColumn]);
1834                                newSolver->setColUpper(iColumn, upper[iColumn]);
1835                                if (value < lower[iColumn])
1836                                    value = lower[iColumn];
1837                                else if (value > upper[iColumn])
1838                                    value = upper[iColumn];
1839
1840                            }
1841                            newSolution[iColumn] = value;
1842                        }
1843                        newSolver->setColSolution(newSolution);
1844                        //#define CLP_INVESTIGATE2
1845#ifdef CLP_INVESTIGATE2
1846                        {
1847                            // check
1848                            // get row activities
1849                            int numberRows = newSolver->getNumRows();
1850                            double * rowActivity = new double[numberRows];
1851                            memset(rowActivity, 0, numberRows*sizeof(double));
1852                            newSolver->getMatrixByCol()->times(newSolution, rowActivity) ;
1853                            double largestInfeasibility = primalTolerance;
1854                            double sumInfeasibility = 0.0;
1855                            int numberBadRows = 0;
1856                            const double * rowLower = newSolver->getRowLower();
1857                            const double * rowUpper = newSolver->getRowUpper();
1858                            for (i = 0 ; i < numberRows ; i++) {
1859                                double value;
1860                                value = rowLower[i] - rowActivity[i];
1861                                if (value > primalTolerance) {
1862                                    numberBadRows++;
1863                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
1864                                    sumInfeasibility += value;
1865                                }
1866                                value = rowActivity[i] - rowUpper[i];
1867                                if (value > primalTolerance) {
1868                                    numberBadRows++;
1869                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
1870                                    sumInfeasibility += value;
1871                                }
1872                            }
1873                            printf("%d bad rows, largest inf %g sum %g\n",
1874                                   numberBadRows, largestInfeasibility, sumInfeasibility);
1875                            delete [] rowActivity;
1876                        }
1877#endif
1878#ifdef COIN_HAS_CLP
1879                        OsiClpSolverInterface * clpSolver
1880                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1881                        if (clpSolver) {
1882                            ClpSimplex * simplex = clpSolver->getModelPtr();
1883                            //simplex->writeBasis("test.bas",true,2);
1884                            //simplex->writeMps("test.mps",2,1);
1885                            if (nFixed*3 > numberColumns*2)
1886                                simplex->allSlackBasis(); // may as well go from all slack
1887                            simplex->primal(1);
1888                            clpSolver->setWarmStart(NULL);
1889                        }
1890#endif
1891                        newSolver->initialSolve();
1892                        if (newSolver->isProvenOptimal()) {
1893                            double value = newSolver->getObjValue() * newSolver->getObjSense();
1894                            if (value < newSolutionValue) {
1895                                //newSolver->writeMps("query","mps");
1896#ifdef JJF_ZERO
1897                                {
1898                                    double saveOffset;
1899                                    newSolver->getDblParam(OsiObjOffset, saveOffset);
1900                                    const double * obj = newSolver->getObjCoefficients();
1901                                    double newTrueSolutionValue = -saveOffset;
1902                                    double newSumInfeas = 0.0;
1903                                    int numberColumns = newSolver->getNumCols();
1904                                    const double * solution = newSolver->getColSolution();
1905                                    for (int  i = 0 ; i < numberColumns ; i++ ) {
1906                                        if (newSolver->isInteger(i)) {
1907                                            double value = solution[i];
1908                                            double nearest = floor(value + 0.5);
1909                                            newSumInfeas += fabs(value - nearest);
1910                                        }
1911                                        if (solution[i])
1912                                            printf("%d obj %g val %g - total %g\n", i, obj[i], solution[i],
1913                                                   newTrueSolutionValue);
1914                                        newTrueSolutionValue += obj[i] * solution[i];
1915                                    }
1916                                    printf("obj %g - inf %g\n", newTrueSolutionValue,
1917                                           newSumInfeas);
1918                                }
1919#endif
1920                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
1921                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1922                                << pumpPrint
1923                                << CoinMessageEol;
1924                                newSolutionValue = value;
1925                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
1926                            }
1927                        } else {
1928                            //newSolver->writeMps("bad3.mps");
1929                            exitAll = true;
1930                            break;
1931                        }
1932                    }
1933                } else {
1934                    sprintf(pumpPrint, "Mini branch and bound did not improve solution (%.2f seconds)",
1935                            model_->getCurrentSeconds());
1936                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1937                    << pumpPrint
1938                    << CoinMessageEol;
1939                    if (returnCode && newSolutionValue < saveValue + 1.0e-3 && nFixC + nFixC2) {
1940                        // may be able to do better
1941                        const double * lower = model_->solver()->getColLower();
1942                        const double * upper = model_->solver()->getColUpper();
1943                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1944                            if (newSolver->isInteger(iColumn)) {
1945                                double value = floor(newSolution[iColumn] + 0.5);
1946                                newSolver->setColLower(iColumn, value);
1947                                newSolver->setColUpper(iColumn, value);
1948                            } else {
1949                                newSolver->setColLower(iColumn, lower[iColumn]);
1950                                newSolver->setColUpper(iColumn, upper[iColumn]);
1951                            }
1952                        }
1953                        newSolver->initialSolve();
1954                        if (newSolver->isProvenOptimal()) {
1955                            double value = newSolver->getObjValue() * newSolver->getObjSense();
1956                            if (value < saveValue) {
1957                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
1958                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1959                                << pumpPrint
1960                                << CoinMessageEol;
1961                                newSolutionValue = value;
1962                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
1963                                gotSolution = true;
1964                            }
1965                        }
1966                    }
1967                }
1968                if (gotSolution) {
1969                    if ((accumulate_&1) != 0) {
1970                        model_->incrementUsed(betterSolution); // for local search
1971                    }
1972                    solutionValue = newSolutionValue;
1973                    solutionFound = true;
1974                    numberSolutions++;
1975                    if (numberSolutions>=maxSolutions)
1976                      exitAll = true;
1977                    if (exitNow(newSolutionValue))
1978                        exitAll = true;
1979                    CoinWarmStartBasis * basis =
1980                        dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
1981                    if (basis) {
1982                        bestBasis = * basis;
1983                        delete basis;
1984                        int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
1985                        if (action == 0) {
1986                            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
1987                            double saveObjectiveValue = model_->getMinimizationObjValue();
1988                            model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
1989                            if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
1990                                model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
1991                            delete [] saveOldSolution;
1992                        }
1993                        if (!action || model_->getCurrentSeconds() > model_->getMaximumSeconds()) {
1994                            exitAll = true; // exit
1995                            break;
1996                        }
1997                    }
1998                }
1999            } // end stopBAB while
2000            delete newSolver;
2001        }
2002        if (solutionFound) finalReturnCode = 1;
2003        cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
2004        if (numberTries >= maximumRetries_ || !solutionFound || exitAll || cutoff < continuousObjectiveValue + 1.0e-7) {
2005            break;
2006        } else {
2007            solutionFound = false;
2008            if (absoluteIncrement_ > 0.0 || relativeIncrement_ > 0.0) {
2009                double gap = relativeIncrement_ * fabs(solutionValue);
2010                double change = CoinMax(gap, absoluteIncrement_);
2011                cutoff = CoinMin(cutoff, solutionValue - change);
2012            } else {
2013                //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
2014                double weights[10] = {0.1, 0.2, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5, 0.6};
2015                cutoff -= weights[CoinMin(numberTries-1, 9)] * (cutoff - continuousObjectiveValue);
2016            }
2017            // But round down
2018            if (exactMultiple)
2019                cutoff = exactMultiple * floor(cutoff / exactMultiple);
2020            if (cutoff < continuousObjectiveValue)
2021                break;
2022            sprintf(pumpPrint, "Round again with cutoff of %g", cutoff);
2023            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2024            << pumpPrint
2025            << CoinMessageEol;
2026            if ((accumulate_&3) < 2 && usedColumn) {
2027                for (int i = 0; i < numberColumns; i++)
2028                    usedColumn[i] = -1;
2029            }
2030            averageIterationsPerTry = totalNumberIterations / numberTries;
2031            numberIterationsLastPass =  totalNumberIterations;
2032            totalNumberPasses += numberPasses - 1;
2033        }
2034    }
2035/*
2036  End of the `exitAll' loop.
2037*/
2038#ifdef RAND_RAND
2039    delete [] randomFactor;
2040#endif
2041    delete solver; // probably NULL but do anyway
2042    if (!finalReturnCode && closestSolution && closestObjectiveValue <= 10.0 &&
2043            usedColumn && !model_->maximumSecondsReached()) {
2044        // try a bit of branch and bound
2045        OsiSolverInterface * newSolver = cloneBut(1); // was model_->continuousSolver()->clone();
2046        const double * colLower = newSolver->getColLower();
2047        const double * colUpper = newSolver->getColUpper();
2048        int i;
2049        double rhs = 0.0;
2050        for (i = 0; i < numberIntegersOrig; i++) {
2051            int iColumn = integerVariableOrig[i];
2052            int direction = closestSolution[i];
2053            closestSolution[i] = iColumn;
2054            if (direction == 0) {
2055                // keep close to LB
2056                rhs += colLower[iColumn];
2057                lastSolution[i] = 1.0;
2058            } else {
2059                // keep close to UB
2060                rhs -= colUpper[iColumn];
2061                lastSolution[i] = -1.0;
2062            }
2063        }
2064        newSolver->addRow(numberIntegersOrig, closestSolution,
2065                          lastSolution, -COIN_DBL_MAX, rhs + 10.0);
2066        //double saveValue = newSolutionValue;
2067        //newSolver->writeMps("sub");
2068        int returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
2069                                             newSolutionValue, "CbcHeuristicLocalAfterFPump");
2070        if (returnCode < 0)
2071            returnCode = 0; // returned on size
2072        if ((returnCode&2) != 0) {
2073            // could add cut
2074            returnCode &= ~2;
2075        }
2076        if (returnCode) {
2077            //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
2078            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2079            //abort();
2080            solutionValue = newSolutionValue;
2081            solutionFound = true;
2082            numberSolutions++;
2083            if (numberSolutions>=maxSolutions)
2084              exitAll = true;
2085            if (exitNow(newSolutionValue))
2086                exitAll = true;
2087        }
2088        delete newSolver;
2089    }
2090    delete clonedSolver;
2091    delete [] roundingSolution;
2092    delete [] usedColumn;
2093    delete [] lastSolution;
2094    delete [] newSolution;
2095    delete [] closestSolution;
2096    delete [] integerVariable;
2097    delete [] firstPerturbedObjective;
2098    delete [] firstPerturbedSolution;
2099    if (solutionValue == incomingObjective)
2100        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
2101                model_->getCurrentSeconds(), CoinCpuTime() - time1);
2102    else if (numberSolutions < maxSolutions)
2103      sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
2104              model_->getCurrentSeconds(), solutionValue, CoinCpuTime() - time1);
2105    else 
2106        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g (stopping after %d solutions) - took %.2f seconds",
2107                model_->getCurrentSeconds(), solutionValue, 
2108                numberSolutions,CoinCpuTime() - time1);
2109    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2110      << pumpPrint
2111      << CoinMessageEol;
2112    if (bestBasis.getNumStructural())
2113        model_->setBestSolutionBasis(bestBasis);
2114    //model_->setMinimizationObjValue(saveBestObjective);
2115    if ((accumulate_&1) != 0 && numberSolutions > 1 && !model_->getSolutionCount()) {
2116        model_->setSolutionCount(1); // for local search
2117        model_->setNumberHeuristicSolutions(1);
2118    }
2119#ifdef COIN_DEVELOP
2120    {
2121        double ncol = model_->solver()->getNumCols();
2122        double nrow = model_->solver()->getNumRows();
2123        printf("XXX total iterations %d ratios - %g %g %g\n",
2124               totalNumberIterations,
2125               static_cast<double> (totalNumberIterations) / nrow,
2126               static_cast<double> (totalNumberIterations) / ncol,
2127               static_cast<double> (totalNumberIterations) / (2*nrow + 2*ncol));
2128    }
2129#endif
2130    return finalReturnCode;
2131}
2132
2133/**************************END MAIN PROCEDURE ***********************************/
2134
2135// update model
2136void CbcHeuristicFPump::setModel(CbcModel * model)
2137{
2138    model_ = model;
2139}
2140
2141/* Rounds solution - down if < downValue
2142   returns 1 if current is a feasible solution
2143*/
2144int
2145CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution,
2146                          //const double * objective,
2147                          int numberIntegers, const int * integerVariable,
2148                          /*char * pumpPrint,*/ int iter,
2149                          /*bool roundExpensive,*/ double downValue, int *flip)
2150{
2151    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2152    double primalTolerance ;
2153    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
2154
2155    int i;
2156
2157    const double * cost = solver->getObjCoefficients();
2158    int flip_up = 0;
2159    int flip_down  = 0;
2160    double  v = randomNumberGenerator_.randomDouble() * 20.0;
2161    int nn = 10 + static_cast<int> (v);
2162    int nnv = 0;
2163    int * list = new int [nn];
2164    double * val = new double [nn];
2165    for (i = 0; i < nn; i++) val[i] = .001;
2166
2167    const double * rowLower = solver->getRowLower();
2168    const double * rowUpper = solver->getRowUpper();
2169    int numberRows = solver->getNumRows();
2170    if (false && (iter&1) != 0) {
2171        // Do set covering variables
2172        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2173        const double * elementByRow = matrixByRow->getElements();
2174        const int * column = matrixByRow->getIndices();
2175        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2176        const int * rowLength = matrixByRow->getVectorLengths();
2177        for (i = 0; i < numberRows; i++) {
2178            if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2179                bool cover = true;
2180                double largest = 0.0;
2181                int jColumn = -1;
2182                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2183                    int iColumn = column[k];
2184                    if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
2185                        cover = false;
2186                        break;
2187                    } else {
2188                        if (solution[iColumn]) {
2189                            double value = solution[iColumn] *
2190                                           (randomNumberGenerator_.randomDouble() + 5.0);
2191                            if (value > largest) {
2192                                largest = value;
2193                                jColumn = iColumn;
2194                            }
2195                        }
2196                    }
2197                }
2198                if (cover) {
2199                    for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2200                        int iColumn = column[k];
2201                        if (iColumn == jColumn)
2202                            solution[iColumn] = 1.0;
2203                        else
2204                            solution[iColumn] = 0.0;
2205                    }
2206                }
2207            }
2208        }
2209    }
2210    int numberColumns = solver->getNumCols();
2211#ifdef JJF_ZERO
2212    // Do set covering variables
2213    const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2214    const double * elementByRow = matrixByRow->getElements();
2215    const int * column = matrixByRow->getIndices();
2216    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2217    const int * rowLength = matrixByRow->getVectorLengths();
2218    double * sortTemp = new double[numberColumns];
2219    int * whichTemp = new int [numberColumns];
2220    char * rowTemp = new char [numberRows];
2221    memset(rowTemp, 0, numberRows);
2222    for (i = 0; i < numberColumns; i++)
2223        whichTemp[i] = -1;
2224    int nSOS = 0;
2225    for (i = 0; i < numberRows; i++) {
2226        if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2227            bool cover = true;
2228            for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2229                int iColumn = column[k];
2230                if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
2231                    cover = false;
2232                    break;
2233                }
2234            }
2235            if (cover) {
2236                rowTemp[i] = 1;
2237                nSOS++;
2238                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2239                    int iColumn = column[k];
2240                    double value = solution[iColumn];
2241                    whichTemp[iColumn] = iColumn;
2242                }
2243            }
2244        }
2245    }
2246    if (nSOS) {
2247        // Column copy
2248        const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol();
2249        //const double * element = matrixByColumn->getElements();
2250        const int * row = matrixByColumn->getIndices();
2251        const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
2252        const int * columnLength = matrixByColumn->getVectorLengths();
2253        int nLook = 0;
2254        for (i = 0; i < numberColumns; i++) {
2255            if (whichTemp[i] >= 0) {
2256                whichTemp[nLook] = i;
2257                double value = solution[i];
2258                if (value < 0.5)
2259                    value *= (0.1 * randomNumberGenerator_.randomDouble() + 0.3);
2260                sortTemp[nLook++] = -value;
2261            }
2262        }
2263        CoinSort_2(sortTemp, sortTemp + nLook, whichTemp);
2264        double smallest = 1.0;
2265        int nFix = 0;
2266        int nOne = 0;
2267        for (int j = 0; j < nLook; j++) {
2268            int jColumn = whichTemp[j];
2269            double thisValue = solution[jColumn];
2270            if (!thisValue)
2271                continue;
2272            if (thisValue == 1.0)
2273                nOne++;
2274            smallest = CoinMin(smallest, thisValue);
2275            solution[jColumn] = 1.0;
2276            double largest = 0.0;
2277            for (CoinBigIndex jEl = columnStart[jColumn];
2278                    jEl < columnStart[jColumn] + columnLength[jColumn]; jEl++) {
2279                int jRow = row[jEl];
2280                if (rowTemp[jRow]) {
2281                    for (CoinBigIndex k = rowStart[jRow]; k < rowStart[jRow] + rowLength[jRow]; k++) {
2282                        int iColumn = column[k];
2283                        if (solution[iColumn]) {
2284                            if (iColumn != jColumn) {
2285                                double value = solution[iColumn];
2286                                if (value > largest)
2287                                    largest = value;
2288                                solution[iColumn] = 0.0;
2289                            }
2290                        }
2291                    }
2292                }
2293            }
2294            if (largest > thisValue)
2295                printf("%d was at %g - chosen over a value of %g\n",
2296                       jColumn, thisValue, largest);
2297            nFix++;
2298        }
2299        printf("%d fixed out of %d (%d at one already)\n",
2300               nFix, nLook, nOne);
2301    }
2302    delete [] sortTemp;
2303    delete [] whichTemp;
2304    delete [] rowTemp;
2305#endif
2306    const double * columnLower = solver->getColLower();
2307    const double * columnUpper = solver->getColUpper();
2308    // Check if valid with current solution (allow for 0.99999999s)
2309    for (i = 0; i < numberIntegers; i++) {
2310        int iColumn = integerVariable[i];
2311        double value = solution[iColumn];
2312        double round = floor(value + 0.5);
2313        if (fabs(value - round) > primalTolerance)
2314            break;
2315    }
2316    if (i == numberIntegers) {
2317        // may be able to use solution even if 0.99999's
2318        double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
2319        double * saveUpper = CoinCopyOfArray(columnUpper, numberColumns);
2320        double * saveSolution = CoinCopyOfArray(solution, numberColumns);
2321        double * tempSolution = CoinCopyOfArray(solution, numberColumns);
2322        CoinWarmStartBasis  * saveBasis =
2323            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
2324        for (i = 0; i < numberIntegers; i++) {
2325            int iColumn = integerVariable[i];
2326            double value = solution[iColumn];
2327            double round = floor(value + 0.5);
2328            solver->setColLower(iColumn, round);
2329            solver->setColUpper(iColumn, round);
2330            tempSolution[iColumn] = round;
2331        }
2332        solver->setColSolution(tempSolution);
2333        delete [] tempSolution;
2334        solver->resolve();
2335        solver->setColLower(saveLower);
2336        solver->setColUpper(saveUpper);
2337        solver->setWarmStart(saveBasis);
2338        delete [] saveLower;
2339        delete [] saveUpper;
2340        delete saveBasis;
2341        if (!solver->isProvenOptimal()) {
2342            solver->setColSolution(saveSolution);
2343        }
2344        delete [] saveSolution;
2345        if (solver->isProvenOptimal()) {
2346            // feasible
2347            delete [] list;
2348            delete [] val;
2349            return 1;
2350        }
2351    }
2352    //double * saveSolution = CoinCopyOfArray(solution,numberColumns);
2353    // return rounded solution
2354    for (i = 0; i < numberIntegers; i++) {
2355        int iColumn = integerVariable[i];
2356        double value = solution[iColumn];
2357        double round = floor(value + primalTolerance);
2358        if (value - round > downValue) round += 1.;
2359#ifndef JJF_ONE
2360        if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2361        if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2362#else
2363        if (round < columnLower[iColumn] + integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2364        if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2365#endif
2366        if (flip_up + flip_down == 0) {
2367            for (int k = 0; k < nn; k++) {
2368                if (fabs(value - round) > val[k]) {
2369                    nnv++;
2370                    for (int j = nn - 2; j >= k; j--) {
2371                        val[j+1] = val[j];
2372                        list[j+1] = list[j];
2373                    }
2374                    val[k] = fabs(value - round);
2375                    list[k] = iColumn;
2376                    break;
2377                }
2378            }
2379        }
2380        solution[iColumn] = round;
2381    }
2382
2383    if (nnv > nn) nnv = nn;
2384    //if (iter != 0)
2385    //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
2386    *flip = flip_up + flip_down;
2387
2388    if (*flip == 0 && iter != 0) {
2389        //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
2390        for (i = 0; i < nnv; i++) {
2391            // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
2392            int index = list[i];
2393            double value = solution[index];
2394            if (value <= 1.0) {
2395                solution[index] = 1.0 - value;
2396            } else if (value < columnLower[index] + integerTolerance) {
2397                solution[index] = value + 1.0;
2398            } else if (value > columnUpper[index] - integerTolerance) {
2399                solution[index] = value - 1.0;
2400            } else {
2401                solution[index] = value - 1.0;
2402            }
2403        }
2404        *flip = nnv;
2405    } else {
2406        //sprintf(pumpPrint+strlen(pumpPrint)," ");
2407    }
2408    delete [] list;
2409    delete [] val;
2410    //iter++;
2411
2412    // get row activities
2413    double * rowActivity = new double[numberRows];
2414    memset(rowActivity, 0, numberRows*sizeof(double));
2415    solver->getMatrixByCol()->times(solution, rowActivity) ;
2416    double largestInfeasibility = primalTolerance;
2417    double sumInfeasibility = 0.0;
2418    int numberBadRows = 0;
2419    for (i = 0 ; i < numberRows ; i++) {
2420        double value;
2421        value = rowLower[i] - rowActivity[i];
2422        if (value > primalTolerance) {
2423            numberBadRows++;
2424            largestInfeasibility = CoinMax(largestInfeasibility, value);
2425            sumInfeasibility += value;
2426        }
2427        value = rowActivity[i] - rowUpper[i];
2428        if (value > primalTolerance) {
2429            numberBadRows++;
2430            largestInfeasibility = CoinMax(largestInfeasibility, value);
2431            sumInfeasibility += value;
2432        }
2433    }
2434#ifdef JJF_ZERO
2435    if (largestInfeasibility > primalTolerance && numberBadRows*10 < numberRows) {
2436        // Can we improve by flipping
2437        for (int iPass = 0; iPass < 10; iPass++) {
2438            int numberColumns = solver->getNumCols();
2439            const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol();
2440            const double * element = matrixByCol->getElements();
2441            const int * row = matrixByCol->getIndices();
2442            const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
2443            const int * columnLength = matrixByCol->getVectorLengths();
2444            double oldSum = sumInfeasibility;
2445            // First improve by moving continuous ones
2446            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2447                if (!solver->isInteger(iColumn)) {
2448                    double solValue = solution[iColumn];
2449                    double thetaUp = columnUpper[iColumn] - solValue;
2450                    double improvementUp = 0.0;
2451                    if (thetaUp > primalTolerance) {
2452                        // can go up
2453                        for (CoinBigIndex j = columnStart[iColumn];
2454                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2455                            int iRow = row[j];
2456                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2457                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2458                            double el = element[j];
2459                            if (el > 0.0) {
2460                                // positive element
2461                                if (distanceUp > 0.0) {
2462                                    if (thetaUp*el > distanceUp)
2463                                        thetaUp = distanceUp / el;
2464                                } else {
2465                                    improvementUp -= el;
2466                                }
2467                                if (distanceDown > 0.0) {
2468                                    if (thetaUp*el > distanceDown)
2469                                        thetaUp = distanceDown / el;
2470                                    improvementUp += el;
2471                                }
2472                            } else {
2473                                // negative element
2474                                if (distanceDown < 0.0) {
2475                                    if (thetaUp*el < distanceDown)
2476                                        thetaUp = distanceDown / el;
2477                                } else {
2478                                    improvementUp += el;
2479                                }
2480                                if (distanceUp < 0.0) {
2481                                    if (thetaUp*el < distanceUp)
2482                                        thetaUp = distanceUp / el;
2483                                    improvementUp -= el;
2484                                }
2485                            }
2486                        }
2487                    }
2488                    double thetaDown = solValue - columnLower[iColumn];
2489                    double improvementDown = 0.0;
2490                    if (thetaDown > primalTolerance) {
2491                        // can go down
2492                        for (CoinBigIndex j = columnStart[iColumn];
2493                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2494                            int iRow = row[j];
2495                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2496                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2497                            double el = -element[j]; // not change in sign form up
2498                            if (el > 0.0) {
2499                                // positive element
2500                                if (distanceUp > 0.0) {
2501                                    if (thetaDown*el > distanceUp)
2502                                        thetaDown = distanceUp / el;
2503                                } else {
2504                                    improvementDown -= el;
2505                                }
2506                                if (distanceDown > 0.0) {
2507                                    if (thetaDown*el > distanceDown)
2508                                        thetaDown = distanceDown / el;
2509                                    improvementDown += el;
2510                                }
2511                            } else {
2512                                // negative element
2513                                if (distanceDown < 0.0) {
2514                                    if (thetaDown*el < distanceDown)
2515                                        thetaDown = distanceDown / el;
2516                                } else {
2517                                    improvementDown += el;
2518                                }
2519                                if (distanceUp < 0.0) {
2520                                    if (thetaDown*el < distanceUp)
2521                                        thetaDown = distanceUp / el;
2522                                    improvementDown -= el;
2523                                }
2524                            }
2525                        }
2526                        if (thetaUp < 1.0e-8)
2527                            improvementUp = 0.0;
2528                        if (thetaDown < 1.0e-8)
2529                            improvementDown = 0.0;
2530                        double theta;
2531                        if (improvementUp >= improvementDown) {
2532                            theta = thetaUp;
2533                        } else {
2534                            improvementUp = improvementDown;
2535                            theta = -thetaDown;
2536                        }
2537                        if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
2538                            // Could move
2539                            double oldSum = 0.0;
2540                            double newSum = 0.0;
2541                            solution[iColumn] += theta;
2542                            for (CoinBigIndex j = columnStart[iColumn];
2543                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2544                                int iRow = row[j];
2545                                double lower = rowLower[iRow];
2546                                double upper = rowUpper[iRow];
2547                                double value = rowActivity[iRow];
2548                                if (value > upper)
2549                                    oldSum += value - upper;
2550                                else if (value < lower)
2551                                    oldSum += lower - value;
2552                                value += theta * element[j];
2553                                rowActivity[iRow] = value;
2554                                if (value > upper)
2555                                    newSum += value - upper;
2556                                else if (value < lower)
2557                                    newSum += lower - value;
2558                            }
2559                            assert (newSum <= oldSum);
2560                            sumInfeasibility += newSum - oldSum;
2561                        }
2562                    }
2563                }
2564            }
2565            // Now flip some integers?
2566#ifdef JJF_ZERO
2567            for (i = 0; i < numberIntegers; i++) {
2568                int iColumn = integerVariable[i];
2569                double solValue = solution[iColumn];
2570                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
2571                double improvementUp = 0.0;
2572                if (columnUpper[iColumn] >= solValue + 1.0) {
2573                    // can go up
2574                    double oldSum = 0.0;
2575                    double newSum = 0.0;
2576                    for (CoinBigIndex j = columnStart[iColumn];
2577                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2578                        int iRow = row[j];
2579                        double lower = rowLower[iRow];
2580                        double upper = rowUpper[iRow];
2581                        double value = rowActivity[iRow];
2582                        if (value > upper)
2583                            oldSum += value - upper;
2584                        else if (value < lower)
2585                            oldSum += lower - value;
2586                        value += element[j];
2587                        if (value > upper)
2588                            newSum += value - upper;
2589                        else if (value < lower)
2590                            newSum += lower - value;
2591                    }
2592                    improvementUp = oldSum - newSum;
2593                }
2594                double improvementDown = 0.0;
2595                if (columnLower[iColumn] <= solValue - 1.0) {
2596                    // can go down
2597                    double oldSum = 0.0;
2598                    double newSum = 0.0;
2599                    for (CoinBigIndex j = columnStart[iColumn];
2600                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2601                        int iRow = row[j];
2602                        double lower = rowLower[iRow];
2603                        double upper = rowUpper[iRow];
2604                        double value = rowActivity[iRow];
2605                        if (value > upper)
2606                            oldSum += value - upper;
2607                        else if (value < lower)
2608                            oldSum += lower - value;
2609                        value -= element[j];
2610                        if (value > upper)
2611                            newSum += value - upper;
2612                        else if (value < lower)
2613                            newSum += lower - value;
2614                    }
2615                    improvementDown = oldSum - newSum;
2616                }
2617                double theta;
2618                if (improvementUp >= improvementDown) {
2619                    theta = 1.0;
2620                } else {
2621                    improvementUp = improvementDown;
2622                    theta = -1.0;
2623                }
2624                if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
2625                    // Could move
2626                    double oldSum = 0.0;
2627                    double newSum = 0.0;
2628                    solution[iColumn] += theta;
2629                    for (CoinBigIndex j = columnStart[iColumn];
2630                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2631                        int iRow = row[j];
2632                        double lower = rowLower[iRow];
2633                        double upper = rowUpper[iRow];
2634                        double value = rowActivity[iRow];
2635                        if (value > upper)
2636                            oldSum += value - upper;
2637                        else if (value < lower)
2638                            oldSum += lower - value;
2639                        value += theta * element[j];
2640                        rowActivity[iRow] = value;
2641                        if (value > upper)
2642                            newSum += value - upper;
2643                        else if (value < lower)
2644                            newSum += lower - value;
2645                    }
2646                    assert (newSum <= oldSum);
2647                    sumInfeasibility += newSum - oldSum;
2648                }
2649            }
2650#else
2651            int bestColumn = -1;
2652            double bestImprovement = primalTolerance;
2653            double theta = 0.0;
2654            for (i = 0; i < numberIntegers; i++) {
2655                int iColumn = integerVariable[i];
2656                double solValue = solution[iColumn];
2657                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
2658                double improvementUp = 0.0;
2659                if (columnUpper[iColumn] >= solValue + 1.0) {
2660                    // can go up
2661                    double oldSum = 0.0;
2662                    double newSum = 0.0;
2663                    for (CoinBigIndex j = columnStart[iColumn];
2664                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2665                        int iRow = row[j];
2666                        double lower = rowLower[iRow];
2667                        double upper = rowUpper[iRow];
2668                        double value = rowActivity[iRow];
2669                        if (value > upper)
2670                            oldSum += value - upper;
2671                        else if (value < lower)
2672                            oldSum += lower - value;
2673                        value += element[j];
2674                        if (value > upper)
2675                            newSum += value - upper;
2676                        else if (value < lower)
2677                            newSum += lower - value;
2678                    }
2679                    improvementUp = oldSum - newSum;
2680                }
2681                double improvementDown = 0.0;
2682                if (columnLower[iColumn] <= solValue - 1.0) {
2683                    // can go down
2684                    double oldSum = 0.0;
2685                    double newSum = 0.0;
2686                    for (CoinBigIndex j = columnStart[iColumn];
2687                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2688                        int iRow = row[j];
2689                        double lower = rowLower[iRow];
2690                        double upper = rowUpper[iRow];
2691                        double value = rowActivity[iRow];
2692                        if (value > upper)
2693                            oldSum += value - upper;
2694                        else if (value < lower)
2695                            oldSum += lower - value;
2696                        value -= element[j];
2697                        if (value > upper)
2698                            newSum += value - upper;
2699                        else if (value < lower)
2700                            newSum += lower - value;
2701                    }
2702                    improvementDown = oldSum - newSum;
2703                }
2704                double improvement = CoinMax(improvementUp, improvementDown);
2705                if (improvement > bestImprovement) {
2706                    bestImprovement = improvement;
2707                    bestColumn = iColumn;
2708                    if (improvementUp > improvementDown)
2709                        theta = 1.0;
2710                    else
2711                        theta = -1.0;
2712                }
2713            }
2714            if (bestColumn >= 0) {
2715                // Could move
2716                int iColumn = bestColumn;
2717                double oldSum = 0.0;
2718                double newSum = 0.0;
2719                solution[iColumn] += theta;
2720                for (CoinBigIndex j = columnStart[iColumn];
2721                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2722                    int iRow = row[j];
2723                    double lower = rowLower[iRow];
2724                    double upper = rowUpper[iRow];
2725                    double value = rowActivity[iRow];
2726                    if (value > upper)
2727                        oldSum += value - upper;
2728                    else if (value < lower)
2729                        oldSum += lower - value;
2730                    value += theta * element[j];
2731                    rowActivity[iRow] = value;
2732                    if (value > upper)
2733                        newSum += value - upper;
2734                    else if (value < lower)
2735                        newSum += lower - value;
2736                }
2737                assert (newSum <= oldSum);
2738                sumInfeasibility += newSum - oldSum;
2739            }
2740#endif
2741            if (oldSum <= sumInfeasibility + primalTolerance)
2742                break; // no good
2743        }
2744    }
2745    //delete [] saveSolution;
2746#endif
2747    delete [] rowActivity;
2748    return (largestInfeasibility > primalTolerance) ? 0 : 1;
2749}
2750// Set maximum Time (default off) - also sets starttime to current
2751void
2752CbcHeuristicFPump::setMaximumTime(double value)
2753{
2754    startTime_ = CoinCpuTime();
2755    maximumTime_ = value;
2756}
2757
2758# ifdef COIN_HAS_CLP
2759
2760//#############################################################################
2761// Constructors / Destructor / Assignment
2762//#############################################################################
2763
2764//-------------------------------------------------------------------
2765// Default Constructor
2766//-------------------------------------------------------------------
2767CbcDisasterHandler::CbcDisasterHandler (CbcModel * model)
2768        : OsiClpDisasterHandler(),
2769        cbcModel_(model)
2770{
2771    if (model) {
2772        osiModel_
2773        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2774        if (osiModel_)
2775            setSimplex(osiModel_->getModelPtr());
2776    }
2777}
2778
2779//-------------------------------------------------------------------
2780// Copy constructor
2781//-------------------------------------------------------------------
2782CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs)
2783        : OsiClpDisasterHandler(rhs),
2784        cbcModel_(rhs.cbcModel_)
2785{
2786}
2787
2788
2789//-------------------------------------------------------------------
2790// Destructor
2791//-------------------------------------------------------------------
2792CbcDisasterHandler::~CbcDisasterHandler ()
2793{
2794}
2795
2796//----------------------------------------------------------------
2797// Assignment operator
2798//-------------------------------------------------------------------
2799CbcDisasterHandler &
2800CbcDisasterHandler::operator=(const CbcDisasterHandler & rhs)
2801{
2802    if (this != &rhs) {
2803        OsiClpDisasterHandler::operator=(rhs);
2804        cbcModel_ = rhs.cbcModel_;
2805    }
2806    return *this;
2807}
2808//-------------------------------------------------------------------
2809// Clone
2810//-------------------------------------------------------------------
2811ClpDisasterHandler * CbcDisasterHandler::clone() const
2812{
2813    return new CbcDisasterHandler(*this);
2814}
2815// Type of disaster 0 can fix, 1 abort
2816int
2817CbcDisasterHandler::typeOfDisaster()
2818{
2819    if (!cbcModel_->parentModel() &&
2820            (cbcModel_->specialOptions()&2048) == 0) {
2821        return 0;
2822    } else {
2823        if (cbcModel_->parentModel())
2824            cbcModel_->setMaximumNodes(0);
2825        return 1;
2826    }
2827}
2828/* set model. */
2829void
2830CbcDisasterHandler::setCbcModel(CbcModel * model)
2831{
2832    cbcModel_ = model;
2833    if (model) {
2834        osiModel_
2835        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2836        if (osiModel_)
2837            setSimplex(osiModel_->getModelPtr());
2838        else
2839            setSimplex(NULL);
2840    }
2841}
2842#endif
2843
Note: See TracBrowser for help on using the repository browser.