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

Last change on this file since 2344 was 2344, checked in by forrest, 2 years ago

change int to CoinBigIndex?

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