source: stable/2.8/Cbc/src/CbcHeuristicFPump.cpp

Last change on this file was 1883, checked in by stefan, 7 years ago

sync with trunk rev 1882

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