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

Last change on this file since 1514 was 1393, checked in by lou, 10 years ago

Mark #if 0 with JJF_ZERO and #if 1 with JJF_ONE as a historical reference
point.

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