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

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

add random seed setting

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