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

Last change on this file since 1955 was 1955, checked in by forrest, 6 years ago

add more twiddles for advanced use of heuristics

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