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

Last change on this file since 2094 was 2094, checked in by forrest, 3 years ago

for memory leaks and heuristics and some experimental stuff

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