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

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

for debug

  • 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 2434 2018-11-23 19:19:58Z 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()-model_->getCutoffIncrement();
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-model_->getCutoffIncrement();
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                    cutoff2 = CoinMin(cutoff2,realCutoff);
2059#if 0
2060                      {
2061                        OsiClpSolverInterface * clpSolver
2062                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
2063                        if (clpSolver) {
2064                            ClpSimplex * simplex = clpSolver->getModelPtr();
2065                            simplex->writeMps("testA.mps",2,1);
2066                        }
2067                      }
2068#endif
2069                    int returnCode2 = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
2070                                                          cutoff2, "CbcHeuristicLocalAfterFPump");
2071                    fractionSmall_ = saveFraction;
2072                    if (returnCode2 < 0) {
2073                        if (returnCode2 == -2) {
2074                            exitAll = true;
2075                            returnCode = 0;
2076                        } else {
2077                            returnCode2 = 0; // returned on size - try changing
2078                            //#define ROUND_AGAIN
2079#ifdef ROUND_AGAIN
2080                            if (numberTries == 1 && numberPasses > 20 && allowedPass < numberPasses - 1) {
2081                                allowedPass = (numberPasses + allowedPass) >> 1;
2082                                sprintf(pumpPrint,
2083                                        "Fixing all variables which were last changed on pass %d and trying again",
2084                                        allowedPass);
2085                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2086                                << pumpPrint
2087                                << CoinMessageEol;
2088                                stopBAB = false;
2089                                continue;
2090                            }
2091#endif
2092                        }
2093                    }
2094                    if ((returnCode2&2) != 0) {
2095                        // could add cut
2096                        returnCode2 &= ~2;
2097                    }
2098                    if (returnCode2) {
2099                        numberBandBsolutions++;
2100                        // may not have got solution earlier
2101                        returnCode |= 1;
2102                    }
2103                } else {
2104                    // no need
2105                    exitAll = true;
2106                    //returnCode=0;
2107                }
2108                // recompute solution value
2109                if (returnCode && true) {
2110#if 0
2111                      {
2112                        OsiClpSolverInterface * clpSolver
2113                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
2114                        if (clpSolver) {
2115                            ClpSimplex * simplex = clpSolver->getModelPtr();
2116                            simplex->writeMps("testB.mps",2,1);
2117                        }
2118                      }
2119#endif
2120                    delete newSolver;
2121                    newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
2122                    newSolutionValue = -saveOffset;
2123                    double newSumInfeas = 0.0;
2124                    const double * obj = newSolver->getObjCoefficients();
2125                    for (int i = 0 ; i < numberColumns ; i++ ) {
2126                        if (isHeuristicInteger(newSolver,i)) {
2127                            double value = newSolution[i];
2128                            double nearest = floor(value + 0.5);
2129                            newSumInfeas += fabs(value - nearest);
2130                        }
2131                        newSolutionValue += obj[i] * newSolution[i];
2132                    }
2133                    newSolutionValue *= direction;
2134                }
2135                bool gotSolution = false;
2136                if (returnCode && newSolutionValue < saveValue) {
2137                    sprintf(pumpPrint, "Mini branch and bound improved solution from %g to %g (%.2f seconds)",
2138                            saveValue, newSolutionValue, model_->getCurrentSeconds());
2139                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2140                    << pumpPrint
2141                    << CoinMessageEol;
2142                    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2143                    gotSolution = true;
2144                    if (fixContinuous && nFixC + nFixC2 > 0) {
2145                        // may be able to do even better
2146                        int nFixed = 0;
2147                        const double * lower = model_->solver()->getColLower();
2148                        const double * upper = model_->solver()->getColUpper();
2149                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2150                            double value = newSolution[iColumn];
2151                            if (isHeuristicInteger(newSolver,iColumn)) {
2152                                value = floor(newSolution[iColumn] + 0.5);
2153                                newSolver->setColLower(iColumn, value);
2154                                newSolver->setColUpper(iColumn, value);
2155                                nFixed++;
2156                            } else {
2157                                newSolver->setColLower(iColumn, lower[iColumn]);
2158                                newSolver->setColUpper(iColumn, upper[iColumn]);
2159                                if (value < lower[iColumn])
2160                                    value = lower[iColumn];
2161                                else if (value > upper[iColumn])
2162                                    value = upper[iColumn];
2163
2164                            }
2165                            newSolution[iColumn] = value;
2166                        }
2167                        newSolver->setColSolution(newSolution);
2168                        //#define CLP_INVESTIGATE2
2169#ifdef CLP_INVESTIGATE2
2170                        {
2171                            // check
2172                            // get row activities
2173                            int numberRows = newSolver->getNumRows();
2174                            double * rowActivity = new double[numberRows];
2175                            memset(rowActivity, 0, numberRows*sizeof(double));
2176                            newSolver->getMatrixByCol()->times(newSolution, rowActivity) ;
2177                            double largestInfeasibility = primalTolerance;
2178                            double sumInfeasibility = 0.0;
2179                            int numberBadRows = 0;
2180                            const double * rowLower = newSolver->getRowLower();
2181                            const double * rowUpper = newSolver->getRowUpper();
2182                            for (i = 0 ; i < numberRows ; i++) {
2183                                double value;
2184                                value = rowLower[i] - rowActivity[i];
2185                                if (value > primalTolerance) {
2186                                    numberBadRows++;
2187                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
2188                                    sumInfeasibility += value;
2189                                }
2190                                value = rowActivity[i] - rowUpper[i];
2191                                if (value > primalTolerance) {
2192                                    numberBadRows++;
2193                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
2194                                    sumInfeasibility += value;
2195                                }
2196                            }
2197                            printf("%d bad rows, largest inf %g sum %g\n",
2198                                   numberBadRows, largestInfeasibility, sumInfeasibility);
2199                            delete [] rowActivity;
2200                        }
2201#endif
2202#ifdef COIN_HAS_CLP
2203                        OsiClpSolverInterface * clpSolver
2204                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
2205                        if (clpSolver) {
2206                            ClpSimplex * simplex = clpSolver->getModelPtr();
2207                            //simplex->writeBasis("test.bas",true,2);
2208                            //simplex->writeMps("test.mps",2,1);
2209                            if (nFixed*3 > numberColumns*2)
2210                                simplex->allSlackBasis(); // may as well go from all slack
2211                            int logLevel=simplex->logLevel();
2212                            if (logLevel<=1)
2213                              simplex->setLogLevel(0);
2214                            simplex->primal(1);
2215                            simplex->setLogLevel(logLevel);
2216                            clpSolver->setWarmStart(NULL);
2217                        }
2218#endif
2219                        newSolver->initialSolve();
2220                        if (newSolver->isProvenOptimal()) {
2221                            double value = newSolver->getObjValue() * newSolver->getObjSense();
2222                            if (value < newSolutionValue) {
2223                              //newSolver->writeMpsNative("query.mps", NULL, NULL, 2);
2224#ifdef JJF_ZERO
2225                                {
2226                                    double saveOffset;
2227                                    newSolver->getDblParam(OsiObjOffset, saveOffset);
2228                                    const double * obj = newSolver->getObjCoefficients();
2229                                    double newTrueSolutionValue = -saveOffset;
2230                                    double newSumInfeas = 0.0;
2231                                    int numberColumns = newSolver->getNumCols();
2232                                    const double * solution = newSolver->getColSolution();
2233                                    for (int  i = 0 ; i < numberColumns ; i++ ) {
2234                                        if (isHeuristicInteger(newSolver,i)) {
2235                                            double value = solution[i];
2236                                            double nearest = floor(value + 0.5);
2237                                            newSumInfeas += fabs(value - nearest);
2238                                        }
2239                                        if (solution[i])
2240                                            printf("%d obj %g val %g - total %g\n", i, obj[i], solution[i],
2241                                                   newTrueSolutionValue);
2242                                        newTrueSolutionValue += obj[i] * solution[i];
2243                                    }
2244                                    printf("obj %g - inf %g\n", newTrueSolutionValue,
2245                                           newSumInfeas);
2246                                }
2247#endif
2248                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
2249                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2250                                << pumpPrint
2251                                << CoinMessageEol;
2252                                newSolutionValue = value;
2253                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
2254                            }
2255                        } else {
2256                            //newSolver->writeMps("bad3.mps");
2257                          sprintf(pumpPrint, "On closer inspection solution is not valid");
2258                          model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2259                            << pumpPrint
2260                            << CoinMessageEol;
2261                            exitAll = true;
2262                            break;
2263                        }
2264                    }
2265                } else {
2266                    sprintf(pumpPrint, "Mini branch and bound did not improve solution (%.2f seconds)",
2267                            model_->getCurrentSeconds());
2268                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2269                    << pumpPrint
2270                    << CoinMessageEol;
2271                    if (returnCode && newSolutionValue < saveValue + 1.0e-3 && nFixC + nFixC2) {
2272                        // may be able to do better
2273                        const double * lower = model_->solver()->getColLower();
2274                        const double * upper = model_->solver()->getColUpper();
2275                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2276                            if (isHeuristicInteger(newSolver,iColumn)) {
2277                                double value = floor(newSolution[iColumn] + 0.5);
2278                                newSolver->setColLower(iColumn, value);
2279                                newSolver->setColUpper(iColumn, value);
2280                            } else {
2281                                newSolver->setColLower(iColumn, lower[iColumn]);
2282                                newSolver->setColUpper(iColumn, upper[iColumn]);
2283                            }
2284                        }
2285                        newSolver->initialSolve();
2286                        if (newSolver->isProvenOptimal()) {
2287                            double value = newSolver->getObjValue() * newSolver->getObjSense();
2288                            if (value < saveValue) {
2289                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
2290                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2291                                << pumpPrint
2292                                << CoinMessageEol;
2293                                //newSolver->writeMpsNative("query2.mps", NULL, NULL, 2);
2294                                newSolutionValue = value;
2295                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
2296                                gotSolution = true;
2297                            }
2298                        }
2299                    }
2300                }
2301                if (gotSolution) {
2302                    if ((accumulate_&1) != 0) {
2303                        model_->incrementUsed(betterSolution); // for local search
2304                    }
2305                    solutionValue = newSolutionValue;
2306                    solutionFound = true;
2307                    numberSolutions++;
2308                    if (numberSolutions>=maxSolutions)
2309                      exitAll = true;
2310                    if (exitNow(newSolutionValue))
2311                        exitAll = true;
2312                    CoinWarmStartBasis * basis =
2313                        dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
2314                    if (basis) {
2315                        bestBasis = * basis;
2316                        delete basis;
2317                        int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
2318                        if (action == 0) {
2319                            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
2320                            double saveObjectiveValue = model_->getMinimizationObjValue();
2321                            model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
2322                            if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
2323                                model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
2324                            delete [] saveOldSolution;
2325                        }
2326                        if (!action || model_->getCurrentSeconds() > model_->getMaximumSeconds()) {
2327                            exitAll = true; // exit
2328                            break;
2329                        }
2330                    }
2331                }
2332            } // end stopBAB while
2333            delete newSolver;
2334        }
2335        if (solutionFound) finalReturnCode = 1;
2336        cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
2337        realCutoff = cutoff;
2338        if (numberTries >= maximumRetries_ || !solutionFound || exitAll || cutoff < continuousObjectiveValue + 1.0e-7) {
2339            break;
2340        } else {
2341            solutionFound = false;
2342            if (absoluteIncrement_ > 0.0 || relativeIncrement_ > 0.0) {
2343                double gap = relativeIncrement_ * fabs(solutionValue);
2344                double change = CoinMax(gap, absoluteIncrement_);
2345                cutoff = CoinMin(cutoff, solutionValue - change);
2346            } else {
2347                //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
2348                double weights[10] = {0.1, 0.2, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5, 0.6};
2349                cutoff -= weights[CoinMin(numberTries-1, 9)] * (cutoff - continuousObjectiveValue);
2350            }
2351            // But round down
2352            if (exactMultiple)
2353                cutoff = exactMultiple * floor(cutoff / exactMultiple);
2354            if (cutoff < continuousObjectiveValue)
2355                break;
2356            sprintf(pumpPrint, "Round again with cutoff of %g", cutoff);
2357            secondMajorPass=true;
2358            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2359            << pumpPrint
2360            << CoinMessageEol;
2361            if ((accumulate_&3) < 2 && usedColumn) {
2362                for (int i = 0; i < numberColumns; i++)
2363                    usedColumn[i] = -1;
2364            }
2365            averageIterationsPerTry = totalNumberIterations / numberTries;
2366            numberIterationsLastPass =  totalNumberIterations;
2367            totalNumberPasses += numberPasses - 1;
2368        }
2369    }
2370/*
2371  End of the `exitAll' loop.
2372*/
2373#ifdef RAND_RAND
2374    delete [] randomFactor;
2375#endif
2376    delete solver; // probably NULL but do anyway
2377    if (!finalReturnCode && closestSolution && closestObjectiveValue <= 10.0 &&
2378            usedColumn && !model_->maximumSecondsReached()) {
2379        // try a bit of branch and bound
2380        OsiSolverInterface * newSolver = cloneBut(1); // was model_->continuousSolver()->clone();
2381        const double * colLower = newSolver->getColLower();
2382        const double * colUpper = newSolver->getColUpper();
2383        int i;
2384        double rhs = 0.0;
2385        for (i = 0; i < numberIntegers; i++) {
2386            int iColumn = integerVariable[i];
2387            int direction = closestSolution[i];
2388            closestSolution[i] = iColumn;
2389            if (direction == 0) {
2390                // keep close to LB
2391                rhs += colLower[iColumn];
2392                lastSolution[i] = 1.0;
2393            } else {
2394                // keep close to UB
2395                rhs -= colUpper[iColumn];
2396                lastSolution[i] = -1.0;
2397            }
2398        }
2399        newSolver->addRow(numberIntegers, closestSolution,
2400                          lastSolution, -COIN_DBL_MAX, rhs + 10.0);
2401        //double saveValue = newSolutionValue;
2402        //newSolver->writeMps("sub");
2403        int returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
2404                                             newSolutionValue, "CbcHeuristicLocalAfterFPump");
2405        if (returnCode < 0)
2406            returnCode = 0; // returned on size
2407        if ((returnCode&2) != 0) {
2408            // could add cut
2409            returnCode &= ~2;
2410        }
2411        if (returnCode) {
2412            //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
2413            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2414            //abort();
2415            solutionValue = newSolutionValue;
2416            solutionFound = true;
2417            numberSolutions++;
2418            if (numberSolutions>=maxSolutions)
2419              exitAll = true;
2420            if (exitNow(newSolutionValue))
2421                exitAll = true;
2422        }
2423        delete newSolver;
2424    }
2425    delete clonedSolver;
2426    delete [] roundingSolution;
2427    delete [] usedColumn;
2428    delete [] lastSolution;
2429    delete [] newSolution;
2430    delete [] closestSolution;
2431    delete [] integerVariable;
2432    delete [] firstPerturbedObjective;
2433    delete [] firstPerturbedSolution;
2434    if (solutionValue == incomingObjective)
2435        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
2436                model_->getCurrentSeconds(), CoinCpuTime() - time1);
2437    else if (numberSolutions < maxSolutions)
2438      sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
2439              model_->getCurrentSeconds(), solutionValue, CoinCpuTime() - time1);
2440    else 
2441        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g (stopping after %d solutions) - took %.2f seconds",
2442                model_->getCurrentSeconds(), solutionValue, 
2443                numberSolutions,CoinCpuTime() - time1);
2444    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2445      << pumpPrint
2446      << CoinMessageEol;
2447    if (bestBasis.getNumStructural())
2448        model_->setBestSolutionBasis(bestBasis);
2449    //model_->setMinimizationObjValue(saveBestObjective);
2450    if ((accumulate_&1) != 0 && numberSolutions > 1 && !model_->getSolutionCount()) {
2451        model_->setSolutionCount(1); // for local search
2452        model_->setNumberHeuristicSolutions(1);
2453    }
2454#ifdef COIN_DEVELOP
2455    {
2456        double ncol = model_->solver()->getNumCols();
2457        double nrow = model_->solver()->getNumRows();
2458        printf("XXX total iterations %d ratios - %g %g %g\n",
2459               totalNumberIterations,
2460               static_cast<double> (totalNumberIterations) / nrow,
2461               static_cast<double> (totalNumberIterations) / ncol,
2462               static_cast<double> (totalNumberIterations) / (2*nrow + 2*ncol));
2463    }
2464#endif
2465    return finalReturnCode;
2466}
2467
2468/**************************END MAIN PROCEDURE ***********************************/
2469/* If general integers then adds variables to turn into binaries round
2470   solution
2471*/
2472int 
2473CbcHeuristicFPump::solution(double & objectiveValue, double * newSolution)
2474{
2475  double * newSolution2=NULL;
2476  double objective2=COIN_DBL_MAX;
2477  int returnCode2=0;
2478  int oddGeneral = (accumulate_&(32|64|128))>>5;
2479  if (oddGeneral) {
2480    int maxAround=2;
2481    bool fixSatisfied=false;
2482    // clone solver and modify
2483    OsiSolverInterface * solver = cloneBut(2); // wasmodel_->solver()->clone();
2484    double cutoff;
2485    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
2486    int numberColumns = model_->getNumCols();
2487    int numberIntegers = model_->numberIntegers();
2488    const int * integerVariableOrig = model_->integerVariable();
2489    int general = 0;
2490    int nAdd=0;
2491    const double * lower = solver->getColLower();
2492    const double * upper = solver->getColUpper();
2493    const double * initialSolution = solver->getColSolution();
2494    // we may be being clever so make sure solver lines up wuth model
2495    for (int i = 0; i < numberColumns; i++) 
2496      solver->setContinuous(i);
2497    for (int i = 0; i < numberIntegers; i++) {
2498      int iColumn = integerVariableOrig[i];
2499#ifdef COIN_HAS_CLP
2500      if (!isHeuristicInteger(solver,iColumn))
2501        continue;
2502#endif
2503      double value = initialSolution[iColumn];
2504      double nearest = floor(value + 0.5);
2505      if (upper[iColumn] - lower[iColumn] > 1.000001) {
2506        if (fabs(value-nearest)<1.0e-6&&fixSatisfied) {
2507          solver->setColLower(iColumn,nearest);
2508          solver->setColUpper(iColumn,nearest);
2509        } else {
2510          general++;
2511          int up=static_cast<int>(upper[iColumn]);
2512          int lo=static_cast<int>(lower[iColumn]);
2513          int near=static_cast<int>(nearest);
2514          up = CoinMin(up,near+maxAround);
2515          lo = CoinMax(lo,near-maxAround);
2516          solver->setColLower(iColumn,lo);
2517          solver->setColUpper(iColumn,up);
2518          int n=up-lo;
2519          // 1 - 1, 2,3 - 2, 4567 - 3
2520          while (n) {
2521            nAdd++;
2522            n = n>>1;
2523          }
2524        }
2525      } else {
2526        solver->setInteger(iColumn);
2527      }
2528    }
2529    if (!general) {
2530      delete solver;
2531    }
2532    if (general) {
2533      CbcModel * saveModel = model_;
2534      CoinBigIndex * addStart = new CoinBigIndex[nAdd+1];
2535      memset(addStart,0,(nAdd+1)*sizeof(CoinBigIndex));
2536      int * addIndex = new int[general+nAdd];
2537      double * addElement = new double[general+nAdd];
2538      double * addLower = new double[nAdd];
2539      double * addUpper = new double[nAdd];
2540      for (int i=0;i<nAdd;i++) {
2541        addLower[i]=0.0;
2542        addUpper[i]=1.0;
2543      }
2544      solver->addCols(nAdd, addStart, NULL, NULL, addLower, addUpper, NULL);
2545      lower = solver->getColLower();
2546      upper = solver->getColUpper();
2547      // now rows
2548      nAdd = 0;
2549      int nEl = 0;
2550      int nAddRow = 0;
2551      for (int i = 0; i < numberIntegers; i++) {
2552        int iColumn = integerVariableOrig[i];
2553#ifdef COIN_HAS_CLP
2554        if (!isHeuristicInteger(solver,iColumn))
2555          continue;
2556#endif
2557        if (upper[iColumn] - lower[iColumn] > 1.000001) {
2558          int up=static_cast<int>(upper[iColumn]);
2559          int lo=static_cast<int>(lower[iColumn]);
2560          addLower[nAddRow] = lo;
2561          addUpper[nAddRow] = lo;
2562          addIndex[nEl] = iColumn;
2563          addElement[nEl++] = 1.0;
2564          int n=up-lo;
2565          // 1 - 1, 2,3 - 2, 4567 - 3
2566          int el=1;
2567          while (n) {
2568            addIndex[nEl] = numberColumns+nAdd;
2569            addElement[nEl++] = -el;
2570            nAdd++;
2571            n = n>>1;
2572            el = el<<1;
2573          }
2574          nAddRow++;
2575          addStart[nAddRow]=nEl;
2576        }
2577      }
2578      for (int i=0;i<nAdd;i++)
2579        solver->setInteger(i+numberColumns);
2580      solver->addRows(nAddRow, addStart, addIndex, addElement, addLower, addUpper);
2581      delete [] addStart;
2582      delete [] addIndex;
2583      delete [] addElement;
2584      delete [] addLower;
2585      delete [] addUpper;
2586      solver->resolve();
2587      solver->writeMps("test");
2588      // new CbcModel
2589      model_=new CbcModel(*solver);
2590      model_->findIntegers(true);
2591      // set cutoff
2592      solver->setDblParam(OsiDualObjectiveLimit, cutoff);
2593      model_->setCutoff(cutoff);
2594      newSolution2 = new double [numberColumns+nAdd]; 
2595      objective2=objectiveValue;
2596      returnCode2=solutionInternal(objective2,newSolution2);
2597      delete solver;
2598      delete model_;
2599      model_=saveModel;
2600    }
2601  }
2602  int returnCode=solutionInternal(objectiveValue,newSolution);
2603  if (returnCode2&&false) {
2604    int numberColumns=model_->getNumCols();
2605    memcpy(newSolution,newSolution2,numberColumns*sizeof(double));
2606  }
2607  delete [] newSolution2;
2608  return returnCode;
2609}
2610
2611// update model
2612void CbcHeuristicFPump::setModel(CbcModel * model)
2613{
2614    model_ = model;
2615}
2616
2617/* Rounds solution - down if < downValue
2618   returns 1 if current is a feasible solution
2619*/
2620int
2621CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution,
2622                          //const double * objective,
2623                          int numberIntegers, const int * integerVariable,
2624                          /*char * pumpPrint,*/ int iter,
2625                          /*bool roundExpensive,*/ double downValue, int *flip)
2626{
2627    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2628    double primalTolerance ;
2629    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
2630
2631    int i;
2632
2633    const double * cost = solver->getObjCoefficients();
2634    int flip_up = 0;
2635    int flip_down  = 0;
2636    double  v = randomNumberGenerator_.randomDouble() * 20.0;
2637    int nn = 10 + static_cast<int> (v);
2638    int nnv = 0;
2639    int * list = new int [nn];
2640    double * val = new double [nn];
2641    for (i = 0; i < nn; i++) val[i] = .001;
2642
2643    const double * rowLower = solver->getRowLower();
2644    const double * rowUpper = solver->getRowUpper();
2645    int numberRows = solver->getNumRows();
2646    if (false && (iter&1) != 0) {
2647        // Do set covering variables
2648        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2649        const double * elementByRow = matrixByRow->getElements();
2650        const int * column = matrixByRow->getIndices();
2651        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2652        const int * rowLength = matrixByRow->getVectorLengths();
2653        for (i = 0; i < numberRows; i++) {
2654            if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2655                bool cover = true;
2656                double largest = 0.0;
2657                int jColumn = -1;
2658                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2659                    int iColumn = column[k];
2660                    if (elementByRow[k] != 1.0 || !isHeuristicInteger(solver,iColumn)) {
2661                        cover = false;
2662                        break;
2663                    } else {
2664                        if (solution[iColumn]) {
2665                            double value = solution[iColumn] *
2666                                           (randomNumberGenerator_.randomDouble() + 5.0);
2667                            if (value > largest) {
2668                                largest = value;
2669                                jColumn = iColumn;
2670                            }
2671                        }
2672                    }
2673                }
2674                if (cover) {
2675                    for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2676                        int iColumn = column[k];
2677                        if (iColumn == jColumn)
2678                            solution[iColumn] = 1.0;
2679                        else
2680                            solution[iColumn] = 0.0;
2681                    }
2682                }
2683            }
2684        }
2685    }
2686    int numberColumns = solver->getNumCols();
2687#ifdef JJF_ZERO
2688    // Do set covering variables
2689    const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2690    const double * elementByRow = matrixByRow->getElements();
2691    const int * column = matrixByRow->getIndices();
2692    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2693    const int * rowLength = matrixByRow->getVectorLengths();
2694    double * sortTemp = new double[numberColumns];
2695    int * whichTemp = new int [numberColumns];
2696    char * rowTemp = new char [numberRows];
2697    memset(rowTemp, 0, numberRows);
2698    for (i = 0; i < numberColumns; i++)
2699        whichTemp[i] = -1;
2700    int nSOS = 0;
2701    for (i = 0; i < numberRows; i++) {
2702        if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2703            bool cover = true;
2704            for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2705                int iColumn = column[k];
2706                if (elementByRow[k] != 1.0 || !isHeuristicInteger(solver,iColumn)) {
2707                    cover = false;
2708                    break;
2709                }
2710            }
2711            if (cover) {
2712                rowTemp[i] = 1;
2713                nSOS++;
2714                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2715                    int iColumn = column[k];
2716                    double value = solution[iColumn];
2717                    whichTemp[iColumn] = iColumn;
2718                }
2719            }
2720        }
2721    }
2722    if (nSOS) {
2723        // Column copy
2724        const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol();
2725        //const double * element = matrixByColumn->getElements();
2726        const int * row = matrixByColumn->getIndices();
2727        const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
2728        const int * columnLength = matrixByColumn->getVectorLengths();
2729        int nLook = 0;
2730        for (i = 0; i < numberColumns; i++) {
2731            if (whichTemp[i] >= 0) {
2732                whichTemp[nLook] = i;
2733                double value = solution[i];
2734                if (value < 0.5)
2735                    value *= (0.1 * randomNumberGenerator_.randomDouble() + 0.3);
2736                sortTemp[nLook++] = -value;
2737            }
2738        }
2739        CoinSort_2(sortTemp, sortTemp + nLook, whichTemp);
2740        double smallest = 1.0;
2741        int nFix = 0;
2742        int nOne = 0;
2743        for (int j = 0; j < nLook; j++) {
2744            int jColumn = whichTemp[j];
2745            double thisValue = solution[jColumn];
2746            if (!thisValue)
2747                continue;
2748            if (thisValue == 1.0)
2749                nOne++;
2750            smallest = CoinMin(smallest, thisValue);
2751            solution[jColumn] = 1.0;
2752            double largest = 0.0;
2753            for (CoinBigIndex jEl = columnStart[jColumn];
2754                    jEl < columnStart[jColumn] + columnLength[jColumn]; jEl++) {
2755                int jRow = row[jEl];
2756                if (rowTemp[jRow]) {
2757                    for (CoinBigIndex k = rowStart[jRow]; k < rowStart[jRow] + rowLength[jRow]; k++) {
2758                        int iColumn = column[k];
2759                        if (solution[iColumn]) {
2760                            if (iColumn != jColumn) {
2761                                double value = solution[iColumn];
2762                                if (value > largest)
2763                                    largest = value;
2764                                solution[iColumn] = 0.0;
2765                            }
2766                        }
2767                    }
2768                }
2769            }
2770            if (largest > thisValue)
2771                printf("%d was at %g - chosen over a value of %g\n",
2772                       jColumn, thisValue, largest);
2773            nFix++;
2774        }
2775        printf("%d fixed out of %d (%d at one already)\n",
2776               nFix, nLook, nOne);
2777    }
2778    delete [] sortTemp;
2779    delete [] whichTemp;
2780    delete [] rowTemp;
2781#endif
2782    const double * columnLower = solver->getColLower();
2783    const double * columnUpper = solver->getColUpper();
2784    // Check if valid with current solution (allow for 0.99999999s)
2785    double newSumInfeas = 0.0;
2786    int newNumberInfeas = 0;
2787    for (i = 0; i < numberIntegers; i++) {
2788        int iColumn = integerVariable[i];
2789        double value = solution[iColumn];
2790        double round = floor(value + 0.5);
2791        if (fabs(value - round) > primalTolerance) {
2792          newSumInfeas += fabs(value-round);
2793          newNumberInfeas++;
2794        }
2795    }
2796    if (!newNumberInfeas) {
2797        // may be able to use solution even if 0.99999's
2798        double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
2799        double * saveUpper = CoinCopyOfArray(columnUpper, numberColumns);
2800        double * saveSolution = CoinCopyOfArray(solution, numberColumns);
2801        double * tempSolution = CoinCopyOfArray(solution, numberColumns);
2802        CoinWarmStartBasis  * saveBasis =
2803            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
2804        for (i = 0; i < numberIntegers; i++) {
2805            int iColumn = integerVariable[i];
2806            double value = solution[iColumn];
2807            double round = floor(value + 0.5);
2808            solver->setColLower(iColumn, round);
2809            solver->setColUpper(iColumn, round);
2810            tempSolution[iColumn] = round;
2811        }
2812        solver->setColSolution(tempSolution);
2813        delete [] tempSolution;
2814        solver->resolve();
2815        solver->setColLower(saveLower);
2816        solver->setColUpper(saveUpper);
2817        solver->setWarmStart(saveBasis);
2818        delete [] saveLower;
2819        delete [] saveUpper;
2820        delete saveBasis;
2821        if (!solver->isProvenOptimal()) {
2822            solver->setColSolution(saveSolution);
2823        }
2824        delete [] saveSolution;
2825        if (solver->isProvenOptimal()) {
2826            // feasible
2827            delete [] list;
2828            delete [] val;
2829            return 1;
2830        }
2831    }
2832    //double * saveSolution = CoinCopyOfArray(solution,numberColumns);
2833    // return rounded solution
2834    for (i = 0; i < numberIntegers; i++) {
2835        int iColumn = integerVariable[i];
2836        double value = solution[iColumn];
2837        double round = floor(value + primalTolerance);
2838        if (value - round > downValue) round += 1.;
2839#ifndef JJF_ONE
2840        if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2841        if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2842#else
2843        if (round < columnLower[iColumn] + integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2844        if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2845#endif
2846        if (flip_up + flip_down == 0) {
2847            for (int k = 0; k < nn; k++) {
2848                if (fabs(value - round) > val[k]) {
2849                    nnv++;
2850                    for (int j = nn - 2; j >= k; j--) {
2851                        val[j+1] = val[j];
2852                        list[j+1] = list[j];
2853                    }
2854                    val[k] = fabs(value - round);
2855                    list[k] = iColumn;
2856                    break;
2857                }
2858            }
2859        }
2860        solution[iColumn] = round;
2861    }
2862
2863    if (nnv > nn) nnv = nn;
2864    //if (iter != 0)
2865    //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
2866    *flip = flip_up + flip_down;
2867
2868    if (*flip == 0 && iter != 0) {
2869        //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
2870        for (i = 0; i < nnv; i++) {
2871            // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
2872            int index = list[i];
2873            double value = solution[index];
2874            if (value <= 1.0) {
2875                solution[index] = 1.0 - value;
2876            } else if (value < columnLower[index] + integerTolerance) {
2877                solution[index] = value + 1.0;
2878            } else if (value > columnUpper[index] - integerTolerance) {
2879                solution[index] = value - 1.0;
2880            } else {
2881                solution[index] = value - 1.0;
2882            }
2883        }
2884        *flip = nnv;
2885    } else {
2886        //sprintf(pumpPrint+strlen(pumpPrint)," ");
2887    }
2888    delete [] list;
2889    delete [] val;
2890    //iter++;
2891
2892    // get row activities
2893    double * rowActivity = new double[numberRows];
2894    memset(rowActivity, 0, numberRows*sizeof(double));
2895    solver->getMatrixByCol()->times(solution, rowActivity) ;
2896    double largestInfeasibility = primalTolerance;
2897    double sumInfeasibility = 0.0;
2898    int numberBadRows = 0;
2899    for (i = 0 ; i < numberRows ; i++) {
2900        double value;
2901        value = rowLower[i] - rowActivity[i];
2902        if (value > primalTolerance) {
2903            numberBadRows++;
2904            largestInfeasibility = CoinMax(largestInfeasibility, value);
2905            sumInfeasibility += value;
2906        }
2907        value = rowActivity[i] - rowUpper[i];
2908        if (value > primalTolerance) {
2909            numberBadRows++;
2910            largestInfeasibility = CoinMax(largestInfeasibility, value);
2911            sumInfeasibility += value;
2912        }
2913    }
2914#ifdef JJF_ZERO
2915    if (largestInfeasibility > primalTolerance && numberBadRows*10 < numberRows) {
2916        // Can we improve by flipping
2917        for (int iPass = 0; iPass < 10; iPass++) {
2918            int numberColumns = solver->getNumCols();
2919            const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol();
2920            const double * element = matrixByCol->getElements();
2921            const int * row = matrixByCol->getIndices();
2922            const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
2923            const int * columnLength = matrixByCol->getVectorLengths();
2924            double oldSum = sumInfeasibility;
2925            // First improve by moving continuous ones
2926            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2927                if (!isHeuristicInteger(solver,iColumn)) {
2928                    double solValue = solution[iColumn];
2929                    double thetaUp = columnUpper[iColumn] - solValue;
2930                    double improvementUp = 0.0;
2931                    if (thetaUp > primalTolerance) {
2932                        // can go up
2933                        for (CoinBigIndex j = columnStart[iColumn];
2934                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2935                            int iRow = row[j];
2936                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2937                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2938                            double el = element[j];
2939                            if (el > 0.0) {
2940                                // positive element
2941                                if (distanceUp > 0.0) {
2942                                    if (thetaUp*el > distanceUp)
2943                                        thetaUp = distanceUp / el;
2944                                } else {
2945                                    improvementUp -= el;
2946                                }
2947                                if (distanceDown > 0.0) {
2948                                    if (thetaUp*el > distanceDown)
2949                                        thetaUp = distanceDown / el;
2950                                    improvementUp += el;
2951                                }
2952                            } else {
2953                                // negative element
2954                                if (distanceDown < 0.0) {
2955                                    if (thetaUp*el < distanceDown)
2956                                        thetaUp = distanceDown / el;
2957                                } else {
2958                                    improvementUp += el;
2959                                }
2960                                if (distanceUp < 0.0) {
2961                                    if (thetaUp*el < distanceUp)
2962                                        thetaUp = distanceUp / el;
2963                                    improvementUp -= el;
2964                                }
2965                            }
2966                        }
2967                    }
2968                    double thetaDown = solValue - columnLower[iColumn];
2969                    double improvementDown = 0.0;
2970                    if (thetaDown > primalTolerance) {
2971                        // can go down
2972                        for (CoinBigIndex j = columnStart[iColumn];
2973                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2974                            int iRow = row[j];
2975                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2976                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2977                            double el = -element[j]; // not change in sign form up
2978                            if (el > 0.0) {
2979                                // positive element
2980                                if (distanceUp > 0.0) {
2981                                    if (thetaDown*el > distanceUp)
2982                                        thetaDown = distanceUp / el;
2983                                } else {
2984                                    improvementDown -= el;
2985                                }
2986                                if (distanceDown > 0.0) {
2987                                    if (thetaDown*el > distanceDown)
2988                                        thetaDown = distanceDown / el;
2989                                    improvementDown += el;
2990                                }
2991                            } else {
2992                                // negative element
2993                                if (distanceDown < 0.0) {
2994                                    if (thetaDown*el < distanceDown)
2995                                        thetaDown = distanceDown / el;
2996                                } else {
2997                                    improvementDown += el;
2998                                }
2999                                if (distanceUp < 0.0) {
3000                                    if (thetaDown*el < distanceUp)
3001                                        thetaDown = distanceUp / el;
3002                                    improvementDown -= el;
3003                                }
3004                            }
3005                        }
3006                        if (thetaUp < 1.0e-8)
3007                            improvementUp = 0.0;
3008                        if (thetaDown < 1.0e-8)
3009                            improvementDown = 0.0;
3010                        double theta;
3011                        if (improvementUp >= improvementDown) {
3012                            theta = thetaUp;
3013                        } else {
3014                            improvementUp = improvementDown;
3015                            theta = -thetaDown;
3016                        }
3017                        if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
3018                            // Could move
3019                            double oldSum = 0.0;
3020                            double newSum = 0.0;
3021                            solution[iColumn] += theta;
3022                            for (CoinBigIndex j = columnStart[iColumn];
3023                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3024                                int iRow = row[j];
3025                                double lower = rowLower[iRow];
3026                                double upper = rowUpper[iRow];
3027                                double value = rowActivity[iRow];
3028                                if (value > upper)
3029                                    oldSum += value - upper;
3030                                else if (value < lower)
3031                                    oldSum += lower - value;
3032                                value += theta * element[j];
3033                                rowActivity[iRow] = value;
3034                                if (value > upper)
3035                                    newSum += value - upper;
3036                                else if (value < lower)
3037                                    newSum += lower - value;
3038                            }
3039                            assert (newSum <= oldSum);
3040                            sumInfeasibility += newSum - oldSum;
3041                        }
3042                    }
3043                }
3044            }
3045            // Now flip some integers?
3046#ifdef JJF_ZERO
3047            for (i = 0; i < numberIntegers; i++) {
3048                int iColumn = integerVariable[i];
3049                double solValue = solution[iColumn];
3050                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
3051                double improvementUp = 0.0;
3052                if (columnUpper[iColumn] >= solValue + 1.0) {
3053                    // can go up
3054                    double oldSum = 0.0;
3055                    double newSum = 0.0;
3056                    for (CoinBigIndex j = columnStart[iColumn];
3057                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3058                        int iRow = row[j];
3059                        double lower = rowLower[iRow];
3060                        double upper = rowUpper[iRow];
3061                        double value = rowActivity[iRow];
3062                        if (value > upper)
3063                            oldSum += value - upper;
3064                        else if (value < lower)
3065                            oldSum += lower - value;
3066                        value += element[j];
3067                        if (value > upper)
3068                            newSum += value - upper;
3069                        else if (value < lower)
3070                            newSum += lower - value;
3071                    }
3072                    improvementUp = oldSum - newSum;
3073                }
3074                double improvementDown = 0.0;
3075                if (columnLower[iColumn] <= solValue - 1.0) {
3076                    // can go down
3077                    double oldSum = 0.0;
3078                    double newSum = 0.0;
3079                    for (CoinBigIndex j = columnStart[iColumn];
3080                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3081                        int iRow = row[j];
3082                        double lower = rowLower[iRow];
3083                        double upper = rowUpper[iRow];
3084                        double value = rowActivity[iRow];
3085                        if (value > upper)
3086                            oldSum += value - upper;
3087                        else if (value < lower)
3088                            oldSum += lower - value;
3089                        value -= element[j];
3090                        if (value > upper)
3091                            newSum += value - upper;
3092                        else if (value < lower)
3093                            newSum += lower - value;
3094                    }
3095                    improvementDown = oldSum - newSum;
3096                }
3097                double theta;
3098                if (improvementUp >= improvementDown) {
3099                    theta = 1.0;
3100                } else {
3101                    improvementUp = improvementDown;
3102                    theta = -1.0;
3103                }
3104                if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
3105                    // Could move
3106                    double oldSum = 0.0;
3107                    double newSum = 0.0;
3108                    solution[iColumn] += theta;
3109                    for (CoinBigIndex j = columnStart[iColumn];
3110                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3111                        int iRow = row[j];
3112                        double lower = rowLower[iRow];
3113                        double upper = rowUpper[iRow];
3114                        double value = rowActivity[iRow];
3115                        if (value > upper)
3116                            oldSum += value - upper;
3117                        else if (value < lower)
3118                            oldSum += lower - value;
3119                        value += theta * element[j];
3120                        rowActivity[iRow] = value;
3121                        if (value > upper)
3122                            newSum += value - upper;
3123                        else if (value < lower)
3124                            newSum += lower - value;
3125                    }
3126                    assert (newSum <= oldSum);
3127                    sumInfeasibility += newSum - oldSum;
3128                }
3129            }
3130#else
3131            int bestColumn = -1;
3132            double bestImprovement = primalTolerance;
3133            double theta = 0.0;
3134            for (i = 0; i < numberIntegers; i++) {
3135                int iColumn = integerVariable[i];
3136                double solValue = solution[iColumn];
3137                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
3138                double improvementUp = 0.0;
3139                if (columnUpper[iColumn] >= solValue + 1.0) {
3140                    // can go up
3141                    double oldSum = 0.0;
3142                    double newSum = 0.0;
3143                    for (CoinBigIndex j = columnStart[iColumn];
3144                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3145                        int iRow = row[j];
3146                        double lower = rowLower[iRow];
3147                        double upper = rowUpper[iRow];
3148                        double value = rowActivity[iRow];
3149                        if (value > upper)
3150                            oldSum += value - upper;
3151                        else if (value < lower)
3152                            oldSum += lower - value;
3153                        value += element[j];
3154                        if (value > upper)
3155                            newSum += value - upper;
3156                        else if (value < lower)
3157                            newSum += lower - value;
3158                    }
3159                    improvementUp = oldSum - newSum;
3160                }
3161                double improvementDown = 0.0;
3162                if (columnLower[iColumn] <= solValue - 1.0) {
3163                    // can go down
3164                    double oldSum = 0.0;
3165                    double newSum = 0.0;
3166                    for (CoinBigIndex j = columnStart[iColumn];
3167                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3168                        int iRow = row[j];
3169                        double lower = rowLower[iRow];
3170                        double upper = rowUpper[iRow];
3171                        double value = rowActivity[iRow];
3172                        if (value > upper)
3173                            oldSum += value - upper;
3174                        else if (value < lower)
3175                            oldSum += lower - value;
3176                        value -= element[j];
3177                        if (value > upper)
3178                            newSum += value - upper;
3179                        else if (value < lower)
3180                            newSum += lower - value;
3181                    }
3182                    improvementDown = oldSum - newSum;
3183                }
3184                double improvement = CoinMax(improvementUp, improvementDown);
3185                if (improvement > bestImprovement) {
3186                    bestImprovement = improvement;
3187                    bestColumn = iColumn;
3188                    if (improvementUp > improvementDown)
3189                        theta = 1.0;
3190                    else
3191                        theta = -1.0;
3192                }
3193            }
3194            if (bestColumn >= 0) {
3195                // Could move
3196                int iColumn = bestColumn;
3197                double oldSum = 0.0;
3198                double newSum = 0.0;
3199                solution[iColumn] += theta;
3200                for (CoinBigIndex j = columnStart[iColumn];
3201                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3202                    int iRow = row[j];
3203                    double lower = rowLower[iRow];
3204                    double upper = rowUpper[iRow];
3205                    double value = rowActivity[iRow];
3206                    if (value > upper)
3207                        oldSum += value - upper;
3208                    else if (value < lower)
3209                        oldSum += lower - value;
3210                    value += theta * element[j];
3211                    rowActivity[iRow] = value;
3212                    if (value > upper)
3213                        newSum += value - upper;
3214                    else if (value < lower)
3215                        newSum += lower - value;
3216                }
3217                assert (newSum <= oldSum);
3218                sumInfeasibility += newSum - oldSum;
3219            }
3220#endif
3221            if (oldSum <= sumInfeasibility + primalTolerance)
3222                break; // no good
3223        }
3224    }
3225    //delete [] saveSolution;
3226#endif
3227    delete [] rowActivity;
3228    return (largestInfeasibility > primalTolerance) ? 0 : 1;
3229}
3230// Set maximum Time (default off) - also sets starttime to current
3231void
3232CbcHeuristicFPump::setMaximumTime(double value)
3233{
3234    startTime_ = CoinCpuTime();
3235    maximumTime_ = value;
3236}
3237
3238# ifdef COIN_HAS_CLP
3239
3240//#############################################################################
3241// Constructors / Destructor / Assignment
3242//#############################################################################
3243
3244//-------------------------------------------------------------------
3245// Default Constructor
3246//-------------------------------------------------------------------
3247CbcDisasterHandler::CbcDisasterHandler (CbcModel * model)
3248        : OsiClpDisasterHandler(),
3249        cbcModel_(model)
3250{
3251    if (model) {
3252        osiModel_
3253        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
3254        if (osiModel_)
3255            setSimplex(osiModel_->getModelPtr());
3256    }
3257}
3258
3259//-------------------------------------------------------------------
3260// Copy constructor
3261//-------------------------------------------------------------------
3262CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs)
3263        : OsiClpDisasterHandler(rhs),
3264        cbcModel_(rhs.cbcModel_)
3265{
3266}
3267
3268
3269//-------------------------------------------------------------------
3270// Destructor
3271//-------------------------------------------------------------------
3272CbcDisasterHandler::~CbcDisasterHandler ()
3273{
3274}
3275
3276//----------------------------------------------------------------
3277// Assignment operator
3278//-------------------------------------------------------------------
3279CbcDisasterHandler &
3280CbcDisasterHandler::operator=(const CbcDisasterHandler & rhs)
3281{
3282    if (this != &rhs) {
3283        OsiClpDisasterHandler::operator=(rhs);
3284        cbcModel_ = rhs.cbcModel_;
3285    }
3286    return *this;
3287}
3288//-------------------------------------------------------------------
3289// Clone
3290//-------------------------------------------------------------------
3291ClpDisasterHandler * CbcDisasterHandler::clone() const
3292{
3293    return new CbcDisasterHandler(*this);
3294}
3295// Type of disaster 0 can fix, 1 abort
3296int
3297CbcDisasterHandler::typeOfDisaster()
3298{
3299    if (!cbcModel_->parentModel() &&
3300            (cbcModel_->specialOptions()&2048) == 0) {
3301        return 0;
3302    } else {
3303        if (cbcModel_->parentModel())
3304            cbcModel_->setMaximumNodes(0);
3305        return 1;
3306    }
3307}
3308/* set model. */
3309void
3310CbcDisasterHandler::setCbcModel(CbcModel * model)
3311{
3312    cbcModel_ = model;
3313    if (model) {
3314        osiModel_
3315        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
3316        if (osiModel_)
3317            setSimplex(osiModel_->getModelPtr());
3318        else
3319            setSimplex(NULL);
3320    }
3321}
3322#endif
3323
Note: See TracBrowser for help on using the repository browser.