source: trunk/Cbc/src/CbcHeuristic.cpp @ 1943

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

more options, copy statistics structure analysis
start coding of "switch" variables i.e. badly scaled ints or hi/lo
changes to allow more influence on small branch and bound
changes to get correct printout with threads

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 115.9 KB
Line 
1/* $Id: CbcHeuristic.cpp 1943 2013-07-21 09:05:45Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10
11#include "CbcConfig.h"
12
13#include <cassert>
14#include <cstdlib>
15#include <cmath>
16#include <cfloat>
17
18//#define PRINT_DEBUG
19#ifdef COIN_HAS_CLP
20#include "OsiClpSolverInterface.hpp"
21#endif
22#include "CbcModel.hpp"
23#include "CbcMessage.hpp"
24#include "CbcHeuristic.hpp"
25#include "CbcHeuristicFPump.hpp"
26#include "CbcEventHandler.hpp"
27#include "CbcStrategy.hpp"
28#include "CglPreProcess.hpp"
29#include "CglGomory.hpp"
30#include "CglProbing.hpp"
31#include "OsiAuxInfo.hpp"
32#include "OsiPresolve.hpp"
33#include "CbcBranchActual.hpp"
34#include "CbcCutGenerator.hpp"
35//==============================================================================
36
37CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode& rhs)
38{
39    numObjects_ = rhs.numObjects_;
40    brObj_ = new CbcBranchingObject*[numObjects_];
41    for (int i = 0; i < numObjects_; ++i) {
42        brObj_[i] = rhs.brObj_[i]->clone();
43    }
44}
45
46void
47CbcHeuristicNodeList::gutsOfDelete()
48{
49    for (int i = (static_cast<int>(nodes_.size())) - 1; i >= 0; --i) {
50        delete nodes_[i];
51    }
52}
53
54void
55CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList& rhs)
56{
57    append(rhs);
58}
59
60CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs)
61{
62    gutsOfCopy(rhs);
63}
64
65CbcHeuristicNodeList& CbcHeuristicNodeList::operator=
66(const CbcHeuristicNodeList & rhs)
67{
68    if (this != &rhs) {
69        gutsOfDelete();
70        gutsOfCopy(rhs);
71    }
72    return *this;
73}
74
75CbcHeuristicNodeList::~CbcHeuristicNodeList()
76{
77    gutsOfDelete();
78}
79
80void
81CbcHeuristicNodeList::append(CbcHeuristicNode*& node)
82{
83    nodes_.push_back(node);
84    node = NULL;
85}
86
87void
88CbcHeuristicNodeList::append(const CbcHeuristicNodeList& nodes)
89{
90    nodes_.reserve(nodes_.size() + nodes.size());
91    for (int i = 0; i < nodes.size(); ++i) {
92        CbcHeuristicNode* node = new CbcHeuristicNode(*nodes.node(i));
93        append(node);
94    }
95}
96
97//==============================================================================
98#define DEFAULT_WHERE ((255-2-16)*(1+256))
99// Default Constructor
100CbcHeuristic::CbcHeuristic() :
101        model_(NULL),
102        when_(2),
103        numberNodes_(200),
104        feasibilityPumpOptions_(-1),
105        fractionSmall_(1.0),
106        heuristicName_("Unknown"),
107        howOften_(1),
108        decayFactor_(0.0),
109        switches_(0),
110        whereFrom_(DEFAULT_WHERE),
111        shallowDepth_(1),
112        howOftenShallow_(1),
113        numInvocationsInShallow_(0),
114        numInvocationsInDeep_(0),
115        lastRunDeep_(0),
116        numRuns_(0),
117        minDistanceToRun_(1),
118        runNodes_(),
119        numCouldRun_(0),
120        numberSolutionsFound_(0),
121        numberNodesDone_(0),
122        inputSolution_(NULL)
123{
124    // As CbcHeuristic virtual need to modify .cpp if above change
125}
126
127// Constructor from model
128CbcHeuristic::CbcHeuristic(CbcModel & model) :
129        model_(&model),
130        when_(2),
131        numberNodes_(200),
132        feasibilityPumpOptions_(-1),
133        fractionSmall_(1.0),
134        heuristicName_("Unknown"),
135        howOften_(1),
136        decayFactor_(0.0),
137        switches_(0),
138        whereFrom_(DEFAULT_WHERE),
139        shallowDepth_(1),
140        howOftenShallow_(1),
141        numInvocationsInShallow_(0),
142        numInvocationsInDeep_(0),
143        lastRunDeep_(0),
144        numRuns_(0),
145        minDistanceToRun_(1),
146        runNodes_(),
147        numCouldRun_(0),
148        numberSolutionsFound_(0),
149        numberNodesDone_(0),
150        inputSolution_(NULL)
151{
152}
153
154void
155CbcHeuristic::gutsOfCopy(const CbcHeuristic & rhs)
156{
157    model_ = rhs.model_;
158    when_ = rhs.when_;
159    numberNodes_ = rhs.numberNodes_;
160    feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
161    fractionSmall_ = rhs.fractionSmall_;
162    randomNumberGenerator_ = rhs.randomNumberGenerator_;
163    heuristicName_ = rhs.heuristicName_;
164    howOften_ = rhs.howOften_;
165    decayFactor_ = rhs.decayFactor_;
166    switches_ = rhs.switches_;
167    whereFrom_ = rhs.whereFrom_;
168    shallowDepth_ = rhs.shallowDepth_;
169    howOftenShallow_ = rhs.howOftenShallow_;
170    numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
171    numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
172    lastRunDeep_ = rhs.lastRunDeep_;
173    numRuns_ = rhs.numRuns_;
174    numCouldRun_ = rhs.numCouldRun_;
175    minDistanceToRun_ = rhs.minDistanceToRun_;
176    runNodes_ = rhs.runNodes_;
177    numberSolutionsFound_ = rhs.numberSolutionsFound_;
178    numberNodesDone_ = rhs.numberNodesDone_;
179    if (rhs.inputSolution_) {
180        int numberColumns = model_->getNumCols();
181        setInputSolution(rhs.inputSolution_, rhs.inputSolution_[numberColumns]);
182    }
183}
184// Copy constructor
185CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
186{
187    inputSolution_ = NULL;
188    gutsOfCopy(rhs);
189}
190
191// Assignment operator
192CbcHeuristic &
193CbcHeuristic::operator=( const CbcHeuristic & rhs)
194{
195    if (this != &rhs) {
196        gutsOfDelete();
197        gutsOfCopy(rhs);
198    }
199    return *this;
200}
201
202void CbcHeurDebugNodes(CbcModel* model_)
203{
204    CbcNode* node = model_->currentNode();
205    CbcNodeInfo* nodeInfo = node->nodeInfo();
206    std::cout << "===============================================================\n";
207    while (nodeInfo) {
208        const CbcNode* node = nodeInfo->owner();
209        printf("nodeinfo: node %i\n", nodeInfo->nodeNumber());
210        {
211            const CbcIntegerBranchingObject* brPrint =
212                dynamic_cast<const CbcIntegerBranchingObject*>(nodeInfo->parentBranch());
213            if (!brPrint) {
214                printf("    parentBranch: NULL\n");
215            } else {
216                const double* downBounds = brPrint->downBounds();
217                const double* upBounds = brPrint->upBounds();
218                int variable = brPrint->variable();
219                int way = brPrint->way();
220                printf("   parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
221                       variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
222                       static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
223            }
224        }
225        if (! node) {
226            printf("    owner: NULL\n");
227        } else {
228            printf("    owner: node %i depth %i onTree %i active %i",
229                   node->nodeNumber(), node->depth(), node->onTree(), node->active());
230            const OsiBranchingObject* osibr =
231                nodeInfo->owner()->branchingObject();
232            const CbcBranchingObject* cbcbr =
233                dynamic_cast<const CbcBranchingObject*>(osibr);
234            const CbcIntegerBranchingObject* brPrint =
235                dynamic_cast<const CbcIntegerBranchingObject*>(cbcbr);
236            if (!brPrint) {
237                printf("        ownerBranch: NULL\n");
238            } else {
239                const double* downBounds = brPrint->downBounds();
240                const double* upBounds = brPrint->upBounds();
241                int variable = brPrint->variable();
242                int way = brPrint->way();
243                printf("        ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
244                       variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
245                       static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
246            }
247        }
248        nodeInfo = nodeInfo->parent();
249    }
250}
251
252void
253CbcHeuristic::debugNodes()
254{
255    CbcHeurDebugNodes(model_);
256}
257
258void
259CbcHeuristic::printDistanceToNodes()
260{
261    const CbcNode* currentNode = model_->currentNode();
262    if (currentNode != NULL) {
263        CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
264        for (int i = runNodes_.size() - 1; i >= 0; --i) {
265            nodeDesc->distance(runNodes_.node(i));
266        }
267        runNodes_.append(nodeDesc);
268    }
269}
270
271bool
272CbcHeuristic::shouldHeurRun(int whereFrom)
273{
274    assert (whereFrom >= 0 && whereFrom < 16);
275    // take off 8 (code - likes new solution)
276    whereFrom &= 7;
277    if ((whereFrom_&(1 << whereFrom)) == 0)
278        return false;
279    // No longer used for original purpose - so use for ever run at all JJF
280#ifndef JJF_ONE
281    // Don't run if hot start
282    if (model_ && model_->hotstartSolution())
283        return false;
284    else
285        return true;
286#else
287#ifdef JJF_ZERO
288    const CbcNode* currentNode = model_->currentNode();
289    if (currentNode == NULL) {
290        return false;
291    }
292
293    debugNodes();
294//   return false;
295
296    const int depth = currentNode->depth();
297#else
298    int depth = model_->currentDepth();
299#endif
300
301    const int nodeCount = model_->getNodeCount();  // FIXME: check that this is
302    // correct in parallel
303
304    if (nodeCount == 0 || depth <= shallowDepth_) {
305        // what to do when we are in the shallow part of the tree
306        if (model_->getCurrentPassNumber() == 1) {
307            // first time in the node...
308            numInvocationsInShallow_ = 0;
309        }
310        ++numInvocationsInShallow_;
311        // Very large howOftenShallow_ will give the original test:
312        // (model_->getCurrentPassNumber() != 1)
313        //    if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
314        if ((numInvocationsInShallow_ % howOftenShallow_) != 0) {
315            return false;
316        }
317        // LL: should we save these nodes in the list of nodes where the heur was
318        // LL: run?
319#ifndef JJF_ONE
320        if (currentNode != NULL) {
321            // Get where we are and create the appropriate CbcHeuristicNode object
322            CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
323            runNodes_.append(nodeDesc);
324        }
325#endif
326    } else {
327        // deeper in the tree
328        if (model_->getCurrentPassNumber() == 1) {
329            // first time in the node...
330            ++numInvocationsInDeep_;
331        }
332        if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
333            return false;
334        }
335        if (model_->getCurrentPassNumber() != 1) {
336            // Run the heuristic only when first entering the node.
337            // LL: I don't think this is right. It should run just before strong
338            // LL: branching, I believe.
339            return false;
340        }
341        // Get where we are and create the appropriate CbcHeuristicNode object
342        CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
343        //#ifdef PRINT_DEBUG
344#ifndef JJF_ONE
345        const double minDistanceToRun = 1.5 * log((double)depth) / log((double)2);
346#else
347    const double minDistanceToRun = minDistanceToRun_;
348#endif
349#ifdef PRINT_DEBUG
350        double minDistance = nodeDesc->minDistance(runNodes_);
351        std::cout << "minDistance = " << minDistance
352                  << ", minDistanceToRun = " << minDistanceToRun << std::endl;
353#endif
354        if (nodeDesc->minDistanceIsSmall(runNodes_, minDistanceToRun)) {
355            delete nodeDesc;
356            return false;
357        }
358        runNodes_.append(nodeDesc);
359        lastRunDeep_ = numInvocationsInDeep_;
360        //    ++lastRunDeep_;
361    }
362    ++numRuns_;
363    return true;
364#endif
365}
366
367bool
368CbcHeuristic::shouldHeurRun_randomChoice()
369{
370    if (!when_)
371        return false;
372    int depth = model_->currentDepth();
373    // when_ -999 is special marker to force to run
374    if (depth != 0 && when_ != -999) {
375        const double numerator = depth * depth;
376        const double denominator = exp(depth * log(2.0));
377        double probability = numerator / denominator;
378        double randomNumber = randomNumberGenerator_.randomDouble();
379        int when = when_ % 100;
380        if (when > 2 && when < 8) {
381            /* JJF adjustments
382            3 only at root and if no solution
383            4 only at root and if this heuristic has not got solution
384            5 as 3 but decay more
385            6 decay
386            7 run up to 2 times if solution found 4 otherwise
387            */
388            switch (when) {
389            case 3:
390            default:
391                if (model_->bestSolution())
392                    probability = -1.0;
393                break;
394            case 4:
395                if (numberSolutionsFound_)
396                    probability = -1.0;
397                break;
398            case 5:
399                assert (decayFactor_);
400                if (model_->bestSolution()) {
401                    probability = -1.0;
402                } else if (numCouldRun_ > 1000) {
403                    decayFactor_ *= 0.99;
404                    probability *= decayFactor_;
405                }
406                break;
407            case 6:
408                if (depth >= 3) {
409                    if ((numCouldRun_ % howOften_) == 0 &&
410                            numberSolutionsFound_*howOften_ < numCouldRun_) {
411                      //#define COIN_DEVELOP
412#ifdef COIN_DEVELOP
413                        int old = howOften_;
414#endif
415                        howOften_ = CoinMin(CoinMax(static_cast<int> (howOften_ * 1.1), howOften_ + 1), 1000000);
416#ifdef COIN_DEVELOP
417                        printf("Howoften changed from %d to %d for %s\n",
418                               old, howOften_, heuristicName_.c_str());
419#endif
420                    }
421                    probability = 1.0 / howOften_;
422                    if (model_->bestSolution())
423                        probability *= 0.5;
424                }
425                break;
426            case 7:
427                if ((model_->bestSolution() && numRuns_ >= 2) || numRuns_ >= 4)
428                    probability = -1.0;
429                break;
430            }
431        }
432        if (randomNumber > probability)
433            return false;
434
435        if (model_->getCurrentPassNumber() > 1)
436            return false;
437#ifdef COIN_DEVELOP
438        printf("Running %s, random %g probability %g\n",
439               heuristicName_.c_str(), randomNumber, probability);
440#endif
441    } else {
442#ifdef COIN_DEVELOP
443        printf("Running %s, depth %d when %d\n",
444               heuristicName_.c_str(), depth, when_);
445#endif
446    }
447    ++numRuns_;
448    return true;
449}
450
451// Resets stuff if model changes
452void
453CbcHeuristic::resetModel(CbcModel * model)
454{
455    model_ = model;
456}
457// Set seed
458void
459CbcHeuristic::setSeed(int value)
460{
461    if (value==0) {
462      double time = fabs(CoinGetTimeOfDay());
463      while (time>=COIN_INT_MAX)
464        time *= 0.5;
465      value = static_cast<int>(time);
466      char printArray[100];
467      sprintf(printArray, "using time of day seed was changed from %d to %d",
468              randomNumberGenerator_.getSeed(), value);
469      if (model_)
470        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
471          << printArray
472          << CoinMessageEol;
473    }
474    randomNumberGenerator_.setSeed(value);
475}
476// Get seed
477int
478CbcHeuristic::getSeed() const
479{
480  return randomNumberGenerator_.getSeed();
481}
482
483// Create C++ lines to get to current state
484void
485CbcHeuristic::generateCpp( FILE * fp, const char * heuristic)
486{
487    // hard coded as CbcHeuristic virtual
488    if (when_ != 2)
489        fprintf(fp, "3  %s.setWhen(%d);\n", heuristic, when_);
490    else
491        fprintf(fp, "4  %s.setWhen(%d);\n", heuristic, when_);
492    if (numberNodes_ != 200)
493        fprintf(fp, "3  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
494    else
495        fprintf(fp, "4  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
496    if (feasibilityPumpOptions_ != -1)
497        fprintf(fp, "3  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
498    else
499        fprintf(fp, "4  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
500    if (fractionSmall_ != 1.0)
501        fprintf(fp, "3  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
502    else
503        fprintf(fp, "4  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
504    if (heuristicName_ != "Unknown")
505        fprintf(fp, "3  %s.setHeuristicName(\"%s\");\n",
506                heuristic, heuristicName_.c_str()) ;
507    else
508        fprintf(fp, "4  %s.setHeuristicName(\"%s\");\n",
509                heuristic, heuristicName_.c_str()) ;
510    if (decayFactor_ != 0.0)
511        fprintf(fp, "3  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
512    else
513        fprintf(fp, "4  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
514    if (switches_ != 0)
515        fprintf(fp, "3  %s.setSwitches(%d);\n", heuristic, switches_);
516    else
517        fprintf(fp, "4  %s.setSwitches(%d);\n", heuristic, switches_);
518    if (whereFrom_ != DEFAULT_WHERE)
519        fprintf(fp, "3  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
520    else
521        fprintf(fp, "4  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
522    if (shallowDepth_ != 1)
523        fprintf(fp, "3  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
524    else
525        fprintf(fp, "4  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
526    if (howOftenShallow_ != 1)
527        fprintf(fp, "3  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
528    else
529        fprintf(fp, "4  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
530    if (minDistanceToRun_ != 1)
531        fprintf(fp, "3  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
532    else
533        fprintf(fp, "4  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
534}
535// Destructor
536CbcHeuristic::~CbcHeuristic ()
537{
538    delete [] inputSolution_;
539}
540
541// update model
542void CbcHeuristic::setModel(CbcModel * model)
543{
544    model_ = model;
545}
546/* Clone but ..
547   type 0 clone solver, 1 clone continuous solver
548   Add 2 to say without integer variables which are at low priority
549   Add 4 to say quite likely infeasible so give up easily.*/
550OsiSolverInterface *
551CbcHeuristic::cloneBut(int type)
552{
553    OsiSolverInterface * solver;
554    if ((type&1) == 0 || !model_->continuousSolver())
555        solver = model_->solver()->clone();
556    else
557        solver = model_->continuousSolver()->clone();
558#ifdef COIN_HAS_CLP
559    OsiClpSolverInterface * clpSolver
560    = dynamic_cast<OsiClpSolverInterface *> (solver);
561#endif
562    if ((type&2) != 0) {
563        int n = model_->numberObjects();
564        int priority = model_->continuousPriority();
565        if (priority < COIN_INT_MAX) {
566            for (int i = 0; i < n; i++) {
567                const OsiObject * obj = model_->object(i);
568                const CbcSimpleInteger * thisOne =
569                    dynamic_cast <const CbcSimpleInteger *> (obj);
570                if (thisOne) {
571                    int iColumn = thisOne->columnNumber();
572                    if (thisOne->priority() >= priority)
573                        solver->setContinuous(iColumn);
574                }
575            }
576        }
577#ifdef COIN_HAS_CLP
578        if (clpSolver) {
579            for (int i = 0; i < n; i++) {
580                const OsiObject * obj = model_->object(i);
581                const CbcSimpleInteger * thisOne =
582                    dynamic_cast <const CbcSimpleInteger *> (obj);
583                if (thisOne) {
584                    int iColumn = thisOne->columnNumber();
585                    if (clpSolver->isOptionalInteger(iColumn))
586                        clpSolver->setContinuous(iColumn);
587                }
588            }
589        }
590#endif
591    }
592#ifdef COIN_HAS_CLP
593    if ((type&4) != 0 && clpSolver) {
594        int options = clpSolver->getModelPtr()->moreSpecialOptions();
595        clpSolver->getModelPtr()->setMoreSpecialOptions(options | 64);
596    }
597#endif
598    return solver;
599}
600// Whether to exit at once on gap
601bool
602CbcHeuristic::exitNow(double bestObjective) const
603{
604    if ((switches_&2048) != 0) {
605        // exit may be forced - but unset for next time
606        switches_ &= ~2048;
607        if ((switches_&1024) != 0)
608            return true;
609    } else if ((switches_&1) == 0) {
610        return false;
611    }
612    // See if can stop on gap
613    OsiSolverInterface * solver = model_->solver();
614    double bestPossibleObjective = solver->getObjValue() * solver->getObjSense();
615    double absGap = CoinMax(model_->getAllowableGap(),
616                            model_->getHeuristicGap());
617    double fracGap = CoinMax(model_->getAllowableFractionGap(),
618                             model_->getHeuristicFractionGap());
619    double testGap = CoinMax(absGap, fracGap *
620                             CoinMax(fabs(bestObjective),
621                                     fabs(bestPossibleObjective)));
622
623    if (bestObjective - bestPossibleObjective < testGap
624            && model_->getCutoffIncrement() >= 0.0) {
625        return true;
626    } else {
627        return false;
628    }
629}
630#ifdef HISTORY_STATISTICS
631extern bool getHistoryStatistics_;
632#endif
633static double sizeRatio(int numberRowsNow, int numberColumnsNow,
634                        int numberRowsStart, int numberColumnsStart)
635{
636    double valueNow;
637    if (numberRowsNow*10 > numberColumnsNow || numberColumnsNow < 200) {
638        valueNow = 2 * numberRowsNow + numberColumnsNow;
639    } else {
640        // long and thin - rows are more important
641        if (numberRowsNow*40 > numberColumnsNow)
642            valueNow = 10 * numberRowsNow + numberColumnsNow;
643        else
644            valueNow = 200 * numberRowsNow + numberColumnsNow;
645    }
646    double valueStart;
647    if (numberRowsStart*10 > numberColumnsStart || numberColumnsStart < 200) {
648        valueStart = 2 * numberRowsStart + numberColumnsStart;
649    } else {
650        // long and thin - rows are more important
651        if (numberRowsStart*40 > numberColumnsStart)
652            valueStart = 10 * numberRowsStart + numberColumnsStart;
653        else
654            valueStart = 200 * numberRowsStart + numberColumnsStart;
655    }
656    //printf("sizeProblem Now %g, %d rows, %d columns\nsizeProblem Start %g, %d rows, %d columns\n",
657    // valueNow,numberRowsNow,numberColumnsNow,
658    // valueStart,numberRowsStart,numberColumnsStart);
659    if (10*numberRowsNow < 8*numberRowsStart || 10*numberColumnsNow < 7*numberColumnsStart)
660        return valueNow / valueStart;
661    else if (10*numberRowsNow < 9*numberRowsStart)
662        return 1.1*(valueNow / valueStart);
663    else if (numberRowsNow < numberRowsStart)
664        return 1.5*(valueNow / valueStart);
665    else
666        return 2.0*(valueNow / valueStart);
667}
668
669//static int saveModel=0;
670// Do mini branch and bound (return 1 if solution)
671int
672CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver, int numberNodes,
673                                  double * newSolution, double & newSolutionValue,
674                                  double cutoff, std::string name) const
675{
676  CbcEventHandler *eventHandler = model_->getEventHandler() ;
677  // Use this fraction
678  double fractionSmall = fractionSmall_;
679  int maximumSolutions =  model_->getMaximumSolutions();
680  int iterationMultiplier = 100;
681  if (eventHandler) {
682    typedef struct {
683      double fractionSmall;
684      double spareDouble[3];
685      OsiSolverInterface * solver;
686      void * sparePointer[2];
687      int numberNodes;
688      int maximumSolutions;
689      int iterationMultiplier;
690      int howOften;
691      int spareInt[3];
692    } SmallMod;
693    SmallMod temp;
694    temp.solver=solver;
695    temp.fractionSmall=fractionSmall;
696    temp.numberNodes=numberNodes;
697    temp.iterationMultiplier=iterationMultiplier;
698    temp.howOften=howOften_;
699    temp.maximumSolutions=maximumSolutions;
700    CbcEventHandler::CbcAction status = 
701      eventHandler->event(CbcEventHandler::smallBranchAndBound,
702                          &temp);
703    if (status==CbcEventHandler::killSolution)
704      return -1;
705    if (status==CbcEventHandler::takeAction) {
706      fractionSmall=temp.fractionSmall;
707      numberNodes=temp.numberNodes;
708      iterationMultiplier=temp.iterationMultiplier;
709      howOften_=temp.howOften;
710      maximumSolutions=temp.maximumSolutions;
711    }
712  }
713#if 0 
714  if (saveModel || model_->getMaximumSolutions()==100) {
715    printf("writing model\n");
716    solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
717  }
718#endif
719    // size before
720    int shiftRows = 0;
721    if (numberNodes < 0)
722        shiftRows = solver->getNumRows() - numberNodes_;
723    int numberRowsStart = solver->getNumRows() - shiftRows;
724    int numberColumnsStart = solver->getNumCols();
725#ifdef CLP_INVESTIGATE
726    printf("%s has %d rows, %d columns\n",
727           name.c_str(), solver->getNumRows(), solver->getNumCols());
728#endif
729    double before = 2 * numberRowsStart + numberColumnsStart;
730    if (before > 40000.0) {
731        // fairly large - be more conservative
732        double multiplier = 1.0 - 0.3 * CoinMin(100000.0, before - 40000.0) / 100000.0;
733        if (multiplier < 1.0) {
734            fractionSmall *= multiplier;
735#ifdef CLP_INVESTIGATE
736            printf("changing fractionSmall from %g to %g for %s\n",
737                   fractionSmall_, fractionSmall, name.c_str());
738#endif
739        }
740    }
741#ifdef COIN_HAS_CLP
742    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
743    if (osiclp && (osiclp->specialOptions()&65536) == 0) {
744        // go faster stripes
745        if (osiclp->getNumRows() < 300 && osiclp->getNumCols() < 500) {
746            osiclp->setupForRepeatedUse(2, 0);
747        } else {
748            osiclp->setupForRepeatedUse(0, 0);
749        }
750        // Turn this off if you get problems
751        // Used to be automatically set
752        osiclp->setSpecialOptions(osiclp->specialOptions() | (128 + 64 - 128));
753        ClpSimplex * lpSolver = osiclp->getModelPtr();
754        lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
755        lpSolver->setSpecialOptions(lpSolver->specialOptions() |
756                                    (/*16384+*/4096 + 512 + 128));
757    }
758#endif
759#ifdef HISTORY_STATISTICS
760    getHistoryStatistics_ = false;
761#endif
762#ifdef COIN_DEVELOP
763    int status = 0;
764#endif
765    int logLevel = model_->logLevel();
766#define LEN_PRINT 250
767    char generalPrint[LEN_PRINT];
768    // Do presolve to see if possible
769    int numberColumns = solver->getNumCols();
770    char * reset = NULL;
771    int returnCode = 1;
772    int saveModelOptions = model_->specialOptions();
773    //assert ((saveModelOptions&2048) == 0);
774    model_->setSpecialOptions(saveModelOptions | 2048);
775    if (fractionSmall<1.0) {
776        int saveLogLevel = solver->messageHandler()->logLevel();
777        if (saveLogLevel == 1) 
778            solver->messageHandler()->setLogLevel(0);
779        OsiPresolve * pinfo = new OsiPresolve();
780        int presolveActions = 0;
781        // Allow dual stuff on integers
782        presolveActions = 1;
783        // Do not allow all +1 to be tampered with
784        //if (allPlusOnes)
785        //presolveActions |= 2;
786        // allow transfer of costs
787        // presolveActions |= 4;
788        pinfo->setPresolveActions(presolveActions);
789        OsiSolverInterface * presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
790        delete pinfo;
791        // see if too big
792
793        if (presolvedModel) {
794            int afterRows = presolvedModel->getNumRows();
795            int afterCols = presolvedModel->getNumCols();
796            //#define COIN_DEVELOP
797#ifdef COIN_DEVELOP_z
798            if (numberNodes < 0) {
799                solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
800                presolvedModel->writeMpsNative("after1.mps", NULL, NULL, 2, 1);
801            }
802#endif
803            delete presolvedModel;
804            double ratio = sizeRatio(afterRows - shiftRows, afterCols,
805                                     numberRowsStart, numberColumnsStart);
806            double after = 2 * afterRows + afterCols;
807            if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
808                // Need code to try again to compress further using used
809                const int * used =  model_->usedInSolution();
810                int maxUsed = 0;
811                int iColumn;
812                const double * lower = solver->getColLower();
813                const double * upper = solver->getColUpper();
814                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
815                    if (upper[iColumn] > lower[iColumn]) {
816                        if (solver->isBinary(iColumn))
817                            maxUsed = CoinMax(maxUsed, used[iColumn]);
818                    }
819                }
820                if (maxUsed) {
821                    reset = new char [numberColumns];
822                    int nFix = 0;
823                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
824                        reset[iColumn] = 0;
825                        if (upper[iColumn] > lower[iColumn]) {
826                            if (solver->isBinary(iColumn) && used[iColumn] == maxUsed) {
827                                bool setValue = true;
828                                if (maxUsed == 1) {
829                                    double randomNumber = randomNumberGenerator_.randomDouble();
830                                    if (randomNumber > 0.3)
831                                        setValue = false;
832                                }
833                                if (setValue) {
834                                    reset[iColumn] = 1;
835                                    solver->setColLower(iColumn, 1.0);
836                                    nFix++;
837                                }
838                            }
839                        }
840                    }
841                    pinfo = new OsiPresolve();
842                    presolveActions = 0;
843                    // Allow dual stuff on integers
844                    presolveActions = 1;
845                    // Do not allow all +1 to be tampered with
846                    //if (allPlusOnes)
847                    //presolveActions |= 2;
848                    // allow transfer of costs
849                    // presolveActions |= 4;
850                    pinfo->setPresolveActions(presolveActions);
851                    presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
852                    delete pinfo;
853                    if (presolvedModel) {
854                        // see if too big
855                        int afterRows2 = presolvedModel->getNumRows();
856                        int afterCols2 = presolvedModel->getNumCols();
857                        delete presolvedModel;
858                        double ratio = sizeRatio(afterRows2 - shiftRows, afterCols2,
859                                                 numberRowsStart, numberColumnsStart);
860                        double after = 2 * afterRows2 + afterCols2;
861                        if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
862                            sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
863                                    solver->getNumRows(), solver->getNumCols(),
864                                    afterRows, afterCols, nFix, afterRows2, afterCols2);
865                            // If much too big - give up
866                            if (ratio > 0.75)
867                                returnCode = -1;
868                        } else {
869                            sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
870                                    solver->getNumRows(), solver->getNumCols(),
871                                    afterRows, afterCols, nFix, afterRows2, afterCols2);
872                        }
873                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
874                        << generalPrint
875                        << CoinMessageEol;
876                    } else {
877                        returnCode = 2; // infeasible
878                    }
879                }
880            } else if (ratio > fractionSmall && after > 300 && numberNodes >=0) {
881                returnCode = -1;
882            }
883        } else {
884            returnCode = 2; // infeasible
885        }
886        solver->messageHandler()->setLogLevel(saveLogLevel);
887    }
888    if (returnCode == 2 || returnCode == -1) {
889        model_->setSpecialOptions(saveModelOptions);
890        delete [] reset;
891#ifdef HISTORY_STATISTICS
892        getHistoryStatistics_ = true;
893#endif
894        //printf("small no good\n");
895        return returnCode;
896    }
897    // Reduce printout
898    bool takeHint;
899    OsiHintStrength strength;
900    solver->getHintParam(OsiDoReducePrint, takeHint, strength);
901    solver->setHintParam(OsiDoReducePrint, true, OsiHintTry);
902    solver->setHintParam(OsiDoPresolveInInitial, false, OsiHintTry);
903    double signedCutoff = cutoff*solver->getObjSense();
904    solver->setDblParam(OsiDualObjectiveLimit, signedCutoff);
905    solver->initialSolve();
906    if (solver->isProvenOptimal()) {
907        CglPreProcess process;
908        OsiSolverInterface * solver2 = NULL;
909        if ((model_->moreSpecialOptions()&65536)!=0)
910          process.setOptions(2+4+8+16); // no cuts
911        else
912          process.setOptions(16); // no complicated dupcol stuff
913        /* Do not try and produce equality cliques and
914           do up to 2 passes (normally) 5 if restart */
915        int numberPasses = 2;
916        if ((model_->moreSpecialOptions2()&16)!=0) {
917          // quick
918          process.setOptions(2+4+8+16); // no cuts
919          numberPasses = 1;
920        }
921        if (numberNodes < 0) {
922          numberPasses = 5;
923          // Say some rows cuts
924          int numberRows = solver->getNumRows();
925          if (numberNodes_ < numberRows && true /* think */) {
926            char * type = new char[numberRows];
927            memset(type, 0, numberNodes_);
928            memset(type + numberNodes_, 1, numberRows - numberNodes_);
929            process.passInRowTypes(type, numberRows);
930            delete [] type;
931          }
932        }
933        if (logLevel <= 1)
934          process.messageHandler()->setLogLevel(0);
935        if (!solver->defaultHandler()&&
936            solver->messageHandler()->logLevel(0)!=-1000)
937          process.passInMessageHandler(solver->messageHandler());
938        solver2 = process.preProcessNonDefault(*solver, false,
939                                               numberPasses);
940          if (!solver2) {
941            if (logLevel > 1)
942              printf("Pre-processing says infeasible\n");
943            returnCode = 2; // so will be infeasible
944          } else {
945#ifdef COIN_DEVELOP_z
946            if (numberNodes < 0) {
947              solver2->writeMpsNative("after2.mps", NULL, NULL, 2, 1);
948            }
949#endif
950            // see if too big
951            double ratio = sizeRatio(solver2->getNumRows() - shiftRows, solver2->getNumCols(),
952                                     numberRowsStart, numberColumnsStart);
953            double after = 2 * solver2->getNumRows() + solver2->getNumCols();
954            if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
955                sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
956                        solver->getNumRows(), solver->getNumCols(),
957                        solver2->getNumRows(), solver2->getNumCols());
958                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
959                << generalPrint
960                << CoinMessageEol;
961                returnCode = -1;
962                //printf("small no good2\n");
963            } else {
964                sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns",
965                        solver->getNumRows(), solver->getNumCols(),
966                        solver2->getNumRows(), solver2->getNumCols());
967                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
968                << generalPrint
969                << CoinMessageEol;
970            }
971            if (returnCode == 1) {
972                solver2->resolve();
973                CbcModel model(*solver2);
974                // move seed across
975                model.randomNumberGenerator()->setSeed(model_->randomNumberGenerator()->getSeed());
976                if (numberNodes >= 0) {
977                    // normal
978                    model.setSpecialOptions(saveModelOptions | 2048);
979                    if (logLevel <= 1 && feasibilityPumpOptions_ != -3)
980                        model.setLogLevel(0);
981                    else
982                        model.setLogLevel(logLevel);
983                    // No small fathoming
984                    model.setFastNodeDepth(-1);
985                    model.setCutoff(signedCutoff);
986                    model.setStrongStrategy(0);
987                    // Don't do if original fraction > 1.0 and too large
988                    if (fractionSmall_>1.0 && fractionSmall_ < 1000000.0) {
989                      /* 1.4 means -1 nodes if >.4
990                         2.4 means -1 nodes if >.5 and 0 otherwise
991                         3.4 means -1 nodes if >.6 and 0 or 5
992                         4.4 means -1 nodes if >.7 and 0, 5 or 10
993                      */
994                      double fraction = fractionSmall_-floor(fractionSmall_);
995                      if (ratio>fraction) {
996                        int type = static_cast<int>(floor(fractionSmall_*0.1));
997                        int over = static_cast<int>(ceil(ratio-fraction));
998                        int maxNodes[]={-1,0,5,10};
999                        if (type>over)
1000                          numberNodes=maxNodes[type-over];
1001                        else
1002                          numberNodes=-1;
1003                      }
1004                    }
1005                    model.setMaximumNodes(numberNodes);
1006                    model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
1007                    if ((saveModelOptions&2048) == 0)
1008                      model.setMoreSpecialOptions(model_->moreSpecialOptions());
1009                      model.setMoreSpecialOptions2(model_->moreSpecialOptions2());
1010                    // off conflict analysis
1011                    model.setMoreSpecialOptions(model.moreSpecialOptions()&~4194304);
1012                   
1013                    // Lightweight
1014                    CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
1015                    model.setStrategy(strategy);
1016                    model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10);
1017                    model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot())));
1018                    model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses()));
1019                } else {
1020                    model.setSpecialOptions(saveModelOptions);
1021                    model_->messageHandler()->message(CBC_RESTART, model_->messages())
1022                    << solver2->getNumRows() << solver2->getNumCols()
1023                    << CoinMessageEol;
1024                    // going for full search and copy across more stuff
1025                    model.gutsOfCopy(*model_, 2);
1026                    assert (!model_->heuristicModel());
1027                    model_->setHeuristicModel(&model);
1028                    for (int i = 0; i < model.numberCutGenerators(); i++) {
1029                        CbcCutGenerator * generator = model.cutGenerator(i);
1030                        CglGomory * gomory = dynamic_cast<CglGomory *>
1031                          (generator->generator());
1032                        if (gomory&&gomory->originalSolver()) 
1033                          gomory->passInOriginalSolver(model.solver());
1034                        generator->setTiming(true);
1035                        // Turn on if was turned on
1036                        int iOften = model_->cutGenerator(i)->howOften();
1037#ifdef CLP_INVESTIGATE
1038                        printf("Gen %d often %d %d\n",
1039                               i, generator->howOften(),
1040                               iOften);
1041#endif
1042                        if (iOften > 0)
1043                            generator->setHowOften(iOften % 1000000);
1044                        if (model_->cutGenerator(i)->howOftenInSub() == -200)
1045                            generator->setHowOften(-100);
1046                    }
1047                    model.setCutoff(signedCutoff);
1048                    // make sure can't do nested search! but allow heuristics
1049                    model.setSpecialOptions((model.specialOptions()&(~(512 + 2048))) | 1024);
1050                    bool takeHint;
1051                    OsiHintStrength strength;
1052                    // Switch off printing if asked to
1053                    model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength);
1054                    model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength);
1055                    // no cut generators if none in parent
1056                    CbcStrategyDefault
1057                      strategy(model_->numberCutGenerators() ? 1 : -1, 
1058                               model_->numberStrong(),
1059                               model_->numberBeforeTrust());
1060                    // Set up pre-processing - no
1061                    strategy.setupPreProcessing(0); // was (4);
1062                    model.setStrategy(strategy);
1063                    //model.solver()->writeMps("crunched");
1064                    int numberCuts = process.cuts().sizeRowCuts();
1065                    if (numberCuts) {
1066                        // add in cuts
1067                        CglStored cuts = process.cuts();
1068                        model.addCutGenerator(&cuts, 1, "Stored from first");
1069                        model.cutGenerator(model.numberCutGenerators()-1)->setGlobalCuts(true);
1070                    }
1071                }
1072                // Do search
1073                if (logLevel > 1)
1074                    model_->messageHandler()->message(CBC_START_SUB, model_->messages())
1075                    << name
1076                    << model.getMaximumNodes()
1077                    << CoinMessageEol;
1078                // probably faster to use a basis to get integer solutions
1079                model.setSpecialOptions(model.specialOptions() | 2);
1080#ifdef CBC_THREAD
1081                if (model_->getNumberThreads() > 0 && (model_->getThreadMode()&4) != 0) {
1082                    // See if at root node
1083                    bool atRoot = model_->getNodeCount() == 0;
1084                    int passNumber = model_->getCurrentPassNumber();
1085                    if (atRoot && passNumber == 1)
1086                        model.setNumberThreads(model_->getNumberThreads());
1087                }
1088#endif
1089                model.setParentModel(*model_);
1090                model.setMaximumSolutions(maximumSolutions); 
1091                model.setOriginalColumns(process.originalColumns());
1092                model.setSearchStrategy(-1);
1093                // If no feasibility pump then insert a lightweight one
1094                if (feasibilityPumpOptions_ >= 0 || feasibilityPumpOptions_ == -2) {
1095                    CbcHeuristicFPump * fpump = NULL;
1096                    for (int i = 0; i < model.numberHeuristics(); i++) {
1097                        CbcHeuristicFPump* pump =
1098                            dynamic_cast<CbcHeuristicFPump*>(model.heuristic(i));
1099                        if (pump) {
1100                            fpump = pump;
1101                            break;
1102                        }
1103                    }
1104                    if (!fpump) {
1105                        CbcHeuristicFPump heuristic4;
1106                        // use any cutoff
1107                        heuristic4.setFakeCutoff(0.5*COIN_DBL_MAX);
1108                        if (fractionSmall_<=1.0) 
1109                          heuristic4.setMaximumPasses(10);
1110                        int pumpTune = feasibilityPumpOptions_;
1111                        if (pumpTune==-2)
1112                          pumpTune = 4; // proximity
1113                        if (pumpTune > 0) {
1114                            /*
1115                            >=10000000 for using obj
1116                            >=1000000 use as accumulate switch
1117                            >=1000 use index+1 as number of large loops
1118                            >=100 use 0.05 objvalue as increment
1119                            %100 == 10,20 etc for experimentation
1120                            1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
1121                            4 and static continuous, 5 as 3 but no internal integers
1122                            6 as 3 but all slack basis!
1123                            */
1124                            double value = solver2->getObjSense() * solver2->getObjValue();
1125                            int w = pumpTune / 10;
1126                            int ix = w % 10;
1127                            w /= 10;
1128                            int c = w % 10;
1129                            w /= 10;
1130                            int r = w;
1131                            int accumulate = r / 1000;
1132                            r -= 1000 * accumulate;
1133                            if (accumulate >= 10) {
1134                                int which = accumulate / 10;
1135                                accumulate -= 10 * which;
1136                                which--;
1137                                // weights and factors
1138                                double weight[] = {0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};
1139                                double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};
1140                                heuristic4.setInitialWeight(weight[which]);
1141                                heuristic4.setWeightFactor(factor[which]);
1142                            }
1143                            // fake cutoff
1144                            if (c) {
1145                                double cutoff;
1146                                solver2->getDblParam(OsiDualObjectiveLimit, cutoff);
1147                                cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c);
1148                                heuristic4.setFakeCutoff(cutoff);
1149                            }
1150                            if (r) {
1151                                // also set increment
1152                                //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
1153                                double increment = 0.0;
1154                                heuristic4.setAbsoluteIncrement(increment);
1155                                heuristic4.setAccumulate(accumulate);
1156                                heuristic4.setMaximumRetries(r + 1);
1157                            }
1158                            pumpTune = pumpTune % 100;
1159                            if (pumpTune == 6)
1160                                pumpTune = 13;
1161                            if (pumpTune != 13)
1162                                pumpTune = pumpTune % 10;
1163                            heuristic4.setWhen(pumpTune);
1164                            if (ix) {
1165                                heuristic4.setFeasibilityPumpOptions(ix*10);
1166                            }
1167                        }
1168                        model.addHeuristic(&heuristic4, "feasibility pump", 0);
1169       
1170                    }
1171                } else if (feasibilityPumpOptions_==-3) {
1172                  // add all (except this)
1173                  for (int i = 0; i < model_->numberHeuristics(); i++) {
1174                    if (strcmp(heuristicName(),model_->heuristic(i)->heuristicName()))
1175                      model.addHeuristic(model_->heuristic(i)); 
1176                  }
1177                }
1178                //printf("sol %x\n",inputSolution_);
1179                if (inputSolution_) {
1180                    // translate and add a serendipity heuristic
1181                    int numberColumns = solver2->getNumCols();
1182                    const int * which = process.originalColumns();
1183                    OsiSolverInterface * solver3 = solver2->clone();
1184                    for (int i = 0; i < numberColumns; i++) {
1185                        if (solver3->isInteger(i)) {
1186                            int k = which[i];
1187                            double value = inputSolution_[k];
1188                            //if (value)
1189                            //printf("orig col %d now %d val %g\n",
1190                            //       k,i,value);
1191                            solver3->setColLower(i, value);
1192                            solver3->setColUpper(i, value);
1193                        }
1194                    }
1195                    solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1196                    solver3->resolve();
1197                    if (!solver3->isProvenOptimal()) {
1198                        // Try just setting nonzeros
1199                        OsiSolverInterface * solver4 = solver2->clone();
1200                        for (int i = 0; i < numberColumns; i++) {
1201                            if (solver4->isInteger(i)) {
1202                                int k = which[i];
1203                                double value = floor(inputSolution_[k] + 0.5);
1204                                if (value) {
1205                                    solver3->setColLower(i, value);
1206                                    solver3->setColUpper(i, value);
1207                                }
1208                            }
1209                        }
1210                        solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1211                        solver4->resolve();
1212                        int nBad = -1;
1213                        if (solver4->isProvenOptimal()) {
1214                            nBad = 0;
1215                            const double * solution = solver4->getColSolution();
1216                            for (int i = 0; i < numberColumns; i++) {
1217                                if (solver4->isInteger(i)) {
1218                                    double value = floor(solution[i] + 0.5);
1219                                    if (fabs(value - solution[i]) > 1.0e-6)
1220                                        nBad++;
1221                                }
1222                            }
1223                        }
1224                        if (nBad) {
1225                            delete solver4;
1226                        } else {
1227                            delete solver3;
1228                            solver3 = solver4;
1229                        }
1230                    }
1231                    if (solver3->isProvenOptimal()) {
1232                        // good
1233                        CbcSerendipity heuristic(model);
1234                        double value = solver3->getObjSense() * solver3->getObjValue();
1235                        heuristic.setInputSolution(solver3->getColSolution(), value);
1236                        value = value + 1.0e-7*(1.0 + fabs(value));
1237                        value *= solver3->getObjSense();
1238                        model.setCutoff(value);
1239                        model.addHeuristic(&heuristic, "Previous solution", 0);
1240                        //printf("added seren\n");
1241                    } else {
1242                        double value = model_->getMinimizationObjValue();
1243                        value = value + 1.0e-7*(1.0 + fabs(value));
1244                        value *= solver3->getObjSense();
1245                        model.setCutoff(value);
1246#ifdef CLP_INVESTIGATE
1247                        printf("NOT added seren\n");
1248                        solver3->writeMps("bad_seren");
1249                        solver->writeMps("orig_seren");
1250#endif
1251                    }
1252                    delete solver3;
1253                }
1254                if (model_->searchStrategy() == 2) {
1255                    model.setNumberStrong(5);
1256                    model.setNumberBeforeTrust(5);
1257                }
1258                if (model.getNumCols()) {
1259                    if (numberNodes >= 0) {
1260                        setCutAndHeuristicOptions(model);
1261                        // not too many iterations
1262                        model.setMaximumNumberIterations(iterationMultiplier*(numberNodes + 10));
1263                        // Not fast stuff
1264                        model.setFastNodeDepth(-1);
1265                    } else if (model.fastNodeDepth() >= 1000000) {
1266                        // already set
1267                        model.setFastNodeDepth(model.fastNodeDepth() - 1000000);
1268                    }
1269                    model.setWhenCuts(999998);
1270#define ALWAYS_DUAL
1271#ifdef ALWAYS_DUAL
1272                    OsiSolverInterface * solverD = model.solver();
1273                    bool takeHint;
1274                    OsiHintStrength strength;
1275                    solverD->getHintParam(OsiDoDualInResolve, takeHint, strength);
1276                    solverD->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
1277#endif
1278                    model.passInEventHandler(model_->getEventHandler());
1279                    // say model_ is sitting there
1280                    int saveOptions = model_->specialOptions();
1281                    model_->setSpecialOptions(saveOptions|1048576);
1282                    // and switch off debugger
1283                    model.setSpecialOptions(model.specialOptions()&(~1));
1284#if 0 //def COIN_HAS_CLP
1285                    OsiClpSolverInterface * clpSolver
1286                      = dynamic_cast<OsiClpSolverInterface *> (model.solver());
1287                    if (clpSolver)
1288                      clpSolver->zapDebugger();
1289#endif
1290                    model.branchAndBound();
1291                    model_->setHeuristicModel(NULL);
1292                    model_->setSpecialOptions(saveOptions);
1293#ifdef ALWAYS_DUAL
1294                    solverD = model.solver();
1295                    solverD->setHintParam(OsiDoDualInResolve, takeHint, strength);
1296#endif
1297                    numberNodesDone_ = model.getNodeCount();
1298#ifdef COIN_DEVELOP
1299                    printf("sub branch %d nodes, %d iterations - max %d\n",
1300                           model.getNodeCount(), model.getIterationCount(),
1301                           100*(numberNodes + 10));
1302#endif
1303                    if (numberNodes < 0) {
1304                        model_->incrementIterationCount(model.getIterationCount());
1305                        model_->incrementNodeCount(model.getNodeCount());
1306                        // update best solution (in case ctrl-c)
1307                        // !!! not a good idea - think a bit harder
1308                        //model_->setMinimizationObjValue(model.getMinimizationObjValue());
1309                        for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) {
1310                            CbcCutGenerator * generator = model.cutGenerator(iGenerator);
1311                            sprintf(generalPrint,
1312                                    "%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)",
1313                                    generator->cutGeneratorName(),
1314                                    generator->numberTimesEntered(),
1315                                    generator->numberCutsInTotal() +
1316                                    generator->numberColumnCuts(),
1317                                    generator->numberCutsActive(),
1318                                    generator->timeInCutGenerator());
1319                            CglStored * stored = dynamic_cast<CglStored*>(generator->generator());
1320                            if (stored && !generator->numberCutsInTotal())
1321                                continue;
1322#ifndef CLP_INVESTIGATE
1323                            CglImplication * implication = dynamic_cast<CglImplication*>(generator->generator());
1324                            if (implication)
1325                                continue;
1326#endif
1327                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1328                            << generalPrint
1329                            << CoinMessageEol;
1330                        }
1331                    }
1332                } else {
1333                    // empty model
1334                    model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
1335                }
1336                if (logLevel > 1)
1337                    model_->messageHandler()->message(CBC_END_SUB, model_->messages())
1338                    << name
1339                    << CoinMessageEol;
1340                if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) {
1341                    // solution
1342                    if (model.getNumCols())
1343                        returnCode = model.isProvenOptimal() ? 3 : 1;
1344                    else
1345                        returnCode = 3;
1346                    // post process
1347#ifdef COIN_HAS_CLP
1348                    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
1349                    if (clpSolver) {
1350                        ClpSimplex * lpSolver = clpSolver->getModelPtr();
1351                        lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
1352                    }
1353#endif
1354                    //if (fractionSmall_ < 1000000.0)
1355                      process.postProcess(*model.solver());
1356                    if (solver->isProvenOptimal() && solver->getObjValue()*solver->getObjSense() < cutoff) {
1357                        // Solution now back in solver
1358                        int numberColumns = solver->getNumCols();
1359                        memcpy(newSolution, solver->getColSolution(),
1360                               numberColumns*sizeof(double));
1361                        newSolutionValue = model.getMinimizationObjValue();
1362                    } else {
1363                        // odd - but no good
1364                        returnCode = 0; // so will be infeasible
1365                    }
1366                } else {
1367                    // no good
1368                    returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
1369                }
1370                int totalNumberIterations = model.getIterationCount() +
1371                                            process.numberIterationsPre() +
1372                                            process.numberIterationsPost();
1373                if (totalNumberIterations > 100*(numberNodes + 10)
1374                    && fractionSmall_ < 1000000.0) {
1375                    // only allow smaller problems
1376                    fractionSmall = fractionSmall_;
1377                    fractionSmall_ *= 0.9;
1378#ifdef CLP_INVESTIGATE
1379                    printf("changing fractionSmall from %g to %g for %s as %d iterations\n",
1380                           fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations);
1381#endif
1382                }
1383                if (model.status() == 5)
1384                   model_->sayEventHappened();
1385#ifdef COIN_DEVELOP
1386                if (model.isProvenInfeasible())
1387                    status = 1;
1388                else if (model.isProvenOptimal())
1389                    status = 2;
1390#endif
1391            }
1392        }
1393    } else {
1394        returnCode = 2; // infeasible finished
1395        printf("Infeasible on initial solve\n");
1396    }
1397    model_->setSpecialOptions(saveModelOptions);
1398    model_->setLogLevel(logLevel);
1399    if (returnCode == 1 || returnCode == 2) {
1400        OsiSolverInterface * solverC = model_->continuousSolver();
1401        if (false && solverC) {
1402            const double * lower = solver->getColLower();
1403            const double * upper = solver->getColUpper();
1404            const double * lowerC = solverC->getColLower();
1405            const double * upperC = solverC->getColUpper();
1406            bool good = true;
1407            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1408                if (solverC->isInteger(iColumn)) {
1409                    if (lower[iColumn] > lowerC[iColumn] &&
1410                            upper[iColumn] < upperC[iColumn]) {
1411                        good = false;
1412                        printf("CUT - can't add\n");
1413                        break;
1414                    }
1415                }
1416            }
1417            if (good) {
1418                double * cut = new double [numberColumns];
1419                int * which = new int [numberColumns];
1420                double rhs = -1.0;
1421                int n = 0;
1422                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1423                    if (solverC->isInteger(iColumn)) {
1424                        if (lower[iColumn] == upperC[iColumn]) {
1425                            rhs += lower[iColumn];
1426                            cut[n] = 1.0;
1427                            which[n++] = iColumn;
1428                        } else if (upper[iColumn] == lowerC[iColumn]) {
1429                            rhs -= upper[iColumn];
1430                            cut[n] = -1.0;
1431                            which[n++] = iColumn;
1432                        }
1433                    }
1434                }
1435                printf("CUT has %d entries\n", n);
1436                OsiRowCut newCut;
1437                newCut.setLb(-COIN_DBL_MAX);
1438                newCut.setUb(rhs);
1439                newCut.setRow(n, which, cut, false);
1440                model_->makeGlobalCut(newCut);
1441                delete [] cut;
1442                delete [] which;
1443            }
1444        }
1445#ifdef COIN_DEVELOP
1446        if (status == 1)
1447            printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str());
1448        else if (status == 2)
1449            printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str());
1450#endif
1451    }
1452    if (reset) {
1453        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1454            if (reset[iColumn])
1455                solver->setColLower(iColumn, 0.0);
1456        }
1457        delete [] reset;
1458    }
1459#ifdef HISTORY_STATISTICS
1460    getHistoryStatistics_ = true;
1461#endif
1462    solver->setHintParam(OsiDoReducePrint, takeHint, strength);
1463    return returnCode;
1464}
1465// Set input solution
1466void
1467CbcHeuristic::setInputSolution(const double * solution, double objValue)
1468{
1469    delete [] inputSolution_;
1470    inputSolution_ = NULL;
1471    if (model_ && solution) {
1472        int numberColumns = model_->getNumCols();
1473        inputSolution_ = new double [numberColumns+1];
1474        memcpy(inputSolution_, solution, numberColumns*sizeof(double));
1475        inputSolution_[numberColumns] = objValue;
1476    }
1477}
1478
1479//##############################################################################
1480
1481inline int compare3BranchingObjects(const CbcBranchingObject* br0,
1482                                    const CbcBranchingObject* br1)
1483{
1484    const int t0 = br0->type();
1485    const int t1 = br1->type();
1486    if (t0 < t1) {
1487        return -1;
1488    }
1489    if (t0 > t1) {
1490        return 1;
1491    }
1492    return br0->compareOriginalObject(br1);
1493}
1494
1495//==============================================================================
1496
1497inline bool compareBranchingObjects(const CbcBranchingObject* br0,
1498                                    const CbcBranchingObject* br1)
1499{
1500    return compare3BranchingObjects(br0, br1) < 0;
1501}
1502
1503//==============================================================================
1504
1505void
1506CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
1507{
1508    //  CbcHeurDebugNodes(&model);
1509    CbcNode* node = model.currentNode();
1510    brObj_ = new CbcBranchingObject*[node->depth()];
1511    CbcNodeInfo* nodeInfo = node->nodeInfo();
1512    int cnt = 0;
1513    while (nodeInfo->parentBranch() != NULL) {
1514        const OsiBranchingObject* br = nodeInfo->parentBranch();
1515        const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br);
1516        if (! cbcbr) {
1517            throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
1518                            "gutsOfConstructor",
1519                            "CbcHeuristicNode",
1520                            __FILE__, __LINE__);
1521        }
1522        brObj_[cnt] = cbcbr->clone();
1523        brObj_[cnt]->previousBranch();
1524        ++cnt;
1525        nodeInfo = nodeInfo->parent();
1526    }
1527    std::sort(brObj_, brObj_ + cnt, compareBranchingObjects);
1528    if (cnt <= 1) {
1529        numObjects_ = cnt;
1530    } else {
1531        numObjects_ = 0;
1532        CbcBranchingObject* br = NULL; // What should this be?
1533        for (int i = 1; i < cnt; ++i) {
1534            if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
1535                int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0);
1536                switch (comp) {
1537                case CbcRangeSame: // the same range
1538                case CbcRangeDisjoint: // disjoint decisions
1539                    // should not happen! we are on a chain!
1540                    abort();
1541                case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
1542                    delete brObj_[i];
1543                    break;
1544                case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
1545                    delete brObj_[numObjects_];
1546                    brObj_[numObjects_] = brObj_[i];
1547                    break;
1548                case CbcRangeOverlap: // overlap
1549                    delete brObj_[i];
1550                    delete brObj_[numObjects_];
1551                    brObj_[numObjects_] = br;
1552                    break;
1553                }
1554                continue;
1555            } else {
1556                brObj_[++numObjects_] = brObj_[i];
1557            }
1558        }
1559        ++numObjects_;
1560    }
1561}
1562
1563//==============================================================================
1564
1565CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
1566{
1567    gutsOfConstructor(model);
1568}
1569
1570//==============================================================================
1571
1572double
1573CbcHeuristicNode::distance(const CbcHeuristicNode* node) const
1574{
1575
1576    const double disjointWeight = 1;
1577    const double overlapWeight = 0.4;
1578    const double subsetWeight = 0.2;
1579    int countDisjointWeight = 0;
1580    int countOverlapWeight = 0;
1581    int countSubsetWeight = 0;
1582    int i = 0;
1583    int j = 0;
1584    double dist = 0.0;
1585#ifdef PRINT_DEBUG
1586    printf(" numObjects_ = %i, node->numObjects_ = %i\n",
1587           numObjects_, node->numObjects_);
1588#endif
1589    while ( i < numObjects_ && j < node->numObjects_) {
1590        CbcBranchingObject* br0 = brObj_[i];
1591        const CbcBranchingObject* br1 = node->brObj_[j];
1592#ifdef PRINT_DEBUG
1593        const CbcIntegerBranchingObject* brPrint0 =
1594            dynamic_cast<const CbcIntegerBranchingObject*>(br0);
1595        const double* downBounds = brPrint0->downBounds();
1596        const double* upBounds = brPrint0->upBounds();
1597        int variable = brPrint0->variable();
1598        int way = brPrint0->way();
1599        printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1600               variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
1601               static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
1602        const CbcIntegerBranchingObject* brPrint1 =
1603            dynamic_cast<const CbcIntegerBranchingObject*>(br1);
1604        downBounds = brPrint1->downBounds();
1605        upBounds = brPrint1->upBounds();
1606        variable = brPrint1->variable();
1607        way = brPrint1->way();
1608        printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1609               variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
1610               static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
1611#endif
1612        const int brComp = compare3BranchingObjects(br0, br1);
1613        if (brComp < 0) {
1614            dist += subsetWeight;
1615            countSubsetWeight++;
1616            ++i;
1617        } else if (brComp > 0) {
1618            dist += subsetWeight;
1619            countSubsetWeight++;
1620            ++j;
1621        } else {
1622            const int comp = br0->compareBranchingObject(br1, false);
1623            switch (comp) {
1624            case CbcRangeSame:
1625                // do nothing
1626                break;
1627            case CbcRangeDisjoint: // disjoint decisions
1628                dist += disjointWeight;
1629                countDisjointWeight++;
1630                break;
1631            case CbcRangeSubset: // subset one way or another
1632            case CbcRangeSuperset:
1633                dist += subsetWeight;
1634                countSubsetWeight++;
1635                break;
1636            case CbcRangeOverlap: // overlap
1637                dist += overlapWeight;
1638                countOverlapWeight++;
1639                break;
1640            }
1641            ++i;
1642            ++j;
1643        }
1644    }
1645    dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
1646    countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
1647    COIN_DETAIL_PRINT(printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
1648                             countOverlapWeight, countDisjointWeight));
1649    return dist;
1650}
1651
1652//==============================================================================
1653
1654CbcHeuristicNode::~CbcHeuristicNode()
1655{
1656    for (int i = 0; i < numObjects_; ++i) {
1657        delete brObj_[i];
1658    }
1659    delete [] brObj_;
1660}
1661
1662//==============================================================================
1663
1664double
1665CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList) const
1666{
1667    double minDist = COIN_DBL_MAX;
1668    for (int i = nodeList.size() - 1; i >= 0; --i) {
1669        minDist = CoinMin(minDist, distance(nodeList.node(i)));
1670    }
1671    return minDist;
1672}
1673
1674//==============================================================================
1675
1676bool
1677CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList& nodeList,
1678                                     const double threshold) const
1679{
1680    for (int i = nodeList.size() - 1; i >= 0; --i) {
1681        if (distance(nodeList.node(i)) >= threshold) {
1682            continue;
1683        } else {
1684            return true;
1685        }
1686    }
1687    return false;
1688}
1689
1690//==============================================================================
1691
1692double
1693CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList) const
1694{
1695    if (nodeList.size() == 0) {
1696        return COIN_DBL_MAX;
1697    }
1698    double sumDist = 0;
1699    for (int i = nodeList.size() - 1; i >= 0; --i) {
1700        sumDist += distance(nodeList.node(i));
1701    }
1702    return sumDist / nodeList.size();
1703}
1704
1705//##############################################################################
1706
1707// Default Constructor
1708CbcRounding::CbcRounding()
1709        : CbcHeuristic()
1710{
1711    // matrix and row copy will automatically be empty
1712    seed_ = 7654321;
1713    down_ = NULL;
1714    up_ = NULL;
1715    equal_ = NULL;
1716    //whereFrom_ |= 16; // allow more often
1717}
1718
1719// Constructor from model
1720CbcRounding::CbcRounding(CbcModel & model)
1721        : CbcHeuristic(model)
1722{
1723    // Get a copy of original matrix (and by row for rounding);
1724    assert(model.solver());
1725    if (model.solver()->getNumRows()) {
1726        matrix_ = *model.solver()->getMatrixByCol();
1727        matrixByRow_ = *model.solver()->getMatrixByRow();
1728        validate();
1729    }
1730    down_ = NULL;
1731    up_ = NULL;
1732    equal_ = NULL;
1733    seed_ = 7654321;
1734    //whereFrom_ |= 16; // allow more often
1735}
1736
1737// Destructor
1738CbcRounding::~CbcRounding ()
1739{
1740    delete [] down_;
1741    delete [] up_;
1742    delete [] equal_;
1743}
1744
1745// Clone
1746CbcHeuristic *
1747CbcRounding::clone() const
1748{
1749    return new CbcRounding(*this);
1750}
1751// Create C++ lines to get to current state
1752void
1753CbcRounding::generateCpp( FILE * fp)
1754{
1755    CbcRounding other;
1756    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
1757    fprintf(fp, "3  CbcRounding rounding(*cbcModel);\n");
1758    CbcHeuristic::generateCpp(fp, "rounding");
1759    if (seed_ != other.seed_)
1760        fprintf(fp, "3  rounding.setSeed(%d);\n", seed_);
1761    else
1762        fprintf(fp, "4  rounding.setSeed(%d);\n", seed_);
1763    fprintf(fp, "3  cbcModel->addHeuristic(&rounding);\n");
1764}
1765//#define NEW_ROUNDING
1766// Copy constructor
1767CbcRounding::CbcRounding(const CbcRounding & rhs)
1768        :
1769        CbcHeuristic(rhs),
1770        matrix_(rhs.matrix_),
1771        matrixByRow_(rhs.matrixByRow_),
1772        seed_(rhs.seed_)
1773{
1774#ifdef NEW_ROUNDING
1775    int numberColumns = matrix_.getNumCols();
1776    down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1777    up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1778    equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
1779#else
1780    down_ = NULL;
1781    up_ = NULL;
1782    equal_ = NULL;
1783#endif
1784}
1785
1786// Assignment operator
1787CbcRounding &
1788CbcRounding::operator=( const CbcRounding & rhs)
1789{
1790    if (this != &rhs) {
1791        CbcHeuristic::operator=(rhs);
1792        matrix_ = rhs.matrix_;
1793        matrixByRow_ = rhs.matrixByRow_;
1794#ifdef NEW_ROUNDING
1795        delete [] down_;
1796        delete [] up_;
1797        delete [] equal_;
1798        int numberColumns = matrix_.getNumCols();
1799        down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1800        up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1801        equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
1802#else
1803        down_ = NULL;
1804        up_ = NULL;
1805        equal_ = NULL;
1806#endif
1807        seed_ = rhs.seed_;
1808    }
1809    return *this;
1810}
1811
1812// Resets stuff if model changes
1813void
1814CbcRounding::resetModel(CbcModel * model)
1815{
1816    model_ = model;
1817    // Get a copy of original matrix (and by row for rounding);
1818    assert(model_->solver());
1819    matrix_ = *model_->solver()->getMatrixByCol();
1820    matrixByRow_ = *model_->solver()->getMatrixByRow();
1821    validate();
1822}
1823// See if rounding will give solution
1824// Sets value of solution
1825// Assumes rhs for original matrix still okay
1826// At present only works with integers
1827// Fix values if asked for
1828// Returns 1 if solution, 0 if not
1829int
1830CbcRounding::solution(double & solutionValue,
1831                      double * betterSolution)
1832{
1833
1834    numCouldRun_++;
1835    // See if to do
1836    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
1837            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
1838        return 0; // switched off
1839    numRuns_++;
1840    OsiSolverInterface * solver = model_->solver();
1841    double direction = solver->getObjSense();
1842    double newSolutionValue = direction * solver->getObjValue();
1843    return solution(solutionValue, betterSolution, newSolutionValue);
1844}
1845// See if rounding will give solution
1846// Sets value of solution
1847// Assumes rhs for original matrix still okay
1848// At present only works with integers
1849// Fix values if asked for
1850// Returns 1 if solution, 0 if not
1851int
1852CbcRounding::solution(double & solutionValue,
1853                      double * betterSolution,
1854                      double newSolutionValue)
1855{
1856
1857    // See if to do
1858    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
1859            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
1860        return 0; // switched off
1861    OsiSolverInterface * solver = model_->solver();
1862    const double * lower = solver->getColLower();
1863    const double * upper = solver->getColUpper();
1864    const double * rowLower = solver->getRowLower();
1865    const double * rowUpper = solver->getRowUpper();
1866    const double * solution = solver->getColSolution();
1867    const double * objective = solver->getObjCoefficients();
1868    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1869    double primalTolerance;
1870    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1871    double useTolerance = primalTolerance;
1872
1873    int numberRows = matrix_.getNumRows();
1874    assert (numberRows <= solver->getNumRows());
1875    int numberIntegers = model_->numberIntegers();
1876    const int * integerVariable = model_->integerVariable();
1877    int i;
1878    double direction = solver->getObjSense();
1879    //double newSolutionValue = direction*solver->getObjValue();
1880    int returnCode = 0;
1881    // Column copy
1882    const double * element = matrix_.getElements();
1883    const int * row = matrix_.getIndices();
1884    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1885    const int * columnLength = matrix_.getVectorLengths();
1886    // Row copy
1887    const double * elementByRow = matrixByRow_.getElements();
1888    const int * column = matrixByRow_.getIndices();
1889    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
1890    const int * rowLength = matrixByRow_.getVectorLengths();
1891
1892    // Get solution array for heuristic solution
1893    int numberColumns = solver->getNumCols();
1894    double * newSolution = new double [numberColumns];
1895    memcpy(newSolution, solution, numberColumns*sizeof(double));
1896
1897    double * rowActivity = new double[numberRows];
1898    memset(rowActivity, 0, numberRows*sizeof(double));
1899    for (i = 0; i < numberColumns; i++) {
1900        int j;
1901        double value = newSolution[i];
1902        if (value < lower[i]) {
1903            value = lower[i];
1904            newSolution[i] = value;
1905        } else if (value > upper[i]) {
1906            value = upper[i];
1907            newSolution[i] = value;
1908        }
1909        if (value) {
1910            for (j = columnStart[i];
1911                    j < columnStart[i] + columnLength[i]; j++) {
1912                int iRow = row[j];
1913                rowActivity[iRow] += value * element[j];
1914            }
1915        }
1916    }
1917    // check was feasible - if not adjust (cleaning may move)
1918    for (i = 0; i < numberRows; i++) {
1919        if (rowActivity[i] < rowLower[i]) {
1920            //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
1921            rowActivity[i] = rowLower[i];
1922        } else if (rowActivity[i] > rowUpper[i]) {
1923            //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
1924            rowActivity[i] = rowUpper[i];
1925        }
1926    }
1927    for (i = 0; i < numberIntegers; i++) {
1928        int iColumn = integerVariable[i];
1929        double value = newSolution[iColumn];
1930        double thisTolerance = integerTolerance;
1931        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
1932            double below = floor(value);
1933            double newValue = newSolution[iColumn];
1934            double cost = direction * objective[iColumn];
1935            double move;
1936            if (cost > 0.0) {
1937                // try up
1938                move = 1.0 - (value - below);
1939            } else if (cost < 0.0) {
1940                // try down
1941                move = below - value;
1942            } else {
1943                // won't be able to move unless we can grab another variable
1944                double randomNumber = randomNumberGenerator_.randomDouble();
1945                // which way?
1946                if (randomNumber < 0.5)
1947                    move = below - value;
1948                else
1949                    move = 1.0 - (value - below);
1950            }
1951            newValue += move;
1952            newSolution[iColumn] = newValue;
1953            newSolutionValue += move * cost;
1954            int j;
1955            for (j = columnStart[iColumn];
1956                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1957                int iRow = row[j];
1958                rowActivity[iRow] += move * element[j];
1959            }
1960        }
1961    }
1962
1963    double penalty = 0.0;
1964    const char * integerType = model_->integerType();
1965    // see if feasible - just using singletons
1966    for (i = 0; i < numberRows; i++) {
1967        double value = rowActivity[i];
1968        double thisInfeasibility = 0.0;
1969        if (value < rowLower[i] - primalTolerance)
1970            thisInfeasibility = value - rowLower[i];
1971        else if (value > rowUpper[i] + primalTolerance)
1972            thisInfeasibility = value - rowUpper[i];
1973        if (thisInfeasibility) {
1974            // See if there are any slacks I can use to fix up
1975            // maybe put in coding for multiple slacks?
1976            double bestCost = 1.0e50;
1977            int k;
1978            int iBest = -1;
1979            double addCost = 0.0;
1980            double newValue = 0.0;
1981            double changeRowActivity = 0.0;
1982            double absInfeasibility = fabs(thisInfeasibility);
1983            for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
1984                int iColumn = column[k];
1985                // See if all elements help
1986                if (columnLength[iColumn] == 1) {
1987                    double currentValue = newSolution[iColumn];
1988                    double elementValue = elementByRow[k];
1989                    double lowerValue = lower[iColumn];
1990                    double upperValue = upper[iColumn];
1991                    double gap = rowUpper[i] - rowLower[i];
1992                    double absElement = fabs(elementValue);
1993                    if (thisInfeasibility*elementValue > 0.0) {
1994                        // we want to reduce
1995                        if ((currentValue - lowerValue)*absElement >= absInfeasibility) {
1996                            // possible - check if integer
1997                            double distance = absInfeasibility / absElement;
1998                            double thisCost = -direction * objective[iColumn] * distance;
1999                            if (integerType[iColumn]) {
2000                                distance = ceil(distance - useTolerance);
2001                                if (currentValue - distance >= lowerValue - useTolerance) {
2002                                    if (absInfeasibility - distance*absElement < -gap - useTolerance)
2003                                        thisCost = 1.0e100; // no good
2004                                    else
2005                                        thisCost = -direction * objective[iColumn] * distance;
2006                                } else {
2007                                    thisCost = 1.0e100; // no good
2008                                }
2009                            }
2010                            if (thisCost < bestCost) {
2011                                bestCost = thisCost;
2012                                iBest = iColumn;
2013                                addCost = thisCost;
2014                                newValue = currentValue - distance;
2015                                changeRowActivity = -distance * elementValue;
2016                            }
2017                        }
2018                    } else {
2019                        // we want to increase
2020                        if ((upperValue - currentValue)*absElement >= absInfeasibility) {
2021                            // possible - check if integer
2022                            double distance = absInfeasibility / absElement;
2023                            double thisCost = direction * objective[iColumn] * distance;
2024                            if (integerType[iColumn]) {
2025                                distance = ceil(distance - 1.0e-7);
2026                                assert (currentValue - distance <= upperValue + useTolerance);
2027                                if (absInfeasibility - distance*absElement < -gap - useTolerance)
2028                                    thisCost = 1.0e100; // no good
2029                                else
2030                                    thisCost = direction * objective[iColumn] * distance;
2031                            }
2032                            if (thisCost < bestCost) {
2033                                bestCost = thisCost;
2034                                iBest = iColumn;
2035                                addCost = thisCost;
2036                                newValue = currentValue + distance;
2037                                changeRowActivity = distance * elementValue;
2038                            }
2039                        }
2040                    }
2041                }
2042            }
2043            if (iBest >= 0) {
2044                /*printf("Infeasibility of %g on row %d cost %g\n",
2045                  thisInfeasibility,i,addCost);*/
2046                newSolution[iBest] = newValue;
2047                thisInfeasibility = 0.0;
2048                newSolutionValue += addCost;
2049                rowActivity[i] += changeRowActivity;
2050            }
2051            penalty += fabs(thisInfeasibility);
2052        }
2053    }
2054    if (penalty) {
2055        // see if feasible using any
2056        // first continuous
2057        double penaltyChange = 0.0;
2058        int iColumn;
2059        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2060            if (integerType[iColumn])
2061                continue;
2062            double currentValue = newSolution[iColumn];
2063            double lowerValue = lower[iColumn];
2064            double upperValue = upper[iColumn];
2065            int j;
2066            int anyBadDown = 0;
2067            int anyBadUp = 0;
2068            double upImprovement = 0.0;
2069            double downImprovement = 0.0;
2070            for (j = columnStart[iColumn];
2071                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2072                int iRow = row[j];
2073                if (rowUpper[iRow] > rowLower[iRow]) {
2074                    double value = element[j];
2075                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2076                        // infeasible above
2077                        downImprovement += value;
2078                        upImprovement -= value;
2079                        if (value > 0.0)
2080                            anyBadUp++;
2081                        else
2082                            anyBadDown++;
2083                    } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2084                        // feasible at ub
2085                        if (value > 0.0) {
2086                            upImprovement -= value;
2087                            anyBadUp++;
2088                        } else {
2089                            downImprovement += value;
2090                            anyBadDown++;
2091                        }
2092                    } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2093                        // feasible in interior
2094                    } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2095                        // feasible at lb
2096                        if (value < 0.0) {
2097                            upImprovement += value;
2098                            anyBadUp++;
2099                        } else {
2100                            downImprovement -= value;
2101                            anyBadDown++;
2102                        }
2103                    } else {
2104                        // infeasible below
2105                        downImprovement -= value;
2106                        upImprovement += value;
2107                        if (value < 0.0)
2108                            anyBadUp++;
2109                        else
2110                            anyBadDown++;
2111                    }
2112                } else {
2113                    // equality row
2114                    double value = element[j];
2115                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2116                        // infeasible above
2117                        downImprovement += value;
2118                        upImprovement -= value;
2119                        if (value > 0.0)
2120                            anyBadUp++;
2121                        else
2122                            anyBadDown++;
2123                    } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2124                        // infeasible below
2125                        downImprovement -= value;
2126                        upImprovement += value;
2127                        if (value < 0.0)
2128                            anyBadUp++;
2129                        else
2130                            anyBadDown++;
2131                    } else {
2132                        // feasible - no good
2133                        anyBadUp = -1;
2134                        anyBadDown = -1;
2135                        break;
2136                    }
2137                }
2138            }
2139            // could change tests for anyBad
2140            if (anyBadUp)
2141                upImprovement = 0.0;
2142            if (anyBadDown)
2143                downImprovement = 0.0;
2144            double way = 0.0;
2145            double improvement = 0.0;
2146            if (downImprovement > 0.0 && currentValue > lowerValue) {
2147                way = -1.0;
2148                improvement = downImprovement;
2149            } else if (upImprovement > 0.0 && currentValue < upperValue) {
2150                way = 1.0;
2151                improvement = upImprovement;
2152            }
2153            if (way) {
2154                // can improve
2155                double distance;
2156                if (way > 0.0)
2157                    distance = upperValue - currentValue;
2158                else
2159                    distance = currentValue - lowerValue;
2160                for (j = columnStart[iColumn];
2161                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2162                    int iRow = row[j];
2163                    double value = element[j] * way;
2164                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2165                        // infeasible above
2166                        assert (value < 0.0);
2167                        double gap = rowActivity[iRow] - rowUpper[iRow];
2168                        if (gap + value*distance < 0.0)
2169                            distance = -gap / value;
2170                    } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2171                        // infeasible below
2172                        assert (value > 0.0);
2173                        double gap = rowActivity[iRow] - rowLower[iRow];
2174                        if (gap + value*distance > 0.0)
2175                            distance = -gap / value;
2176                    } else {
2177                        // feasible
2178                        if (value > 0) {
2179                            double gap = rowActivity[iRow] - rowUpper[iRow];
2180                            if (gap + value*distance > 0.0)
2181                                distance = -gap / value;
2182                        } else {
2183                            double gap = rowActivity[iRow] - rowLower[iRow];
2184                            if (gap + value*distance < 0.0)
2185                                distance = -gap / value;
2186                        }
2187                    }
2188                }
2189                //move
2190                penaltyChange += improvement * distance;
2191                distance *= way;
2192                newSolution[iColumn] += distance;
2193                newSolutionValue += direction * objective[iColumn] * distance;
2194                for (j = columnStart[iColumn];
2195                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2196                    int iRow = row[j];
2197                    double value = element[j];
2198                    rowActivity[iRow] += distance * value;
2199                }
2200            }
2201        }
2202        // and now all if improving
2203        double lastChange = penaltyChange ? 1.0 : 0.0;
2204        int numberPasses=0;
2205        while (lastChange > 1.0e-2 && numberPasses < 1000) {
2206            lastChange = 0;
2207            numberPasses++;
2208            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2209                bool isInteger = (integerType[iColumn] != 0);
2210                double currentValue = newSolution[iColumn];
2211                double lowerValue = lower[iColumn];
2212                double upperValue = upper[iColumn];
2213                int j;
2214                int anyBadDown = 0;
2215                int anyBadUp = 0;
2216                double upImprovement = 0.0;
2217                double downImprovement = 0.0;
2218                for (j = columnStart[iColumn];
2219                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2220                    int iRow = row[j];
2221                    double value = element[j];
2222                    if (isInteger) {
2223                        if (value > 0.0) {
2224                            if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance)
2225                                anyBadUp++;
2226                            if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance)
2227                                anyBadDown++;
2228                        } else {
2229                            if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance)
2230                                anyBadDown++;
2231                            if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance)
2232                                anyBadUp++;
2233                        }
2234                    }
2235                    if (rowUpper[iRow] > rowLower[iRow]) {
2236                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2237                            // infeasible above
2238                            downImprovement += value;
2239                            upImprovement -= value;
2240                            if (value > 0.0)
2241                                anyBadUp++;
2242                            else
2243                                anyBadDown++;
2244                        } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2245                            // feasible at ub
2246                            if (value > 0.0) {
2247                                upImprovement -= value;
2248                                anyBadUp++;
2249                            } else {
2250                                downImprovement += value;
2251                                anyBadDown++;
2252                            }
2253                        } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2254                            // feasible in interior
2255                        } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2256                            // feasible at lb
2257                            if (value < 0.0) {
2258                                upImprovement += value;
2259                                anyBadUp++;
2260                            } else {
2261                                downImprovement -= value;
2262                                anyBadDown++;
2263                            }
2264                        } else {
2265                            // infeasible below
2266                            downImprovement -= value;
2267                            upImprovement += value;
2268                            if (value < 0.0)
2269                                anyBadUp++;
2270                            else
2271                                anyBadDown++;
2272                        }
2273                    } else {
2274                        // equality row
2275                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2276                            // infeasible above
2277                            downImprovement += value;
2278                            upImprovement -= value;
2279                            if (value > 0.0)
2280                                anyBadUp++;
2281                            else
2282                                anyBadDown++;
2283                        } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2284                            // infeasible below
2285                            downImprovement -= value;
2286                            upImprovement += value;
2287                            if (value < 0.0)
2288                                anyBadUp++;
2289                            else
2290                                anyBadDown++;
2291                        } else {
2292                            // feasible - no good
2293                            anyBadUp = -1;
2294                            anyBadDown = -1;
2295                            break;
2296                        }
2297                    }
2298                }
2299                // could change tests for anyBad
2300                if (anyBadUp)
2301                    upImprovement = 0.0;
2302                if (anyBadDown)
2303                    downImprovement = 0.0;
2304                double way = 0.0;
2305                double improvement = 0.0;
2306                if (downImprovement > 0.0 && currentValue > lowerValue) {
2307                    way = -1.0;
2308                    improvement = downImprovement;
2309                } else if (upImprovement > 0.0 && currentValue < upperValue) {
2310                    way = 1.0;
2311                    improvement = upImprovement;
2312                }
2313                if (way) {
2314                    // can improve
2315                    double distance = COIN_DBL_MAX;
2316                    for (j = columnStart[iColumn];
2317                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2318                        int iRow = row[j];
2319                        double value = element[j] * way;
2320                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2321                            // infeasible above
2322                            assert (value < 0.0);
2323                            double gap = rowActivity[iRow] - rowUpper[iRow];
2324                            if (gap + value*distance < 0.0) {
2325                                // If integer then has to move by 1
2326                                if (!isInteger)
2327                                    distance = -gap / value;
2328                                else
2329                                    distance = CoinMax(-gap / value, 1.0);
2330                            }
2331                        } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2332                            // infeasible below
2333                            assert (value > 0.0);
2334                            double gap = rowActivity[iRow] - rowLower[iRow];
2335                            if (gap + value*distance > 0.0) {
2336                                // If integer then has to move by 1
2337                                if (!isInteger)
2338                                    distance = -gap / value;
2339                                else
2340                                    distance = CoinMax(-gap / value, 1.0);
2341                            }
2342                        } else {
2343                            // feasible
2344                            if (value > 0) {
2345                                double gap = rowActivity[iRow] - rowUpper[iRow];
2346                                if (gap + value*distance > 0.0)
2347                                    distance = -gap / value;
2348                            } else {
2349                                double gap = rowActivity[iRow] - rowLower[iRow];
2350                                if (gap + value*distance < 0.0)
2351                                    distance = -gap / value;
2352                            }
2353                        }
2354                    }
2355                    if (isInteger)
2356                        distance = floor(distance + 1.05e-8);
2357                    if (!distance) {
2358                        // should never happen
2359                        //printf("zero distance in CbcRounding - debug\n");
2360                    }
2361                    //move
2362                    lastChange += improvement * distance;
2363                    distance *= way;
2364                    newSolution[iColumn] += distance;
2365                    newSolutionValue += direction * objective[iColumn] * distance;
2366                    for (j = columnStart[iColumn];
2367                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2368                        int iRow = row[j];
2369                        double value = element[j];
2370                        rowActivity[iRow] += distance * value;
2371                    }
2372                }
2373            }
2374            penaltyChange += lastChange;
2375        }
2376        penalty -= penaltyChange;
2377        if (penalty < 1.0e-5*fabs(penaltyChange)) {
2378            // recompute
2379            penalty = 0.0;
2380            for (i = 0; i < numberRows; i++) {
2381                double value = rowActivity[i];
2382                if (value < rowLower[i] - primalTolerance)
2383                    penalty += rowLower[i] - value;
2384                else if (value > rowUpper[i] + primalTolerance)
2385                    penalty += value - rowUpper[i];
2386            }
2387        }
2388    }
2389
2390    // Could also set SOS (using random) and repeat
2391    if (!penalty) {
2392        // See if we can do better
2393        //seed_++;
2394        //CoinSeedRandom(seed_);
2395        // Random number between 0 and 1.
2396        double randomNumber = randomNumberGenerator_.randomDouble();
2397        int iPass;
2398        int start[2];
2399        int end[2];
2400        int iRandom = static_cast<int> (randomNumber * (static_cast<double> (numberIntegers)));
2401        start[0] = iRandom;
2402        end[0] = numberIntegers;
2403        start[1] = 0;
2404        end[1] = iRandom;
2405        for (iPass = 0; iPass < 2; iPass++) {
2406            int i;
2407            for (i = start[iPass]; i < end[iPass]; i++) {
2408                int iColumn = integerVariable[i];
2409#ifndef NDEBUG
2410                double value = newSolution[iColumn];
2411                assert (fabs(floor(value + 0.5) - value) < integerTolerance);
2412#endif
2413                double cost = direction * objective[iColumn];
2414                double move = 0.0;
2415                if (cost > 0.0)
2416                    move = -1.0;
2417                else if (cost < 0.0)
2418                    move = 1.0;
2419                while (move) {
2420                    bool good = true;
2421                    double newValue = newSolution[iColumn] + move;
2422                    if (newValue < lower[iColumn] - useTolerance ||
2423                            newValue > upper[iColumn] + useTolerance) {
2424                        move = 0.0;
2425                    } else {
2426                        // see if we can move
2427                        int j;
2428                        for (j = columnStart[iColumn];
2429                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2430                            int iRow = row[j];
2431                            double newActivity = rowActivity[iRow] + move * element[j];
2432                            if (newActivity < rowLower[iRow] - primalTolerance ||
2433                                    newActivity > rowUpper[iRow] + primalTolerance) {
2434                                good = false;
2435                                break;
2436                            }
2437                        }
2438                        if (good) {
2439                            newSolution[iColumn] = newValue;
2440                            newSolutionValue += move * cost;
2441                            int j;
2442                            for (j = columnStart[iColumn];
2443                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2444                                int iRow = row[j];
2445                                rowActivity[iRow] += move * element[j];
2446                            }
2447                        } else {
2448                            move = 0.0;
2449                        }
2450                    }
2451                }
2452            }
2453        }
2454        // Just in case of some stupidity
2455        double objOffset = 0.0;
2456        solver->getDblParam(OsiObjOffset, objOffset);
2457        newSolutionValue = -objOffset;
2458        for ( i = 0 ; i < numberColumns ; i++ )
2459            newSolutionValue += objective[i] * newSolution[i];
2460        newSolutionValue *= direction;
2461        //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
2462        if (newSolutionValue < solutionValue) {
2463            // paranoid check
2464            memset(rowActivity, 0, numberRows*sizeof(double));
2465            for (i = 0; i < numberColumns; i++) {
2466                int j;
2467                double value = newSolution[i];
2468                if (value) {
2469                    for (j = columnStart[i];
2470                            j < columnStart[i] + columnLength[i]; j++) {
2471                        int iRow = row[j];
2472                        rowActivity[iRow] += value * element[j];
2473                    }
2474                }
2475            }
2476            // check was approximately feasible
2477            bool feasible = true;
2478            for (i = 0; i < numberRows; i++) {
2479                if (rowActivity[i] < rowLower[i]) {
2480                    if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
2481                        feasible = false;
2482                } else if (rowActivity[i] > rowUpper[i]) {
2483                    if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
2484                        feasible = false;
2485                }
2486            }
2487            if (feasible) {
2488                // new solution
2489                memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2490                solutionValue = newSolutionValue;
2491                //printf("** Solution of %g found by rounding\n",newSolutionValue);
2492                returnCode = 1;
2493            } else {
2494                // Can easily happen
2495                //printf("Debug CbcRounding giving bad solution\n");
2496            }
2497        }
2498    }
2499#ifdef NEW_ROUNDING
2500    if (!returnCode) {
2501#ifdef JJF_ZERO
2502        // back to starting point
2503        memcpy(newSolution, solution, numberColumns*sizeof(double));
2504        memset(rowActivity, 0, numberRows*sizeof(double));
2505        for (i = 0; i < numberColumns; i++) {
2506            int j;
2507            double value = newSolution[i];
2508            if (value < lower[i]) {
2509                value = lower[i];
2510                newSolution[i] = value;
2511            } else if (value > upper[i]) {
2512                value = upper[i];
2513                newSolution[i] = value;
2514            }
2515            if (value) {
2516                for (j = columnStart[i];
2517                        j < columnStart[i] + columnLength[i]; j++) {
2518                    int iRow = row[j];
2519                    rowActivity[iRow] += value * element[j];
2520                }
2521            }
2522        }
2523        // check was feasible - if not adjust (cleaning may move)
2524        for (i = 0; i < numberRows; i++) {
2525            if (rowActivity[i] < rowLower[i]) {
2526                //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
2527                rowActivity[i] = rowLower[i];
2528            } else if (rowActivity[i] > rowUpper[i]) {
2529                //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
2530                rowActivity[i] = rowUpper[i];
2531            }
2532        }
2533#endif
2534        int * candidate = new int [numberColumns];
2535        int nCandidate = 0;
2536        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2537            bool isInteger = (integerType[iColumn] != 0);
2538            if (isInteger) {
2539                double currentValue = newSolution[iColumn];
2540                if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8)
2541                    candidate[nCandidate++] = iColumn;
2542            }
2543        }
2544        if (true) {
2545            // Rounding as in Berthold
2546            while (nCandidate) {
2547                double infeasibility = 1.0e-7;
2548                int iRow = -1;
2549                for (i = 0; i < numberRows; i++) {
2550                    double value = 0.0;
2551                    if (rowActivity[i] < rowLower[i]) {
2552                        value = rowLower[i] - rowActivity[i];
2553                    } else if (rowActivity[i] > rowUpper[i]) {
2554                        value = rowActivity[i] - rowUpper[i];
2555                    }
2556                    if (value > infeasibility) {
2557                        infeasibility = value;
2558                        iRow = i;
2559                    }
2560                }
2561                if (iRow >= 0) {
2562                    // infeasible
2563                } else {
2564                    // feasible
2565                }
2566            }
2567        } else {
2568            // Shifting as in Berthold
2569        }
2570        delete [] candidate;
2571    }
2572#endif
2573    delete [] newSolution;
2574    delete [] rowActivity;
2575    return returnCode;
2576}
2577// update model
2578void CbcRounding::setModel(CbcModel * model)
2579{
2580    model_ = model;
2581    // Get a copy of original matrix (and by row for rounding);
2582    assert(model_->solver());
2583    if (model_->solver()->getNumRows()) {
2584        matrix_ = *model_->solver()->getMatrixByCol();
2585        matrixByRow_ = *model_->solver()->getMatrixByRow();
2586        // make sure model okay for heuristic
2587        validate();
2588    }
2589}
2590// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
2591void
2592CbcRounding::validate()
2593{
2594    if (model_ && (when() % 100) < 10) {
2595        if (model_->numberIntegers() !=
2596                model_->numberObjects() && (model_->numberObjects() ||
2597                                            (model_->specialOptions()&1024) == 0)) {
2598            int numberOdd = 0;
2599            for (int i = 0; i < model_->numberObjects(); i++) {
2600                if (!model_->object(i)->canDoHeuristics())
2601                    numberOdd++;
2602            }
2603            if (numberOdd)
2604                setWhen(0);
2605        }
2606    }
2607#ifdef NEW_ROUNDING
2608    int numberColumns = matrix_.getNumCols();
2609    down_ = new unsigned short [numberColumns];
2610    up_ = new unsigned short [numberColumns];
2611    equal_ = new unsigned short [numberColumns];
2612    // Column copy
2613    const double * element = matrix_.getElements();
2614    const int * row = matrix_.getIndices();
2615    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2616    const int * columnLength = matrix_.getVectorLengths();
2617    const double * rowLower = model.solver()->getRowLower();
2618    const double * rowUpper = model.solver()->getRowUpper();
2619    for (int i = 0; i < numberColumns; i++) {
2620        int down = 0;
2621        int up = 0;
2622        int equal = 0;
2623        if (columnLength[i] > 65535) {
2624            equal[0] = 65535;
2625            break; // unlikely to work
2626        }
2627        for (CoinBigIndex j = columnStart[i];
2628                j < columnStart[i] + columnLength[i]; j++) {
2629            int iRow = row[j];
2630            if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
2631                equal++;
2632            } else if (element[j] > 0.0) {
2633                if (rowUpper[iRow] < 1.0e20)
2634                    up++;
2635                else
2636                    down--;
2637            } else {
2638                if (rowLower[iRow] > -1.0e20)
2639                    up++;
2640                else
2641                    down--;
2642            }
2643        }
2644        down_[i] = (unsigned short) down;
2645        up_[i] = (unsigned short) up;
2646        equal_[i] = (unsigned short) equal;
2647    }
2648#else
2649    down_ = NULL;
2650    up_ = NULL;
2651    equal_ = NULL;
2652#endif
2653}
2654
2655// Default Constructor
2656CbcHeuristicPartial::CbcHeuristicPartial()
2657        : CbcHeuristic()
2658{
2659    fixPriority_ = 10000;
2660}
2661
2662// Constructor from model
2663CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
2664        : CbcHeuristic(model)
2665{
2666    fixPriority_ = fixPriority;
2667    setNumberNodes(numberNodes);
2668    validate();
2669}
2670
2671// Destructor
2672CbcHeuristicPartial::~CbcHeuristicPartial ()
2673{
2674}
2675
2676// Clone
2677CbcHeuristic *
2678CbcHeuristicPartial::clone() const
2679{
2680    return new CbcHeuristicPartial(*this);
2681}
2682// Create C++ lines to get to current state
2683void
2684CbcHeuristicPartial::generateCpp( FILE * fp)
2685{
2686    CbcHeuristicPartial other;
2687    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2688    fprintf(fp, "3  CbcHeuristicPartial partial(*cbcModel);\n");
2689    CbcHeuristic::generateCpp(fp, "partial");
2690    if (fixPriority_ != other.fixPriority_)
2691        fprintf(fp, "3  partial.setFixPriority(%d);\n", fixPriority_);
2692    else
2693        fprintf(fp, "4  partial.setFixPriority(%d);\n", fixPriority_);
2694    fprintf(fp, "3  cbcModel->addHeuristic(&partial);\n");
2695}
2696//#define NEW_PARTIAL
2697// Copy constructor
2698CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
2699        :
2700        CbcHeuristic(rhs),
2701        fixPriority_(rhs.fixPriority_)
2702{
2703}
2704
2705// Assignment operator
2706CbcHeuristicPartial &
2707CbcHeuristicPartial::operator=( const CbcHeuristicPartial & rhs)
2708{
2709    if (this != &rhs) {
2710        CbcHeuristic::operator=(rhs);
2711        fixPriority_ = rhs.fixPriority_;
2712    }
2713    return *this;
2714}
2715
2716// Resets stuff if model changes
2717void
2718CbcHeuristicPartial::resetModel(CbcModel * model)
2719{
2720    model_ = model;
2721    // Get a copy of original matrix (and by row for partial);
2722    assert(model_->solver());
2723    validate();
2724}
2725// See if partial will give solution
2726// Sets value of solution
2727// Assumes rhs for original matrix still okay
2728// At present only works with integers
2729// Fix values if asked for
2730// Returns 1 if solution, 0 if not
2731int
2732CbcHeuristicPartial::solution(double & solutionValue,
2733                              double * betterSolution)
2734{
2735    // Return if already done
2736    if (fixPriority_ < 0)
2737        return 0; // switched off
2738    const double * hotstartSolution = model_->hotstartSolution();
2739    const int * hotstartPriorities = model_->hotstartPriorities();
2740    if (!hotstartSolution)
2741        return 0;
2742    OsiSolverInterface * solver = model_->solver();
2743
2744    int numberIntegers = model_->numberIntegers();
2745    const int * integerVariable = model_->integerVariable();
2746
2747    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
2748    const double * colLower = newSolver->getColLower();
2749    const double * colUpper = newSolver->getColUpper();
2750
2751    double primalTolerance;
2752    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
2753
2754    int i;
2755    int numberFixed = 0;
2756    int returnCode = 0;
2757
2758    for (i = 0; i < numberIntegers; i++) {
2759        int iColumn = integerVariable[i];
2760        if (abs(hotstartPriorities[iColumn]) <= fixPriority_) {
2761            double value = hotstartSolution[iColumn];
2762            double lower = colLower[iColumn];
2763            double upper = colUpper[iColumn];
2764            value = CoinMax(value, lower);
2765            value = CoinMin(value, upper);
2766            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
2767                value = floor(value + 0.5);
2768                newSolver->setColLower(iColumn, value);
2769                newSolver->setColUpper(iColumn, value);
2770                numberFixed++;
2771            }
2772        }
2773    }
2774    if (numberFixed > numberIntegers / 5 - 100000000) {
2775#ifdef COIN_DEVELOP
2776        printf("%d integers fixed\n", numberFixed);
2777#endif
2778        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
2779                                         model_->getCutoff(), "CbcHeuristicPartial");
2780        if (returnCode < 0)
2781            returnCode = 0; // returned on size
2782        //printf("return code %d",returnCode);
2783        if ((returnCode&2) != 0) {
2784            // could add cut
2785            returnCode &= ~2;
2786            //printf("could add cut with %d elements (if all 0-1)\n",nFix);
2787        } else {
2788            //printf("\n");
2789        }
2790    }
2791    fixPriority_ = -1; // switch off
2792
2793    delete newSolver;
2794    return returnCode;
2795}
2796// update model
2797void CbcHeuristicPartial::setModel(CbcModel * model)
2798{
2799    model_ = model;
2800    assert(model_->solver());
2801    // make sure model okay for heuristic
2802    validate();
2803}
2804// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
2805void
2806CbcHeuristicPartial::validate()
2807{
2808    if (model_ && (when() % 100) < 10) {
2809        if (model_->numberIntegers() !=
2810                model_->numberObjects())
2811            setWhen(0);
2812    }
2813}
2814bool
2815CbcHeuristicPartial::shouldHeurRun(int /*whereFrom*/)
2816{
2817    return true;
2818}
2819
2820// Default Constructor
2821CbcSerendipity::CbcSerendipity()
2822        : CbcHeuristic()
2823{
2824}
2825
2826// Constructor from model
2827CbcSerendipity::CbcSerendipity(CbcModel & model)
2828        : CbcHeuristic(model)
2829{
2830}
2831
2832// Destructor
2833CbcSerendipity::~CbcSerendipity ()
2834{
2835}
2836
2837// Clone
2838CbcHeuristic *
2839CbcSerendipity::clone() const
2840{
2841    return new CbcSerendipity(*this);
2842}
2843// Create C++ lines to get to current state
2844void
2845CbcSerendipity::generateCpp( FILE * fp)
2846{
2847    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2848    fprintf(fp, "3  CbcSerendipity serendipity(*cbcModel);\n");
2849    CbcHeuristic::generateCpp(fp, "serendipity");
2850    fprintf(fp, "3  cbcModel->addHeuristic(&serendipity);\n");
2851}
2852
2853// Copy constructor
2854CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
2855        :
2856        CbcHeuristic(rhs)
2857{
2858}
2859
2860// Assignment operator
2861CbcSerendipity &
2862CbcSerendipity::operator=( const CbcSerendipity & rhs)
2863{
2864    if (this != &rhs) {
2865        CbcHeuristic::operator=(rhs);
2866    }
2867    return *this;
2868}
2869
2870// Returns 1 if solution, 0 if not
2871int
2872CbcSerendipity::solution(double & solutionValue,
2873                         double * betterSolution)
2874{
2875    if (!model_)
2876        return 0;
2877    if (!inputSolution_) {
2878        // get information on solver type
2879        OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
2880        OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
2881        if (auxiliaryInfo) {
2882            return auxiliaryInfo->solution(solutionValue, betterSolution, model_->solver()->getNumCols());
2883        } else {
2884            return 0;
2885        }
2886    } else {
2887        int numberColumns = model_->getNumCols();
2888        double value = inputSolution_[numberColumns];
2889        int returnCode = 0;
2890        if (value < solutionValue) {
2891            solutionValue = value;
2892            memcpy(betterSolution, inputSolution_, numberColumns*sizeof(double));
2893            returnCode = 1;
2894        }
2895        delete [] inputSolution_;
2896        inputSolution_ = NULL;
2897        model_ = NULL; // switch off
2898        return returnCode;
2899    }
2900}
2901// update model
2902void CbcSerendipity::setModel(CbcModel * model)
2903{
2904    model_ = model;
2905}
2906// Resets stuff if model changes
2907void
2908CbcSerendipity::resetModel(CbcModel * model)
2909{
2910    model_ = model;
2911}
2912
2913
2914// Default Constructor
2915CbcHeuristicJustOne::CbcHeuristicJustOne()
2916        : CbcHeuristic(),
2917        probabilities_(NULL),
2918        heuristic_(NULL),
2919        numberHeuristics_(0)
2920{
2921}
2922
2923// Constructor from model
2924CbcHeuristicJustOne::CbcHeuristicJustOne(CbcModel & model)
2925        : CbcHeuristic(model),
2926        probabilities_(NULL),
2927        heuristic_(NULL),
2928        numberHeuristics_(0)
2929{
2930}
2931
2932// Destructor
2933CbcHeuristicJustOne::~CbcHeuristicJustOne ()
2934{
2935    for (int i = 0; i < numberHeuristics_; i++)
2936        delete heuristic_[i];
2937    delete [] heuristic_;
2938    delete [] probabilities_;
2939}
2940
2941// Clone
2942CbcHeuristicJustOne *
2943CbcHeuristicJustOne::clone() const
2944{
2945    return new CbcHeuristicJustOne(*this);
2946}
2947
2948// Create C++ lines to get to current state
2949void
2950CbcHeuristicJustOne::generateCpp( FILE * fp)
2951{
2952    CbcHeuristicJustOne other;
2953    fprintf(fp, "0#include \"CbcHeuristicJustOne.hpp\"\n");
2954    fprintf(fp, "3  CbcHeuristicJustOne heuristicJustOne(*cbcModel);\n");
2955    CbcHeuristic::generateCpp(fp, "heuristicJustOne");
2956    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicJustOne);\n");
2957}
2958
2959// Copy constructor
2960CbcHeuristicJustOne::CbcHeuristicJustOne(const CbcHeuristicJustOne & rhs)
2961        :
2962        CbcHeuristic(rhs),
2963        probabilities_(NULL),
2964        heuristic_(NULL),
2965        numberHeuristics_(rhs.numberHeuristics_)
2966{
2967    if (numberHeuristics_) {
2968        probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
2969        heuristic_ = new CbcHeuristic * [numberHeuristics_];
2970        for (int i = 0; i < numberHeuristics_; i++)
2971            heuristic_[i] = rhs.heuristic_[i]->clone();
2972    }
2973}
2974
2975// Assignment operator
2976CbcHeuristicJustOne &
2977CbcHeuristicJustOne::operator=( const CbcHeuristicJustOne & rhs)
2978{
2979    if (this != &rhs) {
2980        CbcHeuristic::operator=(rhs);
2981        for (int i = 0; i < numberHeuristics_; i++)
2982            delete heuristic_[i];
2983        delete [] heuristic_;
2984        delete [] probabilities_;
2985        probabilities_ = NULL;
2986        heuristic_ = NULL;
2987        numberHeuristics_ = rhs.numberHeuristics_;
2988        if (numberHeuristics_) {
2989            probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
2990            heuristic_ = new CbcHeuristic * [numberHeuristics_];
2991            for (int i = 0; i < numberHeuristics_; i++)
2992                heuristic_[i] = rhs.heuristic_[i]->clone();
2993        }
2994    }
2995    return *this;
2996}
2997// Sets value of solution
2998// Returns 1 if solution, 0 if not
2999int
3000CbcHeuristicJustOne::solution(double & solutionValue,
3001                              double * betterSolution)
3002{
3003#ifdef DIVE_DEBUG
3004    std::cout << "solutionValue = " << solutionValue << std::endl;
3005#endif
3006    ++numCouldRun_;
3007
3008    // test if the heuristic can run
3009    if (!shouldHeurRun_randomChoice() || !numberHeuristics_)
3010        return 0;
3011    double randomNumber = randomNumberGenerator_.randomDouble();
3012    int i;
3013    for (i = 0; i < numberHeuristics_; i++) {
3014        if (randomNumber < probabilities_[i])
3015            break;
3016    }
3017    assert (i < numberHeuristics_);
3018    int returnCode;
3019    //model_->unsetDivingHasRun();
3020#ifdef COIN_DEVELOP
3021    printf("JustOne running %s\n",
3022           heuristic_[i]->heuristicName());
3023#endif
3024    returnCode = heuristic_[i]->solution(solutionValue, betterSolution);
3025#ifdef COIN_DEVELOP
3026    if (returnCode)
3027        printf("JustOne running %s found solution\n",
3028               heuristic_[i]->heuristicName());
3029#endif
3030    return returnCode;
3031}
3032// Resets stuff if model changes
3033void
3034CbcHeuristicJustOne::resetModel(CbcModel * model)
3035{
3036    CbcHeuristic::resetModel(model);
3037    for (int i = 0; i < numberHeuristics_; i++)
3038        heuristic_[i]->resetModel(model);
3039}
3040// update model (This is needed if cliques update matrix etc)
3041void
3042CbcHeuristicJustOne::setModel(CbcModel * model)
3043{
3044    CbcHeuristic::setModel(model);
3045    for (int i = 0; i < numberHeuristics_; i++)
3046        heuristic_[i]->setModel(model);
3047}
3048// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
3049void
3050CbcHeuristicJustOne::validate()
3051{
3052    CbcHeuristic::validate();
3053    for (int i = 0; i < numberHeuristics_; i++)
3054        heuristic_[i]->validate();
3055}
3056// Adds an heuristic with probability
3057void
3058CbcHeuristicJustOne::addHeuristic(const CbcHeuristic * heuristic, double probability)
3059{
3060    CbcHeuristic * thisOne = heuristic->clone();
3061    thisOne->setWhen(-999);
3062    CbcHeuristic ** tempH = CoinCopyOfArrayPartial(heuristic_, numberHeuristics_ + 1,
3063                            numberHeuristics_);
3064    delete [] heuristic_;
3065    heuristic_ = tempH;
3066    heuristic_[numberHeuristics_] = thisOne;
3067    double * tempP = CoinCopyOfArrayPartial(probabilities_, numberHeuristics_ + 1,
3068                                            numberHeuristics_);
3069    delete [] probabilities_;
3070    probabilities_ = tempP;
3071    probabilities_[numberHeuristics_] = probability;
3072    numberHeuristics_++;
3073}
3074// Normalize probabilities
3075void
3076CbcHeuristicJustOne::normalizeProbabilities()
3077{
3078    double sum = 0.0;
3079    for (int i = 0; i < numberHeuristics_; i++)
3080        sum += probabilities_[i];
3081    double multiplier = 1.0 / sum;
3082    sum = 0.0;
3083    for (int i = 0; i < numberHeuristics_; i++) {
3084        sum += probabilities_[i];
3085        probabilities_[i] = sum * multiplier;
3086    }
3087    assert (fabs(probabilities_[numberHeuristics_-1] - 1.0) < 1.0e-5);
3088    probabilities_[numberHeuristics_-1] = 1.000001;
3089}
3090
Note: See TracBrowser for help on using the repository browser.