source: branches/sandbox/Cbc/src/CbcHeuristicFPump.cpp @ 1286

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

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