source: stable/2.7/Cbc/src/CbcHeuristicFPump.cpp @ 1995

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

fix for minibab

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