source: releases/2.7.1/Cbc/src/CbcHeuristicFPump.cpp @ 1995

Last change on this file since 1995 was 1573, checked in by lou, 8 years ago

Change to EPL license notice.

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