source: branches/sandbox/Cbc/src/CbcBranchActual.cpp @ 1286

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 174.4 KB
Line 
1/* $Id: CbcBranchActual.cpp 1286 2009-11-09 23:33:07Z EdwinStraver $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12//#define CBC_DEBUG
13
14#include "CoinTypes.hpp"
15#include "OsiSolverInterface.hpp"
16#include "OsiSolverBranch.hpp"
17#include "CbcModel.hpp"
18#include "CbcMessage.hpp"
19#include "CbcBranchActual.hpp"
20#include "CoinSort.hpp"
21#include "CoinError.hpp"
22
23//##############################################################################
24
25// Default Constructor
26CbcClique::CbcClique ()
27        : CbcObject(),
28        numberMembers_(0),
29        numberNonSOSMembers_(0),
30        members_(NULL),
31        type_(NULL),
32        cliqueType_(-1),
33        slack_(-1)
34{
35}
36
37// Useful constructor (which are integer indices)
38CbcClique::CbcClique (CbcModel * model, int cliqueType, int numberMembers,
39                      const int * which, const char * type, int identifier, int slack)
40        : CbcObject(model)
41{
42    id_ = identifier;
43    numberMembers_ = numberMembers;
44    if (numberMembers_) {
45        members_ = new int[numberMembers_];
46        memcpy(members_, which, numberMembers_*sizeof(int));
47        type_ = new char[numberMembers_];
48        if (type) {
49            memcpy(type_, type, numberMembers_*sizeof(char));
50        } else {
51            for (int i = 0; i < numberMembers_; i++)
52                type_[i] = 1;
53        }
54    } else {
55        members_ = NULL;
56        type_ = NULL;
57    }
58    // Find out how many non sos
59    int i;
60    numberNonSOSMembers_ = 0;
61    for (i = 0; i < numberMembers_; i++)
62        if (!type_[i])
63            numberNonSOSMembers_++;
64    cliqueType_ = cliqueType;
65    slack_ = slack;
66}
67
68// Copy constructor
69CbcClique::CbcClique ( const CbcClique & rhs)
70        : CbcObject(rhs)
71{
72    numberMembers_ = rhs.numberMembers_;
73    numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
74    if (numberMembers_) {
75        members_ = new int[numberMembers_];
76        memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
77        type_ = new char[numberMembers_];
78        memcpy(type_, rhs.type_, numberMembers_*sizeof(char));
79    } else {
80        members_ = NULL;
81        type_ = NULL;
82    }
83    cliqueType_ = rhs.cliqueType_;
84    slack_ = rhs.slack_;
85}
86
87// Clone
88CbcObject *
89CbcClique::clone() const
90{
91    return new CbcClique(*this);
92}
93
94// Assignment operator
95CbcClique &
96CbcClique::operator=( const CbcClique & rhs)
97{
98    if (this != &rhs) {
99        CbcObject::operator=(rhs);
100        delete [] members_;
101        delete [] type_;
102        numberMembers_ = rhs.numberMembers_;
103        numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
104        if (numberMembers_) {
105            members_ = new int[numberMembers_];
106            memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
107            type_ = new char[numberMembers_];
108            memcpy(type_, rhs.type_, numberMembers_*sizeof(char));
109        } else {
110            members_ = NULL;
111            type_ = NULL;
112        }
113        cliqueType_ = rhs.cliqueType_;
114        slack_ = rhs.slack_;
115    }
116    return *this;
117}
118
119// Destructor
120CbcClique::~CbcClique ()
121{
122    delete [] members_;
123    delete [] type_;
124}
125double
126CbcClique::infeasibility(const OsiBranchingInformation * /*info*/,
127                         int &preferredWay) const
128{
129    int numberUnsatis = 0, numberFree = 0;
130    int j;
131    const int * integer = model_->integerVariable();
132    OsiSolverInterface * solver = model_->solver();
133    const double * solution = model_->testSolution();
134    const double * lower = solver->getColLower();
135    const double * upper = solver->getColUpper();
136    double largestValue = 0.0;
137    double integerTolerance =
138        model_->getDblParam(CbcModel::CbcIntegerTolerance);
139    double * sort = new double[numberMembers_];
140
141    double slackValue = 0.0;
142    for (j = 0; j < numberMembers_; j++) {
143        int sequence = members_[j];
144        int iColumn = integer[sequence];
145        double value = solution[iColumn];
146        value = CoinMax(value, lower[iColumn]);
147        value = CoinMin(value, upper[iColumn]);
148        double nearest = floor(value + 0.5);
149        double distance = fabs(value - nearest);
150        if (distance > integerTolerance) {
151            if (!type_[j])
152                value = 1.0 - value; // non SOS
153            // if slack then choose that
154            if (j == slack_ && value > 0.05)
155                slackValue = value;
156            largestValue = CoinMax(value, largestValue);
157            sort[numberUnsatis++] = -value;
158        } else if (upper[iColumn] > lower[iColumn]) {
159            numberFree++;
160        }
161    }
162    preferredWay = 1;
163    double otherWay = 0.0;
164    if (numberUnsatis) {
165        // sort
166        std::sort(sort, sort + numberUnsatis);
167        for (j = 0; j < numberUnsatis; j++) {
168            if ((j&1) != 0)
169                otherWay += -sort[j];
170        }
171        // Need to think more
172        double value = 0.2 * numberUnsatis + 0.01 * (numberMembers_ - numberFree);
173        if (fabs(largestValue - 0.5) < 0.1) {
174            // close to half
175            value += 0.1;
176        }
177        if (slackValue) {
178            // branching on slack
179            value += slackValue;
180        }
181        // scale other way
182        otherWay *= value / (1.0 - otherWay);
183        delete [] sort;
184        return value;
185    } else {
186        delete [] sort;
187        return 0.0; // satisfied
188    }
189}
190
191// This looks at solution and sets bounds to contain solution
192void
193CbcClique::feasibleRegion()
194{
195    int j;
196    const int * integer = model_->integerVariable();
197    OsiSolverInterface * solver = model_->solver();
198    const double * solution = model_->testSolution();
199    const double * lower = solver->getColLower();
200    const double * upper = solver->getColUpper();
201#ifndef NDEBUG
202    double integerTolerance =
203        model_->getDblParam(CbcModel::CbcIntegerTolerance);
204#endif
205    for (j = 0; j < numberMembers_; j++) {
206        int sequence = members_[j];
207        int iColumn = integer[sequence];
208        double value = solution[iColumn];
209        value = CoinMax(value, lower[iColumn]);
210        value = CoinMin(value, upper[iColumn]);
211        double nearest = floor(value + 0.5);
212#ifndef NDEBUG
213        double distance = fabs(value - nearest);
214        assert(distance <= integerTolerance);
215#endif
216        solver->setColLower(iColumn, nearest);
217        solver->setColUpper(iColumn, nearest);
218    }
219}
220// Redoes data when sequence numbers change
221void
222CbcClique::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
223{
224    model_ = model;
225    int n2 = 0;
226    for (int j = 0; j < numberMembers_; j++) {
227        int iColumn = members_[j];
228        int i;
229        for (i = 0; i < numberColumns; i++) {
230            if (originalColumns[i] == iColumn)
231                break;
232        }
233        if (i < numberColumns) {
234            members_[n2] = i;
235            type_[n2++] = type_[j];
236        }
237    }
238    if (n2 < numberMembers_) {
239        //printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
240        numberMembers_ = n2;
241    }
242    // Find out how many non sos
243    int i;
244    numberNonSOSMembers_ = 0;
245    for (i = 0; i < numberMembers_; i++)
246        if (!type_[i])
247            numberNonSOSMembers_++;
248}
249CbcBranchingObject *
250CbcClique::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way)
251{
252    int numberUnsatis = 0;
253    int j;
254    int nUp = 0;
255    int nDown = 0;
256    int numberFree = numberMembers_;
257    const int * integer = model_->integerVariable();
258    //OsiSolverInterface * solver = model_->solver();
259    const double * solution = model_->testSolution();
260    const double * lower = solver->getColLower();
261    const double * upper = solver->getColUpper();
262    int * upList = new int[numberMembers_];
263    int * downList = new int[numberMembers_];
264    double * sort = new double[numberMembers_];
265    double integerTolerance =
266        model_->getDblParam(CbcModel::CbcIntegerTolerance);
267
268    double slackValue = 0.0;
269    for (j = 0; j < numberMembers_; j++) {
270        int sequence = members_[j];
271        int iColumn = integer[sequence];
272        double value = solution[iColumn];
273        value = CoinMax(value, lower[iColumn]);
274        value = CoinMin(value, upper[iColumn]);
275        double nearest = floor(value + 0.5);
276        double distance = fabs(value - nearest);
277        if (distance > integerTolerance) {
278            if (!type_[j])
279                value = 1.0 - value; // non SOS
280            // if slack then choose that
281            if (j == slack_ && value > 0.05)
282                slackValue = value;
283            value = -value; // for sort
284            upList[numberUnsatis] = j;
285            sort[numberUnsatis++] = value;
286        } else if (upper[iColumn] > lower[iColumn]) {
287            upList[--numberFree] = j;
288        }
289    }
290    assert (numberUnsatis);
291    if (!slackValue) {
292        // sort
293        CoinSort_2(sort, sort + numberUnsatis, upList);
294        // put first in up etc
295        int kWay = 1;
296        for (j = 0; j < numberUnsatis; j++) {
297            if (kWay > 0)
298                upList[nUp++] = upList[j];
299            else
300                downList[nDown++] = upList[j];
301            kWay = -kWay;
302        }
303        for (j = numberFree; j < numberMembers_; j++) {
304            if (kWay > 0)
305                upList[nUp++] = upList[j];
306            else
307                downList[nDown++] = upList[j];
308            kWay = -kWay;
309        }
310    } else {
311        // put slack to 0 in first way
312        nUp = 1;
313        upList[0] = slack_;
314        for (j = 0; j < numberUnsatis; j++) {
315            downList[nDown++] = upList[j];
316        }
317        for (j = numberFree; j < numberMembers_; j++) {
318            downList[nDown++] = upList[j];
319        }
320    }
321    // create object
322    CbcBranchingObject * branch;
323    if (numberMembers_ <= 64)
324        branch = new CbcCliqueBranchingObject(model_, this, way,
325                                              nDown, downList, nUp, upList);
326    else
327        branch = new CbcLongCliqueBranchingObject(model_, this, way,
328                nDown, downList, nUp, upList);
329    delete [] upList;
330    delete [] downList;
331    delete [] sort;
332    return branch;
333}
334
335//##############################################################################
336
337// Default Constructor
338CbcSOS::CbcSOS ()
339        : CbcObject(),
340        members_(NULL),
341        weights_(NULL),
342        shadowEstimateDown_(1.0),
343        shadowEstimateUp_(1.0),
344        downDynamicPseudoRatio_(0.0),
345        upDynamicPseudoRatio_(0.0),
346        numberTimesDown_(0),
347        numberTimesUp_(0),
348        numberMembers_(0),
349        sosType_(-1),
350        integerValued_(false)
351{
352}
353
354// Useful constructor (which are indices)
355CbcSOS::CbcSOS (CbcModel * model,  int numberMembers,
356                const int * which, const double * weights, int identifier, int type)
357        : CbcObject(model),
358        shadowEstimateDown_(1.0),
359        shadowEstimateUp_(1.0),
360        downDynamicPseudoRatio_(0.0),
361        upDynamicPseudoRatio_(0.0),
362        numberTimesDown_(0),
363        numberTimesUp_(0),
364        numberMembers_(numberMembers),
365        sosType_(type)
366{
367    id_ = identifier;
368    integerValued_ = type == 1;
369    if (integerValued_) {
370        // check all members integer
371        OsiSolverInterface * solver = model->solver();
372        if (solver) {
373            for (int i = 0; i < numberMembers_; i++) {
374                if (!solver->isInteger(which[i]))
375                    integerValued_ = false;
376            }
377        } else {
378            // can't tell
379            integerValued_ = false;
380        }
381    }
382    if (numberMembers_) {
383        members_ = new int[numberMembers_];
384        weights_ = new double[numberMembers_];
385        memcpy(members_, which, numberMembers_*sizeof(int));
386        if (weights) {
387            memcpy(weights_, weights, numberMembers_*sizeof(double));
388        } else {
389            for (int i = 0; i < numberMembers_; i++)
390                weights_[i] = i;
391        }
392        // sort so weights increasing
393        CoinSort_2(weights_, weights_ + numberMembers_, members_);
394        double last = -COIN_DBL_MAX;
395        int i;
396        for (i = 0; i < numberMembers_; i++) {
397            double possible = CoinMax(last + 1.0e-10, weights_[i]);
398            weights_[i] = possible;
399            last = possible;
400        }
401    } else {
402        members_ = NULL;
403        weights_ = NULL;
404    }
405    assert (sosType_ > 0 && sosType_ < 3);
406}
407
408// Copy constructor
409CbcSOS::CbcSOS ( const CbcSOS & rhs)
410        : CbcObject(rhs)
411{
412    shadowEstimateDown_ = rhs.shadowEstimateDown_;
413    shadowEstimateUp_ = rhs.shadowEstimateUp_;
414    downDynamicPseudoRatio_ = rhs.downDynamicPseudoRatio_;
415    upDynamicPseudoRatio_ = rhs.upDynamicPseudoRatio_;
416    numberTimesDown_ = rhs.numberTimesDown_;
417    numberTimesUp_ = rhs.numberTimesUp_;
418    numberMembers_ = rhs.numberMembers_;
419    sosType_ = rhs.sosType_;
420    integerValued_ = rhs.integerValued_;
421    if (numberMembers_) {
422        members_ = new int[numberMembers_];
423        weights_ = new double[numberMembers_];
424        memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
425        memcpy(weights_, rhs.weights_, numberMembers_*sizeof(double));
426    } else {
427        members_ = NULL;
428        weights_ = NULL;
429    }
430}
431
432// Clone
433CbcObject *
434CbcSOS::clone() const
435{
436    return new CbcSOS(*this);
437}
438
439// Assignment operator
440CbcSOS &
441CbcSOS::operator=( const CbcSOS & rhs)
442{
443    if (this != &rhs) {
444        CbcObject::operator=(rhs);
445        delete [] members_;
446        delete [] weights_;
447        shadowEstimateDown_ = rhs.shadowEstimateDown_;
448        shadowEstimateUp_ = rhs.shadowEstimateUp_;
449        downDynamicPseudoRatio_ = rhs.downDynamicPseudoRatio_;
450        upDynamicPseudoRatio_ = rhs.upDynamicPseudoRatio_;
451        numberTimesDown_ = rhs.numberTimesDown_;
452        numberTimesUp_ = rhs.numberTimesUp_;
453        numberMembers_ = rhs.numberMembers_;
454        sosType_ = rhs.sosType_;
455        integerValued_ = rhs.integerValued_;
456        if (numberMembers_) {
457            members_ = new int[numberMembers_];
458            weights_ = new double[numberMembers_];
459            memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
460            memcpy(weights_, rhs.weights_, numberMembers_*sizeof(double));
461        } else {
462            members_ = NULL;
463            weights_ = NULL;
464        }
465    }
466    return *this;
467}
468
469// Destructor
470CbcSOS::~CbcSOS ()
471{
472    delete [] members_;
473    delete [] weights_;
474}
475double
476CbcSOS::infeasibility(const OsiBranchingInformation * info,
477                      int &preferredWay) const
478{
479    int j;
480    int firstNonZero = -1;
481    int lastNonZero = -1;
482    OsiSolverInterface * solver = model_->solver();
483    const double * solution = model_->testSolution();
484    //const double * lower = solver->getColLower();
485    const double * upper = solver->getColUpper();
486    //double largestValue=0.0;
487    double integerTolerance =
488        model_->getDblParam(CbcModel::CbcIntegerTolerance);
489    double weight = 0.0;
490    double sum = 0.0;
491
492    // check bounds etc
493    double lastWeight = -1.0e100;
494    for (j = 0; j < numberMembers_; j++) {
495        int iColumn = members_[j];
496        if (lastWeight >= weights_[j] - 1.0e-7)
497            throw CoinError("Weights too close together in SOS", "infeasibility", "CbcSOS");
498        double value = CoinMax(0.0, solution[iColumn]);
499        sum += value;
500        if (value > integerTolerance && upper[iColumn]) {
501            // Possibly due to scaling a fixed variable might slip through
502            if (value > upper[iColumn]) {
503                value = upper[iColumn];
504                // Could change to #ifdef CBC_DEBUG
505#ifndef NDEBUG
506                if (model_->messageHandler()->logLevel() > 2)
507                    printf("** Variable %d (%d) has value %g and upper bound of %g\n",
508                           iColumn, j, value, upper[iColumn]);
509#endif
510            }
511            weight += weights_[j] * value;
512            if (firstNonZero < 0)
513                firstNonZero = j;
514            lastNonZero = j;
515        }
516    }
517    preferredWay = 1;
518    if (lastNonZero - firstNonZero >= sosType_) {
519        // find where to branch
520        assert (sum > 0.0);
521        weight /= sum;
522        if (info->defaultDual_ >= 0.0 && info->usefulRegion_ && info->columnStart_) {
523            assert (sosType_ == 1);
524            int iWhere;
525            for (iWhere = firstNonZero; iWhere < lastNonZero - 1; iWhere++) {
526                if (weight < weights_[iWhere+1]) {
527                    break;
528                }
529            }
530            int jColumnDown = members_[iWhere];
531            int jColumnUp = members_[iWhere+1];
532            int n = 0;
533            CoinBigIndex j;
534            double objMove = info->objective_[jColumnDown];
535            for (j = info->columnStart_[jColumnDown];
536                    j < info->columnStart_[jColumnDown] + info->columnLength_[jColumnDown]; j++) {
537                double value = info->elementByColumn_[j];
538                int iRow = info->row_[j];
539                info->indexRegion_[n++] = iRow;
540                info->usefulRegion_[iRow] = value;
541            }
542            for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++) {
543                int jColumn = members_[iWhere];
544                double solValue = info->solution_[jColumn];
545                if (!solValue)
546                    continue;
547                objMove -= info->objective_[jColumn] * solValue;
548                for (j = info->columnStart_[jColumn];
549                        j < info->columnStart_[jColumn] + info->columnLength_[jColumn]; j++) {
550                    double value = -info->elementByColumn_[j] * solValue;
551                    int iRow = info->row_[j];
552                    double oldValue = info->usefulRegion_[iRow];
553                    if (!oldValue) {
554                        info->indexRegion_[n++] = iRow;
555                    } else {
556                        value += oldValue;
557                        if (!value)
558                            value = 1.0e-100;
559                    }
560                    info->usefulRegion_[iRow] = value;
561                }
562            }
563            const double * pi = info->pi_;
564            const double * activity = info->rowActivity_;
565            const double * lower = info->rowLower_;
566            const double * upper = info->rowUpper_;
567            double tolerance = info->primalTolerance_;
568            double direction = info->direction_;
569            shadowEstimateDown_ = objMove * direction;
570            bool infeasible = false;
571            for (int k = 0; k < n; k++) {
572                int iRow = info->indexRegion_[k];
573                double movement = info->usefulRegion_[iRow];
574                // not this time info->usefulRegion_[iRow]=0.0;
575                if (lower[iRow] < -1.0e20)
576                    assert (pi[iRow] <= 1.0e-3);
577                if (upper[iRow] > 1.0e20)
578                    assert (pi[iRow] >= -1.0e-3);
579                double valueP = pi[iRow] * direction;
580                // if move makes infeasible then make at least default
581                double newValue = activity[iRow] + movement;
582                if (newValue > upper[iRow] + tolerance || newValue < lower[iRow] - tolerance) {
583                    shadowEstimateDown_ += fabs(movement) * CoinMax(fabs(valueP), info->defaultDual_);
584                    infeasible = true;
585                }
586            }
587            if (shadowEstimateDown_ < info->integerTolerance_) {
588                if (!infeasible) {
589                    shadowEstimateDown_ = 1.0e-10;
590#ifdef COIN_DEVELOP
591                    printf("zero pseudoShadowPrice\n");
592#endif
593                } else
594                    shadowEstimateDown_ = info->integerTolerance_;
595            }
596            // And other way
597            // take off
598            objMove -= info->objective_[jColumnDown];
599            for (j = info->columnStart_[jColumnDown];
600                    j < info->columnStart_[jColumnDown] + info->columnLength_[jColumnDown]; j++) {
601                double value = -info->elementByColumn_[j];
602                int iRow = info->row_[j];
603                double oldValue = info->usefulRegion_[iRow];
604                if (!oldValue) {
605                    info->indexRegion_[n++] = iRow;
606                } else {
607                    value += oldValue;
608                    if (!value)
609                        value = 1.0e-100;
610                }
611                info->usefulRegion_[iRow] = value;
612            }
613            // add on
614            objMove += info->objective_[jColumnUp];
615            for (j = info->columnStart_[jColumnUp];
616                    j < info->columnStart_[jColumnUp] + info->columnLength_[jColumnUp]; j++) {
617                double value = info->elementByColumn_[j];
618                int iRow = info->row_[j];
619                double oldValue = info->usefulRegion_[iRow];
620                if (!oldValue) {
621                    info->indexRegion_[n++] = iRow;
622                } else {
623                    value += oldValue;
624                    if (!value)
625                        value = 1.0e-100;
626                }
627                info->usefulRegion_[iRow] = value;
628            }
629            shadowEstimateUp_ = objMove * direction;
630            infeasible = false;
631            for (int k = 0; k < n; k++) {
632                int iRow = info->indexRegion_[k];
633                double movement = info->usefulRegion_[iRow];
634                info->usefulRegion_[iRow] = 0.0;
635                if (lower[iRow] < -1.0e20)
636                    assert (pi[iRow] <= 1.0e-3);
637                if (upper[iRow] > 1.0e20)
638                    assert (pi[iRow] >= -1.0e-3);
639                double valueP = pi[iRow] * direction;
640                // if move makes infeasible then make at least default
641                double newValue = activity[iRow] + movement;
642                if (newValue > upper[iRow] + tolerance || newValue < lower[iRow] - tolerance) {
643                    shadowEstimateUp_ += fabs(movement) * CoinMax(fabs(valueP), info->defaultDual_);
644                    infeasible = true;
645                }
646            }
647            if (shadowEstimateUp_ < info->integerTolerance_) {
648                if (!infeasible) {
649                    shadowEstimateUp_ = 1.0e-10;
650#ifdef COIN_DEVELOP
651                    printf("zero pseudoShadowPrice\n");
652#endif
653                } else
654                    shadowEstimateUp_ = info->integerTolerance_;
655            }
656            // adjust
657            double downCost = shadowEstimateDown_;
658            double upCost = shadowEstimateUp_;
659            if (numberTimesDown_)
660                downCost *= downDynamicPseudoRatio_ /
661                            static_cast<double> (numberTimesDown_);
662            if (numberTimesUp_)
663                upCost *= upDynamicPseudoRatio_ /
664                          static_cast<double> (numberTimesUp_);
665#define WEIGHT_AFTER 0.7
666#define WEIGHT_BEFORE 0.1
667            int stateOfSearch = model_->stateOfSearch() % 10;
668            double returnValue = 0.0;
669            double minValue = CoinMin(downCost, upCost);
670            double maxValue = CoinMax(downCost, upCost);
671            if (stateOfSearch <= 2) {
672                // no branching solution
673                returnValue = WEIGHT_BEFORE * minValue + (1.0 - WEIGHT_BEFORE) * maxValue;
674            } else {
675                returnValue = WEIGHT_AFTER * minValue + (1.0 - WEIGHT_AFTER) * maxValue;
676            }
677#ifdef PRINT_SHADOW
678            printf("%d id - down %d %g up %d %g shadow %g, %g returned %g\n",
679                   id_, numberTimesDown_, downDynamicPseudoRatio_,
680                   numberTimesUp_, upDynamicPseudoRatio_, shadowEstimateDown_,
681                   shadowEstimateUp_, returnValue);
682#endif
683            return returnValue;
684        } else {
685            double value = lastNonZero - firstNonZero + 1;
686            value *= 0.5 / static_cast<double> (numberMembers_);
687            return value;
688        }
689    } else {
690        return 0.0; // satisfied
691    }
692}
693
694// This looks at solution and sets bounds to contain solution
695void
696CbcSOS::feasibleRegion()
697{
698    int j;
699    int firstNonZero = -1;
700    int lastNonZero = -1;
701    OsiSolverInterface * solver = model_->solver();
702    const double * solution = model_->testSolution();
703    //const double * lower = solver->getColLower();
704    const double * upper = solver->getColUpper();
705    double integerTolerance =
706        model_->getDblParam(CbcModel::CbcIntegerTolerance);
707    double weight = 0.0;
708    double sum = 0.0;
709
710    for (j = 0; j < numberMembers_; j++) {
711        int iColumn = members_[j];
712        double value = CoinMax(0.0, solution[iColumn]);
713        sum += value;
714        if (value > integerTolerance && upper[iColumn]) {
715            weight += weights_[j] * value;
716            if (firstNonZero < 0)
717                firstNonZero = j;
718            lastNonZero = j;
719        }
720    }
721    assert (lastNonZero - firstNonZero < sosType_) ;
722    for (j = 0; j < firstNonZero; j++) {
723        int iColumn = members_[j];
724        solver->setColUpper(iColumn, 0.0);
725    }
726    for (j = lastNonZero + 1; j < numberMembers_; j++) {
727        int iColumn = members_[j];
728        solver->setColUpper(iColumn, 0.0);
729    }
730}
731// Redoes data when sequence numbers change
732void
733CbcSOS::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
734{
735    model_ = model;
736    int n2 = 0;
737    for (int j = 0; j < numberMembers_; j++) {
738        int iColumn = members_[j];
739        int i;
740        for (i = 0; i < numberColumns; i++) {
741            if (originalColumns[i] == iColumn)
742                break;
743        }
744        if (i < numberColumns) {
745            members_[n2] = i;
746            weights_[n2++] = weights_[j];
747        }
748    }
749    if (n2 < numberMembers_) {
750        //printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
751        numberMembers_ = n2;
752    }
753}
754CbcBranchingObject *
755CbcSOS::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way)
756{
757    int j;
758    const double * solution = model_->testSolution();
759    double integerTolerance =
760        model_->getDblParam(CbcModel::CbcIntegerTolerance);
761    //OsiSolverInterface * solver = model_->solver();
762    const double * upper = solver->getColUpper();
763    int firstNonFixed = -1;
764    int lastNonFixed = -1;
765    int firstNonZero = -1;
766    int lastNonZero = -1;
767    double weight = 0.0;
768    double sum = 0.0;
769    for (j = 0; j < numberMembers_; j++) {
770        int iColumn = members_[j];
771        if (upper[iColumn]) {
772            double value = CoinMax(0.0, solution[iColumn]);
773            sum += value;
774            if (firstNonFixed < 0)
775                firstNonFixed = j;
776            lastNonFixed = j;
777            if (value > integerTolerance) {
778                weight += weights_[j] * value;
779                if (firstNonZero < 0)
780                    firstNonZero = j;
781                lastNonZero = j;
782            }
783        }
784    }
785    assert (lastNonZero - firstNonZero >= sosType_) ;
786    // find where to branch
787    assert (sum > 0.0);
788    weight /= sum;
789    int iWhere;
790    double separator = 0.0;
791    for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
792        if (weight < weights_[iWhere+1])
793            break;
794    if (sosType_ == 1) {
795        // SOS 1
796        separator = 0.5 * (weights_[iWhere] + weights_[iWhere+1]);
797    } else {
798        // SOS 2
799        if (iWhere == firstNonFixed)
800            iWhere++;;
801        if (iWhere == lastNonFixed - 1)
802            iWhere = lastNonFixed - 2;
803        separator = weights_[iWhere+1];
804    }
805    // create object
806    CbcBranchingObject * branch;
807    branch = new CbcSOSBranchingObject(model_, this, way, separator);
808    branch->setOriginalObject(this);
809    return branch;
810}
811/* Pass in information on branch just done and create CbcObjectUpdateData instance.
812   If object does not need data then backward pointer will be NULL.
813   Assumes can get information from solver */
814CbcObjectUpdateData
815CbcSOS::createUpdateInformation(const OsiSolverInterface * solver,
816                                const CbcNode * node,
817                                const CbcBranchingObject * branchingObject)
818{
819    double originalValue = node->objectiveValue();
820    int originalUnsatisfied = node->numberUnsatisfied();
821    double objectiveValue = solver->getObjValue() * solver->getObjSense();
822    int unsatisfied = 0;
823    int i;
824    //might be base model - doesn't matter
825    int numberIntegers = model_->numberIntegers();;
826    const double * solution = solver->getColSolution();
827    //const double * lower = solver->getColLower();
828    //const double * upper = solver->getColUpper();
829    double change = CoinMax(0.0, objectiveValue - originalValue);
830    int iStatus;
831    if (solver->isProvenOptimal())
832        iStatus = 0; // optimal
833    else if (solver->isIterationLimitReached()
834             && !solver->isDualObjectiveLimitReached())
835        iStatus = 2; // unknown
836    else
837        iStatus = 1; // infeasible
838
839    bool feasible = iStatus != 1;
840    if (feasible) {
841        double integerTolerance =
842            model_->getDblParam(CbcModel::CbcIntegerTolerance);
843        const int * integerVariable = model_->integerVariable();
844        for (i = 0; i < numberIntegers; i++) {
845            int j = integerVariable[i];
846            double value = solution[j];
847            double nearest = floor(value + 0.5);
848            if (fabs(value - nearest) > integerTolerance)
849                unsatisfied++;
850        }
851    }
852    int way = branchingObject->way();
853    way = - way; // because after branch so moved on
854    double value = branchingObject->value();
855    CbcObjectUpdateData newData (this, way,
856                                 change, iStatus,
857                                 originalUnsatisfied - unsatisfied, value);
858    newData.originalObjective_ = originalValue;
859    // Solvers know about direction
860    double direction = solver->getObjSense();
861    solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
862    newData.cutoff_ *= direction;
863    return newData;
864}
865// Update object by CbcObjectUpdateData
866void
867CbcSOS::updateInformation(const CbcObjectUpdateData & data)
868{
869    bool feasible = data.status_ != 1;
870    int way = data.way_;
871    //double value = data.branchingValue_;
872    double originalValue = data.originalObjective_;
873    double change = data.change_;
874    if (way < 0) {
875        // down
876        if (!feasible) {
877            double distanceToCutoff = 0.0;
878            //double objectiveValue = model_->getCurrentMinimizationObjValue();
879            distanceToCutoff =  model_->getCutoff()  - originalValue;
880            if (distanceToCutoff < 1.0e20)
881                change = distanceToCutoff * 2.0;
882            else
883                change = (downDynamicPseudoRatio_ * shadowEstimateDown_ + 1.0e-3) * 10.0;
884        }
885        change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
886#ifdef PRINT_SHADOW
887        if (numberTimesDown_)
888            printf("Updating id %d - down change %g (true %g) - ndown %d estimated change %g - raw shadow estimate %g\n",
889                   id_, change, data.change_, numberTimesDown_, shadowEstimateDown_*
890                   (downDynamicPseudoRatio_ / ((double) numberTimesDown_)),
891                   shadowEstimateDown_);
892        else
893            printf("Updating id %d - down change %g (true %g) - shadow estimate %g\n",
894                   id_, change, data.change_, shadowEstimateDown_);
895#endif
896        numberTimesDown_++;
897        downDynamicPseudoRatio_ += change / shadowEstimateDown_;
898    } else {
899        // up
900        if (!feasible) {
901            double distanceToCutoff = 0.0;
902            //double objectiveValue = model_->getCurrentMinimizationObjValue();
903            distanceToCutoff =  model_->getCutoff()  - originalValue;
904            if (distanceToCutoff < 1.0e20)
905                change = distanceToCutoff * 2.0;
906            else
907                change = (upDynamicPseudoRatio_ * shadowEstimateUp_ + 1.0e-3) * 10.0;
908        }
909        change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
910#ifdef PRINT_SHADOW
911        if (numberTimesUp_)
912            printf("Updating id %d - up change %g (true %g) - nup %d estimated change %g - raw shadow estimate %g\n",
913                   id_, change, data.change_, numberTimesUp_, shadowEstimateUp_*
914                   (upDynamicPseudoRatio_ / ((double) numberTimesUp_)),
915                   shadowEstimateUp_);
916        else
917            printf("Updating id %d - up change %g (true %g) - shadow estimate %g\n",
918                   id_, change, data.change_, shadowEstimateUp_);
919#endif
920        numberTimesUp_++;
921        upDynamicPseudoRatio_ += change / shadowEstimateUp_;
922    }
923}
924
925/* Create an OsiSolverBranch object
926
927This returns NULL if branch not represented by bound changes
928*/
929OsiSolverBranch *
930CbcSOS::solverBranch() const
931{
932    int j;
933    const double * solution = model_->testSolution();
934    double integerTolerance =
935        model_->getDblParam(CbcModel::CbcIntegerTolerance);
936    OsiSolverInterface * solver = model_->solver();
937    const double * upper = solver->getColUpper();
938    int firstNonFixed = -1;
939    int lastNonFixed = -1;
940    int firstNonZero = -1;
941    int lastNonZero = -1;
942    double weight = 0.0;
943    double sum = 0.0;
944    double * fix = new double[numberMembers_];
945    int * which = new int[numberMembers_];
946    for (j = 0; j < numberMembers_; j++) {
947        int iColumn = members_[j];
948        // fix all on one side or other (even if fixed)
949        fix[j] = 0.0;
950        which[j] = iColumn;
951        if (upper[iColumn]) {
952            double value = CoinMax(0.0, solution[iColumn]);
953            sum += value;
954            if (firstNonFixed < 0)
955                firstNonFixed = j;
956            lastNonFixed = j;
957            if (value > integerTolerance) {
958                weight += weights_[j] * value;
959                if (firstNonZero < 0)
960                    firstNonZero = j;
961                lastNonZero = j;
962            }
963        }
964    }
965    assert (lastNonZero - firstNonZero >= sosType_) ;
966    // find where to branch
967    assert (sum > 0.0);
968    weight /= sum;
969    // down branch fixes ones above weight to 0
970    int iWhere;
971    int iDownStart = 0;
972    int iUpEnd = 0;
973    for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
974        if (weight < weights_[iWhere+1])
975            break;
976    if (sosType_ == 1) {
977        // SOS 1
978        iUpEnd = iWhere + 1;
979        iDownStart = iUpEnd;
980    } else {
981        // SOS 2
982        if (iWhere == firstNonFixed)
983            iWhere++;;
984        if (iWhere == lastNonFixed - 1)
985            iWhere = lastNonFixed - 2;
986        iUpEnd = iWhere + 1;
987        iDownStart = iUpEnd + 1;
988    }
989    //
990    OsiSolverBranch * branch = new OsiSolverBranch();
991    branch->addBranch(-1, 0, NULL, NULL, numberMembers_ - iDownStart, which + iDownStart, fix);
992    branch->addBranch(1, 0, NULL, NULL, iUpEnd, which, fix);
993    delete [] fix;
994    delete [] which;
995    return branch;
996}
997// Construct an OsiSOS object
998OsiSOS *
999CbcSOS::osiObject(const OsiSolverInterface * solver) const
1000{
1001    OsiSOS * obj = new OsiSOS(solver, numberMembers_, members_, weights_, sosType_);
1002    obj->setPriority(priority());
1003    return obj;
1004}
1005
1006//##############################################################################
1007
1008/** Default Constructor
1009
1010  Equivalent to an unspecified binary variable.
1011*/
1012CbcSimpleInteger::CbcSimpleInteger ()
1013        : CbcObject(),
1014        originalLower_(0.0),
1015        originalUpper_(1.0),
1016        breakEven_(0.5),
1017        columnNumber_(-1),
1018        preferredWay_(0)
1019{
1020}
1021
1022/** Useful constructor
1023
1024  Loads actual upper & lower bounds for the specified variable.
1025*/
1026CbcSimpleInteger::CbcSimpleInteger ( CbcModel * model, int iColumn, double breakEven)
1027        : CbcObject(model)
1028{
1029    columnNumber_ = iColumn ;
1030    originalLower_ = model->solver()->getColLower()[columnNumber_] ;
1031    originalUpper_ = model->solver()->getColUpper()[columnNumber_] ;
1032    breakEven_ = breakEven;
1033    assert (breakEven_ > 0.0 && breakEven_ < 1.0);
1034    preferredWay_ = 0;
1035}
1036
1037
1038// Copy constructor
1039CbcSimpleInteger::CbcSimpleInteger ( const CbcSimpleInteger & rhs)
1040        : CbcObject(rhs)
1041
1042{
1043    columnNumber_ = rhs.columnNumber_;
1044    originalLower_ = rhs.originalLower_;
1045    originalUpper_ = rhs.originalUpper_;
1046    breakEven_ = rhs.breakEven_;
1047    preferredWay_ = rhs.preferredWay_;
1048}
1049
1050// Clone
1051CbcObject *
1052CbcSimpleInteger::clone() const
1053{
1054    return new CbcSimpleInteger(*this);
1055}
1056
1057// Assignment operator
1058CbcSimpleInteger &
1059CbcSimpleInteger::operator=( const CbcSimpleInteger & rhs)
1060{
1061    if (this != &rhs) {
1062        CbcObject::operator=(rhs);
1063        columnNumber_ = rhs.columnNumber_;
1064        originalLower_ = rhs.originalLower_;
1065        originalUpper_ = rhs.originalUpper_;
1066        breakEven_ = rhs.breakEven_;
1067        preferredWay_ = rhs.preferredWay_;
1068    }
1069    return *this;
1070}
1071
1072// Destructor
1073CbcSimpleInteger::~CbcSimpleInteger ()
1074{
1075}
1076// Construct an OsiSimpleInteger object
1077OsiSimpleInteger *
1078CbcSimpleInteger::osiObject() const
1079{
1080    OsiSimpleInteger * obj = new OsiSimpleInteger(columnNumber_,
1081            originalLower_, originalUpper_);
1082    obj->setPriority(priority());
1083    return obj;
1084}
1085double
1086CbcSimpleInteger::infeasibility(const OsiBranchingInformation * info,
1087                                int &preferredWay) const
1088{
1089    double value = info->solution_[columnNumber_];
1090    value = CoinMax(value, info->lower_[columnNumber_]);
1091    value = CoinMin(value, info->upper_[columnNumber_]);
1092    double nearest = floor(value + (1.0 - breakEven_));
1093    assert (breakEven_ > 0.0 && breakEven_ < 1.0);
1094    if (nearest > value)
1095        preferredWay = 1;
1096    else
1097        preferredWay = -1;
1098    if (preferredWay_)
1099        preferredWay = preferredWay_;
1100    double weight = fabs(value - nearest);
1101    // normalize so weight is 0.5 at break even
1102    if (nearest < value)
1103        weight = (0.5 / breakEven_) * weight;
1104    else
1105        weight = (0.5 / (1.0 - breakEven_)) * weight;
1106    if (fabs(value - nearest) <= info->integerTolerance_)
1107        return 0.0;
1108    else
1109        return weight;
1110}
1111double
1112CbcSimpleInteger::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
1113{
1114    double value = info->solution_[columnNumber_];
1115#ifdef COIN_DEVELOP
1116    if (fabs(value - floor(value + 0.5)) > 1.0e-5)
1117        printf("value for %d away from integer %g\n", columnNumber_, value);
1118#endif
1119    double newValue = CoinMax(value, info->lower_[columnNumber_]);
1120    newValue = CoinMin(newValue, info->upper_[columnNumber_]);
1121    newValue = floor(newValue + 0.5);
1122    solver->setColLower(columnNumber_, newValue);
1123    solver->setColUpper(columnNumber_, newValue);
1124    return fabs(value - newValue);
1125}
1126
1127/* Create an OsiSolverBranch object
1128
1129This returns NULL if branch not represented by bound changes
1130*/
1131OsiSolverBranch *
1132CbcSimpleInteger::solverBranch(OsiSolverInterface * /*solver*/,
1133                               const OsiBranchingInformation * info) const
1134{
1135    double value = info->solution_[columnNumber_];
1136    value = CoinMax(value, info->lower_[columnNumber_]);
1137    value = CoinMin(value, info->upper_[columnNumber_]);
1138    assert (info->upper_[columnNumber_] > info->lower_[columnNumber_]);
1139#ifndef NDEBUG
1140    double nearest = floor(value + 0.5);
1141    assert (fabs(value - nearest) > info->integerTolerance_);
1142#endif
1143    OsiSolverBranch * branch = new OsiSolverBranch();
1144    branch->addBranch(columnNumber_, value);
1145    return branch;
1146}
1147// Creates a branching object
1148CbcBranchingObject *
1149CbcSimpleInteger::createCbcBranch(OsiSolverInterface * /*solver*/,
1150                                  const OsiBranchingInformation * info, int way)
1151{
1152    CbcIntegerBranchingObject * branch = new CbcIntegerBranchingObject(model_, 0, -1, 0.5);
1153    fillCreateBranch(branch, info, way);
1154    return branch;
1155}
1156// Fills in a created branching object
1157void
1158CbcSimpleInteger::fillCreateBranch(CbcIntegerBranchingObject * branch, const OsiBranchingInformation * info, int way)
1159{
1160    branch->setOriginalObject(this);
1161    double value = info->solution_[columnNumber_];
1162    value = CoinMax(value, info->lower_[columnNumber_]);
1163    value = CoinMin(value, info->upper_[columnNumber_]);
1164    assert (info->upper_[columnNumber_] > info->lower_[columnNumber_]);
1165    if (!info->hotstartSolution_ && priority_ != -999) {
1166#ifndef NDEBUG
1167        double nearest = floor(value + 0.5);
1168        assert (fabs(value - nearest) > info->integerTolerance_);
1169#endif
1170    } else if (info->hotstartSolution_) {
1171        double targetValue = info->hotstartSolution_[columnNumber_];
1172        if (way > 0)
1173            value = targetValue - 0.1;
1174        else
1175            value = targetValue + 0.1;
1176    } else {
1177        if (value <= info->lower_[columnNumber_])
1178            value += 0.1;
1179        else if (value >= info->upper_[columnNumber_])
1180            value -= 0.1;
1181    }
1182    assert (value >= info->lower_[columnNumber_] &&
1183            value <= info->upper_[columnNumber_]);
1184    branch->fillPart(columnNumber_, way, value);
1185}
1186/* Column number if single column object -1 otherwise,
1187   so returns >= 0
1188   Used by heuristics
1189*/
1190int
1191CbcSimpleInteger::columnNumber() const
1192{
1193    return columnNumber_;
1194}
1195/* Reset variable bounds to their original values.
1196
1197    Bounds may be tightened, so it may be good to be able to set this info in object.
1198*/
1199void
1200CbcSimpleInteger::resetBounds(const OsiSolverInterface * solver)
1201{
1202    originalLower_ = solver->getColLower()[columnNumber_] ;
1203    originalUpper_ = solver->getColUpper()[columnNumber_] ;
1204}
1205
1206/*  Change column numbers after preprocessing
1207 */
1208void
1209CbcSimpleInteger::resetSequenceEtc(int /*numberColumns*/,
1210                                   const int * originalColumns)
1211{
1212    //assert (numberColumns>0);
1213    int iColumn;
1214#if 0
1215    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1216        if (columnNumber_ == originalColumns[iColumn])
1217            break;
1218    }
1219    assert (iColumn < numberColumns);
1220#else
1221    iColumn = originalColumns[columnNumber_];
1222    assert (iColumn >= 0);
1223#endif
1224    columnNumber_ = iColumn;
1225}
1226// This looks at solution and sets bounds to contain solution
1227/** More precisely: it first forces the variable within the existing
1228    bounds, and then tightens the bounds to fix the variable at the
1229    nearest integer value.
1230*/
1231void
1232CbcSimpleInteger::feasibleRegion()
1233{
1234    abort();
1235}
1236
1237//##############################################################################
1238
1239// Default Constructor
1240CbcIntegerBranchingObject::CbcIntegerBranchingObject()
1241        : CbcBranchingObject()
1242{
1243    down_[0] = 0.0;
1244    down_[1] = 0.0;
1245    up_[0] = 0.0;
1246    up_[1] = 0.0;
1247#ifdef FUNNY_BRANCHING
1248    variables_ = NULL;
1249    newBounds_ = NULL;
1250    numberExtraChangedBounds_ = 0;
1251#endif
1252}
1253// Useful constructor
1254CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model,
1255        int variable, int way , double value)
1256        : CbcBranchingObject(model, variable, way, value)
1257{
1258    int iColumn = variable;
1259    assert (model_->solver()->getNumCols() > 0);
1260    down_[0] = model_->solver()->getColLower()[iColumn];
1261    down_[1] = floor(value_);
1262    up_[0] = ceil(value_);
1263    up_[1] = model->getColUpper()[iColumn];
1264#ifdef FUNNY_BRANCHING
1265    variables_ = NULL;
1266    newBounds_ = NULL;
1267    numberExtraChangedBounds_ = 0;
1268#endif
1269}
1270// Does part of constructor
1271void
1272CbcIntegerBranchingObject::fillPart (int variable,
1273                                     int way , double value)
1274{
1275    //originalObject_=NULL;
1276    branchIndex_ = 0;
1277    value_ = value;
1278    numberBranches_ = 2;
1279    //model_= model;
1280    //originalCbcObject_=NULL;
1281    variable_ = variable;
1282    way_ = way;
1283    int iColumn = variable;
1284    down_[0] = model_->solver()->getColLower()[iColumn];
1285    down_[1] = floor(value_);
1286    up_[0] = ceil(value_);
1287    up_[1] = model_->getColUpper()[iColumn];
1288}
1289// Useful constructor for fixing
1290CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model,
1291        int variable, int way,
1292        double lowerValue,
1293        double upperValue)
1294        : CbcBranchingObject(model, variable, way, lowerValue)
1295{
1296    setNumberBranchesLeft(1);
1297    down_[0] = lowerValue;
1298    down_[1] = upperValue;
1299    up_[0] = lowerValue;
1300    up_[1] = upperValue;
1301#ifdef FUNNY_BRANCHING
1302    variables_ = NULL;
1303    newBounds_ = NULL;
1304    numberExtraChangedBounds_ = 0;
1305#endif
1306}
1307
1308
1309// Copy constructor
1310CbcIntegerBranchingObject::CbcIntegerBranchingObject ( const CbcIntegerBranchingObject & rhs) : CbcBranchingObject(rhs)
1311{
1312    down_[0] = rhs.down_[0];
1313    down_[1] = rhs.down_[1];
1314    up_[0] = rhs.up_[0];
1315    up_[1] = rhs.up_[1];
1316#ifdef FUNNY_BRANCHING
1317    numberExtraChangedBounds_ = rhs.numberExtraChangedBounds_;
1318    int size = numberExtraChangedBounds_ * (sizeof(double) + sizeof(int));
1319    char * temp = new char [size];
1320    newBounds_ = (double *) temp;
1321    variables_ = (int *) (newBounds_ + numberExtraChangedBounds_);
1322
1323    int i ;
1324    for (i = 0; i < numberExtraChangedBounds_; i++) {
1325        variables_[i] = rhs.variables_[i];
1326        newBounds_[i] = rhs.newBounds_[i];
1327    }
1328#endif
1329}
1330
1331// Assignment operator
1332CbcIntegerBranchingObject &
1333CbcIntegerBranchingObject::operator=( const CbcIntegerBranchingObject & rhs)
1334{
1335    if (this != &rhs) {
1336        CbcBranchingObject::operator=(rhs);
1337        down_[0] = rhs.down_[0];
1338        down_[1] = rhs.down_[1];
1339        up_[0] = rhs.up_[0];
1340        up_[1] = rhs.up_[1];
1341#ifdef FUNNY_BRANCHING
1342        delete [] newBounds_;
1343        numberExtraChangedBounds_ = rhs.numberExtraChangedBounds_;
1344        int size = numberExtraChangedBounds_ * (sizeof(double) + sizeof(int));
1345        char * temp = new char [size];
1346        newBounds_ = (double *) temp;
1347        variables_ = (int *) (newBounds_ + numberExtraChangedBounds_);
1348
1349        int i ;
1350        for (i = 0; i < numberExtraChangedBounds_; i++) {
1351            variables_[i] = rhs.variables_[i];
1352            newBounds_[i] = rhs.newBounds_[i];
1353        }
1354#endif
1355    }
1356    return *this;
1357}
1358CbcBranchingObject *
1359CbcIntegerBranchingObject::clone() const
1360{
1361    return (new CbcIntegerBranchingObject(*this));
1362}
1363
1364
1365// Destructor
1366CbcIntegerBranchingObject::~CbcIntegerBranchingObject ()
1367{
1368    // for debugging threads
1369    way_ = -23456789;
1370#ifdef FUNNY_BRANCHING
1371    delete [] newBounds_;
1372#endif
1373}
1374
1375/*
1376  Perform a branch by adjusting the bounds of the specified variable. Note
1377  that each arm of the branch advances the object to the next arm by
1378  advancing the value of way_.
1379
1380  Providing new values for the variable's lower and upper bounds for each
1381  branching direction gives a little bit of additional flexibility and will
1382  be easily extensible to multi-way branching.
1383  Returns change in guessed objective on next branch
1384*/
1385double
1386CbcIntegerBranchingObject::branch()
1387{
1388    // for debugging threads
1389    if (way_ < -1 || way_ > 100000) {
1390        printf("way %d, left %d, iCol %d, variable %d\n",
1391               way_, numberBranchesLeft(),
1392               originalCbcObject_->columnNumber(), variable_);
1393        assert (way_ != -23456789);
1394    }
1395    decrementNumberBranchesLeft();
1396    if (down_[1] == -COIN_DBL_MAX)
1397        return 0.0;
1398    int iColumn = originalCbcObject_->columnNumber();
1399    assert (variable_ == iColumn);
1400    double olb, oub ;
1401    olb = model_->solver()->getColLower()[iColumn] ;
1402    oub = model_->solver()->getColUpper()[iColumn] ;
1403    //#define TIGHTEN_BOUNDS
1404#ifndef TIGHTEN_BOUNDS
1405#ifdef COIN_DEVELOP
1406    if (olb != down_[0] || oub != up_[1]) {
1407        if (way_ > 0)
1408            printf("branching up on var %d: [%g,%g] => [%g,%g] - other [%g,%g]\n",
1409                   iColumn, olb, oub, up_[0], up_[1], down_[0], down_[1]) ;
1410        else
1411            printf("branching down on var %d: [%g,%g] => [%g,%g] - other [%g,%g]\n",
1412                   iColumn, olb, oub, down_[0], down_[1], up_[0], up_[1]) ;
1413    }
1414#endif
1415#endif
1416    if (way_ < 0) {
1417#ifdef CBC_DEBUG
1418        { double olb, oub ;
1419            olb = model_->solver()->getColLower()[iColumn] ;
1420            oub = model_->solver()->getColUpper()[iColumn] ;
1421            printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
1422                   iColumn, olb, oub, down_[0], down_[1]) ;
1423        }
1424#endif
1425#ifndef TIGHTEN_BOUNDS
1426        model_->solver()->setColLower(iColumn, down_[0]);
1427#else
1428        model_->solver()->setColLower(iColumn, CoinMax(down_[0], olb));
1429#endif
1430        model_->solver()->setColUpper(iColumn, down_[1]);
1431        //#define CBC_PRINT2
1432#ifdef CBC_PRINT2
1433        printf("%d branching down has bounds %g %g", iColumn, down_[0], down_[1]);
1434#endif
1435#ifdef FUNNY_BRANCHING
1436        // branch - do extra bounds
1437        for (int i = 0; i < numberExtraChangedBounds_; i++) {
1438            int variable = variables_[i];
1439            if ((variable&0x40000000) != 0) {
1440                // for going down
1441                int k = variable & 0x3fffffff;
1442                assert (k != iColumn);
1443                if ((variable&0x80000000) == 0) {
1444                    // lower bound changing
1445#ifdef CBC_PRINT2
1446                    printf(" extra for %d changes lower from %g to %g",
1447                           k, model_->solver()->getColLower()[k], newBounds_[i]);
1448#endif
1449                    model_->solver()->setColLower(k, newBounds_[i]);
1450                } else {
1451                    // upper bound changing
1452#ifdef CBC_PRINT2
1453                    printf(" extra for %d changes upper from %g to %g",
1454                           k, model_->solver()->getColUpper()[k], newBounds_[i]);
1455#endif
1456                    model_->solver()->setColUpper(k, newBounds_[i]);
1457                }
1458            }
1459        }
1460#endif
1461#ifdef CBC_PRINT2
1462        printf("\n");
1463#endif
1464        way_ = 1;
1465    } else {
1466#ifdef CBC_DEBUG
1467        { double olb, oub ;
1468            olb = model_->solver()->getColLower()[iColumn] ;
1469            oub = model_->solver()->getColUpper()[iColumn] ;
1470            printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
1471                   iColumn, olb, oub, up_[0], up_[1]) ;
1472        }
1473#endif
1474        model_->solver()->setColLower(iColumn, up_[0]);
1475#ifndef TIGHTEN_BOUNDS
1476        model_->solver()->setColUpper(iColumn, up_[1]);
1477#else
1478        model_->solver()->setColUpper(iColumn, CoinMin(up_[1], oub));
1479#endif
1480#ifdef CBC_PRINT2
1481        printf("%d branching up has bounds %g %g", iColumn, up_[0], up_[1]);
1482#endif
1483#ifdef FUNNY_BRANCHING
1484        // branch - do extra bounds
1485        for (int i = 0; i < numberExtraChangedBounds_; i++) {
1486            int variable = variables_[i];
1487            if ((variable&0x40000000) == 0) {
1488                // for going up
1489                int k = variable & 0x3fffffff;
1490                assert (k != iColumn);
1491                if ((variable&0x80000000) == 0) {
1492                    // lower bound changing
1493#ifdef CBC_PRINT2
1494                    printf(" extra for %d changes lower from %g to %g",
1495                           k, model_->solver()->getColLower()[k], newBounds_[i]);
1496#endif
1497                    model_->solver()->setColLower(k, newBounds_[i]);
1498                } else {
1499                    // upper bound changing
1500#ifdef CBC_PRINT2
1501                    printf(" extra for %d changes upper from %g to %g",
1502                           k, model_->solver()->getColUpper()[k], newBounds_[i]);
1503#endif
1504                    model_->solver()->setColUpper(k, newBounds_[i]);
1505                }
1506            }
1507        }
1508#endif
1509#ifdef CBC_PRINT2
1510        printf("\n");
1511#endif
1512        way_ = -1;        // Swap direction
1513    }
1514    double nlb = model_->solver()->getColLower()[iColumn];
1515    double nub = model_->solver()->getColUpper()[iColumn];
1516    if (nlb < olb) {
1517#ifndef NDEBUG
1518        printf("bad lb change for column %d from %g to %g\n", iColumn, olb, nlb);
1519#endif
1520        model_->solver()->setColLower(iColumn, CoinMin(olb, nub));
1521        nlb = olb;
1522    }
1523    if (nub > oub) {
1524#ifndef NDEBUG
1525        printf("bad ub change for column %d from %g to %g\n", iColumn, oub, nub);
1526#endif
1527        model_->solver()->setColUpper(iColumn, CoinMax(oub, nlb));
1528    }
1529#ifndef NDEBUG
1530    if (nlb < olb + 1.0e-8 && nub > oub - 1.0e-8 && false)
1531        printf("bad null change for column %d - bounds %g,%g\n", iColumn, olb, oub);
1532#endif
1533    return 0.0;
1534}
1535/* Update bounds in solver as in 'branch' and update given bounds.
1536   branchState is -1 for 'down' +1 for 'up' */
1537void
1538CbcIntegerBranchingObject::fix(OsiSolverInterface * /*solver*/,
1539                               double * lower, double * upper,
1540                               int branchState) const
1541{
1542    int iColumn = originalCbcObject_->columnNumber();
1543    assert (variable_ == iColumn);
1544    if (branchState < 0) {
1545        model_->solver()->setColLower(iColumn, down_[0]);
1546        lower[iColumn] = down_[0];
1547        model_->solver()->setColUpper(iColumn, down_[1]);
1548        upper[iColumn] = down_[1];
1549    } else {
1550        model_->solver()->setColLower(iColumn, up_[0]);
1551        lower[iColumn] = up_[0];
1552        model_->solver()->setColUpper(iColumn, up_[1]);
1553        upper[iColumn] = up_[1];
1554    }
1555}
1556#ifdef FUNNY_BRANCHING
1557// Deactivate bounds for branching
1558void
1559CbcIntegerBranchingObject::deactivate()
1560{
1561    down_[1] = -COIN_DBL_MAX;
1562}
1563int
1564CbcIntegerBranchingObject::applyExtraBounds(int iColumn, double lower, double upper, int way)
1565{
1566    // branch - do bounds
1567
1568    int i;
1569    int found = 0;
1570    if (variable_ == iColumn) {
1571        printf("odd applyExtra %d\n", iColumn);
1572        if (way < 0) {
1573            down_[0] = CoinMax(lower, down_[0]);
1574            down_[1] = CoinMin(upper, down_[1]);
1575            assert (down_[0] <= down_[1]);
1576        } else {
1577            up_[0] = CoinMax(lower, up_[0]);
1578            up_[1] = CoinMin(upper, up_[1]);
1579            assert (up_[0] <= up_[1]);
1580        }
1581        return 0;
1582    }
1583    int check = (way < 0) ? 0x40000000 : 0;
1584    double newLower = lower;
1585    double newUpper = upper;
1586    for (i = 0; i < numberExtraChangedBounds_; i++) {
1587        int variable = variables_[i];
1588        if ((variable&0x40000000) == check) {
1589            int k = variable & 0x3fffffff;
1590            if (k == iColumn) {
1591                if ((variable&0x80000000) == 0) {
1592                    // lower bound changing
1593                    found |= 1;
1594                    newBounds_[i] = CoinMax(lower, newBounds_[i]);
1595                    newLower = newBounds_[i];
1596                } else {
1597                    // upper bound changing
1598                    found |= 2;
1599                    newBounds_[i] = CoinMin(upper, newBounds_[i]);
1600                    newUpper = newBounds_[i];
1601                }
1602            }
1603        }
1604    }
1605    int nAdd = 0;
1606    if ((found&2) == 0) {
1607        // need to add new upper
1608        nAdd++;
1609    }
1610    if ((found&1) == 0) {
1611        // need to add new lower
1612        nAdd++;
1613    }
1614    if (nAdd) {
1615        int size = (numberExtraChangedBounds_ + nAdd) * (sizeof(double) + sizeof(int));
1616        char * temp = new char [size];
1617        double * newBounds = (double *) temp;
1618        int * variables = (int *) (newBounds + numberExtraChangedBounds_ + nAdd);
1619
1620        int i ;
1621        for (i = 0; i < numberExtraChangedBounds_; i++) {
1622            variables[i] = variables_[i];
1623            newBounds[i] = newBounds_[i];
1624        }
1625        delete [] newBounds_;
1626        newBounds_ = newBounds;
1627        variables_ = variables;
1628        if ((found&2) == 0) {
1629            // need to add new upper
1630            int variable = iColumn | 0x80000000;
1631            variables_[numberExtraChangedBounds_] = variable;
1632            newBounds_[numberExtraChangedBounds_++] = newUpper;
1633        }
1634        if ((found&1) == 0) {
1635            // need to add new lower
1636            int variable = iColumn;
1637            variables_[numberExtraChangedBounds_] = variable;
1638            newBounds_[numberExtraChangedBounds_++] = newLower;
1639        }
1640    }
1641
1642    return (newUpper >= newLower) ? 0 : 1;
1643}
1644#endif
1645// Print what would happen
1646void
1647CbcIntegerBranchingObject::print()
1648{
1649    int iColumn = originalCbcObject_->columnNumber();
1650    assert (variable_ == iColumn);
1651    if (way_ < 0) {
1652        {
1653            double olb, oub ;
1654            olb = model_->solver()->getColLower()[iColumn] ;
1655            oub = model_->solver()->getColUpper()[iColumn] ;
1656            printf("CbcInteger would branch down on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1657                   iColumn, variable_, olb, oub, down_[0], down_[1]) ;
1658        }
1659    } else {
1660        {
1661            double olb, oub ;
1662            olb = model_->solver()->getColLower()[iColumn] ;
1663            oub = model_->solver()->getColUpper()[iColumn] ;
1664            printf("CbcInteger would branch up on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1665                   iColumn, variable_, olb, oub, up_[0], up_[1]) ;
1666        }
1667    }
1668}
1669
1670/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1671    same type and must have the same original object, but they may have
1672    different feasible regions.
1673    Return the appropriate CbcRangeCompare value (first argument being the
1674    sub/superset if that's the case). In case of overlap (and if \c
1675    replaceIfOverlap is true) replace the current branching object with one
1676    whose feasible region is the overlap.
1677   */
1678CbcRangeCompare
1679CbcIntegerBranchingObject::compareBranchingObject
1680(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1681{
1682    const CbcIntegerBranchingObject* br =
1683        dynamic_cast<const CbcIntegerBranchingObject*>(brObj);
1684    assert(br);
1685    double* thisBd = way_ < 0 ? down_ : up_;
1686    const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
1687    return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
1688}
1689
1690//##############################################################################
1691
1692/** Default Constructor
1693
1694  Equivalent to an unspecified binary variable.
1695*/
1696CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ()
1697        : CbcSimpleInteger(),
1698        downPseudoCost_(1.0e-5),
1699        upPseudoCost_(1.0e-5),
1700        upDownSeparator_(-1.0),
1701        method_(0)
1702{
1703}
1704
1705/** Useful constructor
1706
1707  Loads actual upper & lower bounds for the specified variable.
1708*/
1709CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1710        int iColumn, double breakEven)
1711        : CbcSimpleInteger(model, iColumn, breakEven)
1712{
1713    const double * cost = model->getObjCoefficients();
1714    double costValue = CoinMax(1.0e-5, fabs(cost[iColumn]));
1715    // treat as if will cost what it says up
1716    upPseudoCost_ = costValue;
1717    // and balance at breakeven
1718    downPseudoCost_ = ((1.0 - breakEven_) * upPseudoCost_) / breakEven_;
1719    upDownSeparator_ = -1.0;
1720    method_ = 0;
1721}
1722
1723/** Useful constructor
1724
1725  Loads actual upper & lower bounds for the specified variable.
1726*/
1727CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1728        int iColumn, double downPseudoCost,
1729        double upPseudoCost)
1730        : CbcSimpleInteger(model, iColumn)
1731{
1732    downPseudoCost_ = CoinMax(1.0e-10, downPseudoCost);
1733    upPseudoCost_ = CoinMax(1.0e-10, upPseudoCost);
1734    breakEven_ = upPseudoCost_ / (upPseudoCost_ + downPseudoCost_);
1735    upDownSeparator_ = -1.0;
1736    method_ = 0;
1737}
1738// Useful constructor - passed and model index and pseudo costs
1739CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1740        int /*dummy*/,
1741        int iColumn,
1742        double downPseudoCost, double upPseudoCost)
1743{
1744    *this = CbcSimpleIntegerPseudoCost(model, iColumn, downPseudoCost, upPseudoCost);
1745    columnNumber_ = iColumn;
1746}
1747
1748// Copy constructor
1749CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ( const CbcSimpleIntegerPseudoCost & rhs)
1750        : CbcSimpleInteger(rhs),
1751        downPseudoCost_(rhs.downPseudoCost_),
1752        upPseudoCost_(rhs.upPseudoCost_),
1753        upDownSeparator_(rhs.upDownSeparator_),
1754        method_(rhs.method_)
1755
1756{
1757}
1758
1759// Clone
1760CbcObject *
1761CbcSimpleIntegerPseudoCost::clone() const
1762{
1763    return new CbcSimpleIntegerPseudoCost(*this);
1764}
1765
1766// Assignment operator
1767CbcSimpleIntegerPseudoCost &
1768CbcSimpleIntegerPseudoCost::operator=( const CbcSimpleIntegerPseudoCost & rhs)
1769{
1770    if (this != &rhs) {
1771        CbcSimpleInteger::operator=(rhs);
1772        downPseudoCost_ = rhs.downPseudoCost_;
1773        upPseudoCost_ = rhs.upPseudoCost_;
1774        upDownSeparator_ = rhs.upDownSeparator_;
1775        method_ = rhs.method_;
1776    }
1777    return *this;
1778}
1779
1780// Destructor
1781CbcSimpleIntegerPseudoCost::~CbcSimpleIntegerPseudoCost ()
1782{
1783}
1784CbcBranchingObject *
1785CbcSimpleIntegerPseudoCost::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way)
1786{
1787    //OsiSolverInterface * solver = model_->solver();
1788    const double * solution = model_->testSolution();
1789    const double * lower = solver->getColLower();
1790    const double * upper = solver->getColUpper();
1791    double value = solution[columnNumber_];
1792    value = CoinMax(value, lower[columnNumber_]);
1793    value = CoinMin(value, upper[columnNumber_]);
1794#ifndef NDEBUG
1795    double nearest = floor(value + 0.5);
1796    double integerTolerance =
1797        model_->getDblParam(CbcModel::CbcIntegerTolerance);
1798    assert (upper[columnNumber_] > lower[columnNumber_]);
1799#endif
1800    if (!model_->hotstartSolution()) {
1801        assert (fabs(value - nearest) > integerTolerance);
1802    } else {
1803        const double * hotstartSolution = model_->hotstartSolution();
1804        double targetValue = hotstartSolution[columnNumber_];
1805        if (way > 0)
1806            value = targetValue - 0.1;
1807        else
1808            value = targetValue + 0.1;
1809    }
1810    CbcIntegerPseudoCostBranchingObject * newObject =
1811        new CbcIntegerPseudoCostBranchingObject(model_, columnNumber_, way,
1812                                                value);
1813    double up =  upPseudoCost_ * (ceil(value) - value);
1814    double down =  downPseudoCost_ * (value - floor(value));
1815    double changeInGuessed = up - down;
1816    if (way > 0)
1817        changeInGuessed = - changeInGuessed;
1818    changeInGuessed = CoinMax(0.0, changeInGuessed);
1819    //if (way>0)
1820    //changeInGuessed += 1.0e8; // bias to stay up
1821    newObject->setChangeInGuessed(changeInGuessed);
1822    newObject->setOriginalObject(this);
1823    return newObject;
1824}
1825double
1826CbcSimpleIntegerPseudoCost::infeasibility(const OsiBranchingInformation * /*info*/,
1827        int &preferredWay) const
1828{
1829    OsiSolverInterface * solver = model_->solver();
1830    const double * solution = model_->testSolution();
1831    const double * lower = solver->getColLower();
1832    const double * upper = solver->getColUpper();
1833    if (upper[columnNumber_] == lower[columnNumber_]) {
1834        // fixed
1835        preferredWay = 1;
1836        return 0.0;
1837    }
1838    double value = solution[columnNumber_];
1839    value = CoinMax(value, lower[columnNumber_]);
1840    value = CoinMin(value, upper[columnNumber_]);
1841    /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
1842      solution[columnNumber_],upper[columnNumber_]);*/
1843    double nearest = floor(value + 0.5);
1844    double integerTolerance =
1845        model_->getDblParam(CbcModel::CbcIntegerTolerance);
1846    double below = floor(value + integerTolerance);
1847    double above = below + 1.0;
1848    if (above > upper[columnNumber_]) {
1849        above = below;
1850        below = above - 1;
1851    }
1852    double downCost = CoinMax((value - below) * downPseudoCost_, 0.0);
1853    double upCost = CoinMax((above - value) * upPseudoCost_, 0.0);
1854    if (downCost >= upCost)
1855        preferredWay = 1;
1856    else
1857        preferredWay = -1;
1858    // See if up down choice set
1859    if (upDownSeparator_ > 0.0) {
1860        preferredWay = (value - below >= upDownSeparator_) ? 1 : -1;
1861    }
1862    if (preferredWay_)
1863        preferredWay = preferredWay_;
1864    if (fabs(value - nearest) <= integerTolerance) {
1865        return 0.0;
1866    } else {
1867        // can't get at model so 1,2 don't make sense
1868        assert(method_ < 1 || method_ > 2);
1869        if (!method_)
1870            return CoinMin(downCost, upCost);
1871        else
1872            return CoinMax(downCost, upCost);
1873    }
1874}
1875
1876// Return "up" estimate
1877double
1878CbcSimpleIntegerPseudoCost::upEstimate() const
1879{
1880    OsiSolverInterface * solver = model_->solver();
1881    const double * solution = model_->testSolution();
1882    const double * lower = solver->getColLower();
1883    const double * upper = solver->getColUpper();
1884    double value = solution[columnNumber_];
1885    value = CoinMax(value, lower[columnNumber_]);
1886    value = CoinMin(value, upper[columnNumber_]);
1887    if (upper[columnNumber_] == lower[columnNumber_]) {
1888        // fixed
1889        return 0.0;
1890    }
1891    double integerTolerance =
1892        model_->getDblParam(CbcModel::CbcIntegerTolerance);
1893    double below = floor(value + integerTolerance);
1894    double above = below + 1.0;
1895    if (above > upper[columnNumber_]) {
1896        above = below;
1897        below = above - 1;
1898    }
1899    double upCost = CoinMax((above - value) * upPseudoCost_, 0.0);
1900    return upCost;
1901}
1902// Return "down" estimate
1903double
1904CbcSimpleIntegerPseudoCost::downEstimate() const
1905{
1906    OsiSolverInterface * solver = model_->solver();
1907    const double * solution = model_->testSolution();
1908    const double * lower = solver->getColLower();
1909    const double * upper = solver->getColUpper();
1910    double value = solution[columnNumber_];
1911    value = CoinMax(value, lower[columnNumber_]);
1912    value = CoinMin(value, upper[columnNumber_]);
1913    if (upper[columnNumber_] == lower[columnNumber_]) {
1914        // fixed
1915        return 0.0;
1916    }
1917    double integerTolerance =
1918        model_->getDblParam(CbcModel::CbcIntegerTolerance);
1919    double below = floor(value + integerTolerance);
1920    double above = below + 1.0;
1921    if (above > upper[columnNumber_]) {
1922        above = below;
1923        below = above - 1;
1924    }
1925    double downCost = CoinMax((value - below) * downPseudoCost_, 0.0);
1926    return downCost;
1927}
1928
1929//##############################################################################
1930
1931// Default Constructor
1932CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1933        : CbcIntegerBranchingObject()
1934{
1935    changeInGuessed_ = 1.0e-5;
1936}
1937
1938// Useful constructor
1939CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model,
1940        int variable, int way , double value)
1941        : CbcIntegerBranchingObject(model, variable, way, value)
1942{
1943}
1944// Useful constructor for fixing
1945CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model,
1946        int variable, int way,
1947        double lowerValue,
1948        double /*upperValue*/)
1949        : CbcIntegerBranchingObject(model, variable, way, lowerValue)
1950{
1951    changeInGuessed_ = 1.0e100;
1952}
1953
1954
1955// Copy constructor
1956CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (
1957    const CbcIntegerPseudoCostBranchingObject & rhs)
1958        : CbcIntegerBranchingObject(rhs)
1959{
1960    changeInGuessed_ = rhs.changeInGuessed_;
1961}
1962
1963// Assignment operator
1964CbcIntegerPseudoCostBranchingObject &
1965CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject & rhs)
1966{
1967    if (this != &rhs) {
1968        CbcIntegerBranchingObject::operator=(rhs);
1969        changeInGuessed_ = rhs.changeInGuessed_;
1970    }
1971    return *this;
1972}
1973CbcBranchingObject *
1974CbcIntegerPseudoCostBranchingObject::clone() const
1975{
1976    return (new CbcIntegerPseudoCostBranchingObject(*this));
1977}
1978
1979
1980// Destructor
1981CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
1982{
1983}
1984
1985/*
1986  Perform a branch by adjusting the bounds of the specified variable. Note
1987  that each arm of the branch advances the object to the next arm by
1988  advancing the value of way_.
1989
1990  Providing new values for the variable's lower and upper bounds for each
1991  branching direction gives a little bit of additional flexibility and will
1992  be easily extensible to multi-way branching.
1993  Returns change in guessed objective on next branch
1994*/
1995double
1996CbcIntegerPseudoCostBranchingObject::branch()
1997{
1998    CbcIntegerBranchingObject::branch();
1999    return changeInGuessed_;
2000}
2001
2002/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2003    same type and must have the same original object, but they may have
2004    different feasible regions.
2005    Return the appropriate CbcRangeCompare value (first argument being the
2006    sub/superset if that's the case). In case of overlap (and if \c
2007    replaceIfOverlap is true) replace the current branching object with one
2008    whose feasible region is the overlap.
2009*/
2010CbcRangeCompare
2011CbcIntegerPseudoCostBranchingObject::compareBranchingObject
2012(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
2013{
2014    const CbcIntegerPseudoCostBranchingObject* br =
2015        dynamic_cast<const CbcIntegerPseudoCostBranchingObject*>(brObj);
2016    assert(br);
2017    double* thisBd = way_ < 0 ? down_ : up_;
2018    const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
2019    return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
2020}
2021
2022
2023//##############################################################################
2024
2025// Default Constructor
2026CbcCliqueBranchingObject::CbcCliqueBranchingObject()
2027        : CbcBranchingObject()
2028{
2029    clique_ = NULL;
2030    downMask_[0] = 0;
2031    downMask_[1] = 0;
2032    upMask_[0] = 0;
2033    upMask_[1] = 0;
2034}
2035
2036// Useful constructor
2037CbcCliqueBranchingObject::CbcCliqueBranchingObject (CbcModel * model,
2038        const CbcClique * clique,
2039        int way ,
2040        int numberOnDownSide, const int * down,
2041        int numberOnUpSide, const int * up)
2042        : CbcBranchingObject(model, clique->id(), way, 0.5)
2043{
2044    clique_ = clique;
2045    downMask_[0] = 0;
2046    downMask_[1] = 0;
2047    upMask_[0] = 0;
2048    upMask_[1] = 0;
2049    int i;
2050    for (i = 0; i < numberOnDownSide; i++) {
2051        int sequence = down[i];
2052        int iWord = sequence >> 5;
2053        int iBit = sequence - 32 * iWord;
2054        unsigned int k = 1 << iBit;
2055        downMask_[iWord] |= k;
2056    }
2057    for (i = 0; i < numberOnUpSide; i++) {
2058        int sequence = up[i];
2059        int iWord = sequence >> 5;
2060        int iBit = sequence - 32 * iWord;
2061        unsigned int k = 1 << iBit;
2062        upMask_[iWord] |= k;
2063    }
2064}
2065
2066// Copy constructor
2067CbcCliqueBranchingObject::CbcCliqueBranchingObject ( const CbcCliqueBranchingObject & rhs) : CbcBranchingObject(rhs)
2068{
2069    clique_ = rhs.clique_;
2070    downMask_[0] = rhs.downMask_[0];
2071    downMask_[1] = rhs.downMask_[1];
2072    upMask_[0] = rhs.upMask_[0];
2073    upMask_[1] = rhs.upMask_[1];
2074}
2075
2076// Assignment operator
2077CbcCliqueBranchingObject &
2078CbcCliqueBranchingObject::operator=( const CbcCliqueBranchingObject & rhs)
2079{
2080    if (this != &rhs) {
2081        CbcBranchingObject::operator=(rhs);
2082        clique_ = rhs.clique_;
2083        downMask_[0] = rhs.downMask_[0];
2084        downMask_[1] = rhs.downMask_[1];
2085        upMask_[0] = rhs.upMask_[0];
2086        upMask_[1] = rhs.upMask_[1];
2087    }
2088    return *this;
2089}
2090CbcBranchingObject *
2091CbcCliqueBranchingObject::clone() const
2092{
2093    return (new CbcCliqueBranchingObject(*this));
2094}
2095
2096
2097// Destructor
2098CbcCliqueBranchingObject::~CbcCliqueBranchingObject ()
2099{
2100}
2101double
2102CbcCliqueBranchingObject::branch()
2103{
2104    decrementNumberBranchesLeft();
2105    int iWord;
2106    int numberMembers = clique_->numberMembers();
2107    const int * which = clique_->members();
2108    const int * integerVariables = model_->integerVariable();
2109    int numberWords = (numberMembers + 31) >> 5;
2110    // *** for way - up means fix all those in down section
2111    if (way_ < 0) {
2112#ifdef FULL_PRINT
2113        printf("Down Fix ");
2114#endif
2115        for (iWord = 0; iWord < numberWords; iWord++) {
2116            int i;
2117            for (i = 0; i < 32; i++) {
2118                unsigned int k = 1 << i;
2119                if ((upMask_[iWord]&k) != 0) {
2120                    int iColumn = which[i+32*iWord];
2121#ifdef FULL_PRINT
2122                    printf("%d ", i + 32*iWord);
2123#endif
2124                    // fix weak way
2125                    if (clique_->type(i + 32*iWord))
2126                        model_->solver()->setColUpper(integerVariables[iColumn], 0.0);
2127                    else
2128                        model_->solver()->setColLower(integerVariables[iColumn], 1.0);
2129                }
2130            }
2131        }
2132        way_ = 1;         // Swap direction
2133    } else {
2134#ifdef FULL_PRINT
2135        printf("Up Fix ");
2136#endif
2137        for (iWord = 0; iWord < numberWords; iWord++) {
2138            int i;
2139            for (i = 0; i < 32; i++) {
2140                unsigned int k = 1 << i;
2141                if ((downMask_[iWord]&k) != 0) {
2142                    int iColumn = which[i+32*iWord];
2143#ifdef FULL_PRINT
2144                    printf("%d ", i + 32*iWord);
2145#endif
2146                    // fix weak way
2147                    if (clique_->type(i + 32*iWord))
2148                        model_->solver()->setColUpper(integerVariables[iColumn], 0.0);
2149                    else
2150                        model_->solver()->setColLower(integerVariables[iColumn], 1.0);
2151                }
2152            }
2153        }
2154        way_ = -1;        // Swap direction
2155    }
2156#ifdef FULL_PRINT
2157    printf("\n");
2158#endif
2159    return 0.0;
2160}
2161// Print what would happen
2162void
2163CbcCliqueBranchingObject::print()
2164{
2165    int iWord;
2166    int numberMembers = clique_->numberMembers();
2167    const int * which = clique_->members();
2168    const int * integerVariables = model_->integerVariable();
2169    int numberWords = (numberMembers + 31) >> 5;
2170    // *** for way - up means fix all those in down section
2171    if (way_ < 0) {
2172        printf("Clique - Down Fix ");
2173        for (iWord = 0; iWord < numberWords; iWord++) {
2174            int i;
2175            for (i = 0; i < 32; i++) {
2176                unsigned int k = 1 << i;
2177                if ((upMask_[iWord]&k) != 0) {
2178                    int iColumn = which[i+32*iWord];
2179                    printf("%d ", integerVariables[iColumn]);
2180                }
2181            }
2182        }
2183    } else {
2184        printf("Clique - Up Fix ");
2185        for (iWord = 0; iWord < numberWords; iWord++) {
2186            int i;
2187            for (i = 0; i < 32; i++) {
2188                unsigned int k = 1 << i;
2189                if ((downMask_[iWord]&k) != 0) {
2190                    int iColumn = which[i+32*iWord];
2191                    printf("%d ", integerVariables[iColumn]);
2192                }
2193            }
2194        }
2195    }
2196    printf("\n");
2197}
2198
2199static inline int
2200CbcCompareCliques(const CbcClique* cl0, const CbcClique* cl1)
2201{
2202    if (cl0->cliqueType() < cl1->cliqueType()) {
2203        return -1;
2204    }
2205    if (cl0->cliqueType() > cl1->cliqueType()) {
2206        return 1;
2207    }
2208    if (cl0->numberMembers() != cl1->numberMembers()) {
2209        return cl0->numberMembers() - cl1->numberMembers();
2210    }
2211    if (cl0->numberNonSOSMembers() != cl1->numberNonSOSMembers()) {
2212        return cl0->numberNonSOSMembers() - cl1->numberNonSOSMembers();
2213    }
2214    return memcmp(cl0->members(), cl1->members(),
2215                  cl0->numberMembers() * sizeof(int));
2216}
2217
2218/** Compare the original object of \c this with the original object of \c
2219    brObj. Assumes that there is an ordering of the original objects.
2220    This method should be invoked only if \c this and brObj are of the same
2221    type.
2222    Return negative/0/positive depending on whether \c this is
2223    smaller/same/larger than the argument.
2224*/
2225int
2226CbcCliqueBranchingObject::compareOriginalObject
2227(const CbcBranchingObject* brObj) const
2228{
2229    const CbcCliqueBranchingObject* br =
2230        dynamic_cast<const CbcCliqueBranchingObject*>(brObj);
2231    assert(br);
2232    return CbcCompareCliques(clique_, br->clique_);
2233}
2234
2235/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2236    same type and must have the same original object, but they may have
2237    different feasible regions.
2238    Return the appropriate CbcRangeCompare value (first argument being the
2239    sub/superset if that's the case). In case of overlap (and if \c
2240    replaceIfOverlap is true) replace the current branching object with one
2241    whose feasible region is the overlap.
2242*/
2243CbcRangeCompare
2244CbcCliqueBranchingObject::compareBranchingObject
2245(const CbcBranchingObject* brObj, const bool /*replaceIfOverlap*/)
2246{
2247    const CbcCliqueBranchingObject* br =
2248        dynamic_cast<const CbcCliqueBranchingObject*>(brObj);
2249    assert(br);
2250    unsigned int* thisMask = way_ < 0 ? upMask_ : downMask_;
2251    const unsigned int* otherMask = br->way_ < 0 ? br->upMask_ : br->downMask_;
2252    const CoinUInt64 cl0 =
2253        (static_cast<CoinUInt64>(thisMask[0]) << 32) | thisMask[1];
2254    const CoinUInt64 cl1 =
2255        (static_cast<CoinUInt64>(otherMask[0]) << 32) | otherMask[1];
2256    if (cl0 == cl1) {
2257        return CbcRangeSame;
2258    }
2259    const CoinUInt64 cl_intersection = (cl0 & cl1);
2260    if (cl_intersection == cl0) {
2261        return CbcRangeSuperset;
2262    }
2263    if (cl_intersection == cl1) {
2264        return CbcRangeSubset;
2265    }
2266    const CoinUInt64 cl_xor = (cl0 ^ cl1);
2267    if (cl_intersection == 0 && cl_xor == 0) {
2268        return CbcRangeDisjoint;
2269    }
2270    const CoinUInt64 cl_union = (cl0 | cl1);
2271    thisMask[0] = static_cast<unsigned int>(cl_union >> 32);
2272    thisMask[1] = static_cast<unsigned int>(cl_union & 0xffffffff);
2273    return CbcRangeOverlap;
2274}
2275
2276//##############################################################################
2277
2278// Default Constructor
2279CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject()
2280        : CbcBranchingObject()
2281{
2282    clique_ = NULL;
2283    downMask_ = NULL;
2284    upMask_ = NULL;
2285}
2286
2287// Useful constructor
2288CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject (CbcModel * model,
2289        const CbcClique * clique,
2290        int way ,
2291        int numberOnDownSide, const int * down,
2292        int numberOnUpSide, const int * up)
2293        : CbcBranchingObject(model, clique->id(), way, 0.5)
2294{
2295    clique_ = clique;
2296    int numberMembers = clique_->numberMembers();
2297    int numberWords = (numberMembers + 31) >> 5;
2298    downMask_ = new unsigned int [numberWords];
2299    upMask_ = new unsigned int [numberWords];
2300    memset(downMask_, 0, numberWords*sizeof(unsigned int));
2301    memset(upMask_, 0, numberWords*sizeof(unsigned int));
2302    int i;
2303    for (i = 0; i < numberOnDownSide; i++) {
2304        int sequence = down[i];
2305        int iWord = sequence >> 5;
2306        int iBit = sequence - 32 * iWord;
2307        unsigned int k = 1 << iBit;
2308        downMask_[iWord] |= k;
2309    }
2310    for (i = 0; i < numberOnUpSide; i++) {
2311        int sequence = up[i];
2312        int iWord = sequence >> 5;
2313        int iBit = sequence - 32 * iWord;
2314        unsigned int k = 1 << iBit;
2315        upMask_[iWord] |= k;
2316    }
2317}
2318
2319// Copy constructor
2320CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject ( const CbcLongCliqueBranchingObject & rhs) : CbcBranchingObject(rhs)
2321{
2322    clique_ = rhs.clique_;
2323    if (rhs.downMask_) {
2324        int numberMembers = clique_->numberMembers();
2325        int numberWords = (numberMembers + 31) >> 5;
2326        downMask_ = new unsigned int [numberWords];
2327        memcpy(downMask_, rhs.downMask_, numberWords*sizeof(unsigned int));
2328        upMask_ = new unsigned int [numberWords];
2329        memcpy(upMask_, rhs.upMask_, numberWords*sizeof(unsigned int));
2330    } else {
2331        downMask_ = NULL;
2332        upMask_ = NULL;
2333    }
2334}
2335
2336// Assignment operator
2337CbcLongCliqueBranchingObject &
2338CbcLongCliqueBranchingObject::operator=( const CbcLongCliqueBranchingObject & rhs)
2339{
2340    if (this != &rhs) {
2341        CbcBranchingObject::operator=(rhs);
2342        clique_ = rhs.clique_;
2343        delete [] downMask_;
2344        delete [] upMask_;
2345        if (rhs.downMask_) {
2346            int numberMembers = clique_->numberMembers();
2347            int numberWords = (numberMembers + 31) >> 5;
2348            downMask_ = new unsigned int [numberWords];
2349            memcpy(downMask_, rhs.downMask_, numberWords*sizeof(unsigned int));
2350            upMask_ = new unsigned int [numberWords];
2351            memcpy(upMask_, rhs.upMask_, numberWords*sizeof(unsigned int));
2352        } else {
2353            downMask_ = NULL;
2354            upMask_ = NULL;
2355        }
2356    }
2357    return *this;
2358}
2359CbcBranchingObject *
2360CbcLongCliqueBranchingObject::clone() const
2361{
2362    return (new CbcLongCliqueBranchingObject(*this));
2363}
2364
2365
2366// Destructor
2367CbcLongCliqueBranchingObject::~CbcLongCliqueBranchingObject ()
2368{
2369    delete [] downMask_;
2370    delete [] upMask_;
2371}
2372double
2373CbcLongCliqueBranchingObject::branch()
2374{
2375    decrementNumberBranchesLeft();
2376    int iWord;
2377    int numberMembers = clique_->numberMembers();
2378    const int * which = clique_->members();
2379    const int * integerVariables = model_->integerVariable();
2380    int numberWords = (numberMembers + 31) >> 5;
2381    // *** for way - up means fix all those in down section
2382    if (way_ < 0) {
2383#ifdef FULL_PRINT
2384        printf("Down Fix ");
2385#endif
2386        for (iWord = 0; iWord < numberWords; iWord++) {
2387            int i;
2388            for (i = 0; i < 32; i++) {
2389                unsigned int k = 1 << i;
2390                if ((upMask_[iWord]&k) != 0) {
2391                    int iColumn = which[i+32*iWord];
2392#ifdef FULL_PRINT
2393                    printf("%d ", i + 32*iWord);
2394#endif
2395                    // fix weak way
2396                    if (clique_->type(i + 32*iWord))
2397                        model_->solver()->setColUpper(integerVariables[iColumn], 0.0);
2398                    else
2399                        model_->solver()->setColLower(integerVariables[iColumn], 1.0);
2400                }
2401            }
2402        }
2403        way_ = 1;         // Swap direction
2404    } else {
2405#ifdef FULL_PRINT
2406        printf("Up Fix ");
2407#endif
2408        for (iWord = 0; iWord < numberWords; iWord++) {
2409            int i;
2410            for (i = 0; i < 32; i++) {
2411                unsigned int k = 1 << i;
2412                if ((downMask_[iWord]&k) != 0) {
2413                    int iColumn = which[i+32*iWord];
2414#ifdef FULL_PRINT
2415                    printf("%d ", i + 32*iWord);
2416#endif
2417                    // fix weak way
2418                    if (clique_->type(i + 32*iWord))
2419                        model_->solver()->setColUpper(integerVariables[iColumn], 0.0);
2420                    else
2421                        model_->solver()->setColLower(integerVariables[iColumn], 1.0);
2422                }
2423            }
2424        }
2425        way_ = -1;        // Swap direction
2426    }
2427#ifdef FULL_PRINT
2428    printf("\n");
2429#endif
2430    return 0.0;
2431}
2432void
2433CbcLongCliqueBranchingObject::print()
2434{
2435    int iWord;
2436    int numberMembers = clique_->numberMembers();
2437    const int * which = clique_->members();
2438    const int * integerVariables = model_->integerVariable();
2439    int numberWords = (numberMembers + 31) >> 5;
2440    // *** for way - up means fix all those in down section
2441    if (way_ < 0) {
2442        printf("Clique - Down Fix ");
2443        for (iWord = 0; iWord < numberWords; iWord++) {
2444            int i;
2445            for (i = 0; i < 32; i++) {
2446                unsigned int k = 1 << i;
2447                if ((upMask_[iWord]&k) != 0) {
2448                    int iColumn = which[i+32*iWord];
2449                    printf("%d ", integerVariables[iColumn]);
2450                }
2451            }
2452        }
2453    } else {
2454        printf("Clique - Up Fix ");
2455        for (iWord = 0; iWord < numberWords; iWord++) {
2456            int i;
2457            for (i = 0; i < 32; i++) {
2458                unsigned int k = 1 << i;
2459                if ((downMask_[iWord]&k) != 0) {
2460                    int iColumn = which[i+32*iWord];
2461                    printf("%d ", integerVariables[iColumn]);
2462                }
2463            }
2464        }
2465    }
2466    printf("\n");
2467}
2468
2469/** Compare the original object of \c this with the original object of \c
2470    brObj. Assumes that there is an ordering of the original objects.
2471    This method should be invoked only if \c this and brObj are of the same
2472    type.
2473    Return negative/0/positive depending on whether \c this is
2474    smaller/same/larger than the argument.
2475*/
2476int
2477CbcLongCliqueBranchingObject::compareOriginalObject
2478(const CbcBranchingObject* brObj) const
2479{
2480    const CbcLongCliqueBranchingObject* br =
2481        dynamic_cast<const CbcLongCliqueBranchingObject*>(brObj);
2482    assert(br);
2483    return CbcCompareCliques(clique_, br->clique_);
2484}
2485
2486/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2487    same type and must have the same original object, but they may have
2488    different feasible regions.
2489    Return the appropriate CbcRangeCompare value (first argument being the
2490    sub/superset if that's the case). In case of overlap (and if \c
2491    replaceIfOverlap is true) replace the current branching object with one
2492    whose feasible region is the overlap.
2493*/
2494CbcRangeCompare
2495CbcLongCliqueBranchingObject::compareBranchingObject
2496(const CbcBranchingObject* brObj, const bool /*replaceIfOverlap*/)
2497{
2498    const CbcLongCliqueBranchingObject* br =
2499        dynamic_cast<const CbcLongCliqueBranchingObject*>(brObj);
2500    assert(br);
2501    const int numberMembers = clique_->numberMembers();
2502    const int numberWords = (numberMembers + 31) >> 5;
2503    unsigned int* thisMask = way_ < 0 ? upMask_ : downMask_;
2504    const unsigned int* otherMask = br->way_ < 0 ? br->upMask_ : br->downMask_;
2505
2506    if (memcmp(thisMask, otherMask, numberWords * sizeof(unsigned int)) == 0) {
2507        return CbcRangeSame;
2508    }
2509    bool canBeSuperset = true;
2510    bool canBeSubset = true;
2511    int i;
2512    for (i = numberWords - 1; i >= 0 && (canBeSuperset || canBeSubset); --i) {
2513        const unsigned int both = (thisMask[i] & otherMask[i]);
2514        canBeSuperset &= (both == thisMask[i]);
2515        canBeSubset &= (both == otherMask[i]);
2516    }
2517    if (canBeSuperset) {
2518        return CbcRangeSuperset;
2519    }
2520    if (canBeSubset) {
2521        return CbcRangeSubset;
2522    }
2523
2524    for (i = numberWords - 1; i >= 0; --i) {
2525        if ((thisMask[i] ^ otherMask[i]) != 0) {
2526            break;
2527        }
2528    }
2529    if (i == -1) { // complement
2530        return CbcRangeDisjoint;
2531    }
2532    // must be overlap
2533    for (i = numberWords - 1; i >= 0; --i) {
2534        thisMask[i] |= otherMask[i];
2535    }
2536    return CbcRangeOverlap;
2537}
2538
2539//##############################################################################
2540
2541// Default Constructor
2542CbcSOSBranchingObject::CbcSOSBranchingObject()
2543        : CbcBranchingObject(),
2544        firstNonzero_(-1),
2545        lastNonzero_(-1)
2546{
2547    set_ = NULL;
2548    separator_ = 0.0;
2549}
2550
2551// Useful constructor
2552CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
2553        const CbcSOS * set,
2554        int way ,
2555        double separator)
2556        : CbcBranchingObject(model, set->id(), way, 0.5)
2557{
2558    set_ = set;
2559    separator_ = separator;
2560    computeNonzeroRange();
2561}
2562
2563// Copy constructor
2564CbcSOSBranchingObject::CbcSOSBranchingObject (const CbcSOSBranchingObject & rhs)
2565        : CbcBranchingObject(rhs),
2566        firstNonzero_(rhs.firstNonzero_),
2567        lastNonzero_(rhs.lastNonzero_)
2568{
2569    set_ = rhs.set_;
2570    separator_ = rhs.separator_;
2571}
2572
2573// Assignment operator
2574CbcSOSBranchingObject &
2575CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject & rhs)
2576{
2577    if (this != &rhs) {
2578        CbcBranchingObject::operator=(rhs);
2579        set_ = rhs.set_;
2580        separator_ = rhs.separator_;
2581        firstNonzero_ = rhs.firstNonzero_;
2582        lastNonzero_ = rhs.lastNonzero_;
2583    }
2584    return *this;
2585}
2586CbcBranchingObject *
2587CbcSOSBranchingObject::clone() const
2588{
2589    return (new CbcSOSBranchingObject(*this));
2590}
2591
2592
2593// Destructor
2594CbcSOSBranchingObject::~CbcSOSBranchingObject ()
2595{
2596}
2597
2598void
2599CbcSOSBranchingObject::computeNonzeroRange()
2600{
2601    const int numberMembers = set_->numberMembers();
2602    const double * weights = set_->weights();
2603    int i = 0;
2604    if (way_ < 0) {
2605        for ( i = 0; i < numberMembers; i++) {
2606            if (weights[i] > separator_)
2607                break;
2608        }
2609        assert (i < numberMembers);
2610        firstNonzero_ = 0;
2611        lastNonzero_ = i;
2612    } else {
2613        for ( i = 0; i < numberMembers; i++) {
2614            if (weights[i] >= separator_)
2615                break;
2616        }
2617        assert (i < numberMembers);
2618        firstNonzero_ = i;
2619        lastNonzero_ = numberMembers;
2620    }
2621}
2622
2623double
2624CbcSOSBranchingObject::branch()
2625{
2626    decrementNumberBranchesLeft();
2627    int numberMembers = set_->numberMembers();
2628    const int * which = set_->members();
2629    const double * weights = set_->weights();
2630    OsiSolverInterface * solver = model_->solver();
2631    //const double * lower = solver->getColLower();
2632    //const double * upper = solver->getColUpper();
2633    // *** for way - up means fix all those in down section
2634    if (way_ < 0) {
2635        int i;
2636        for ( i = 0; i < numberMembers; i++) {
2637            if (weights[i] > separator_)
2638                break;
2639        }
2640        assert (i < numberMembers);
2641        for (; i < numberMembers; i++)
2642            solver->setColUpper(which[i], 0.0);
2643        way_ = 1;         // Swap direction
2644    } else {
2645        int i;
2646        for ( i = 0; i < numberMembers; i++) {
2647            if (weights[i] >= separator_)
2648                break;
2649            else
2650                solver->setColUpper(which[i], 0.0);
2651        }
2652        assert (i < numberMembers);
2653        way_ = -1;        // Swap direction
2654    }
2655    computeNonzeroRange();
2656    return 0.0;
2657}
2658/* Update bounds in solver as in 'branch' and update given bounds.
2659   branchState is -1 for 'down' +1 for 'up' */
2660void
2661CbcSOSBranchingObject::fix(OsiSolverInterface * solver,
2662                           double * /*lower*/, double * upper,
2663                           int branchState) const
2664{
2665    int numberMembers = set_->numberMembers();
2666    const int * which = set_->members();
2667    const double * weights = set_->weights();
2668    //const double * lower = solver->getColLower();
2669    //const double * upper = solver->getColUpper();
2670    // *** for way - up means fix all those in down section
2671    if (branchState < 0) {
2672        int i;
2673        for ( i = 0; i < numberMembers; i++) {
2674            if (weights[i] > separator_)
2675                break;
2676        }
2677        assert (i < numberMembers);
2678        for (; i < numberMembers; i++) {
2679            solver->setColUpper(which[i], 0.0);
2680            upper[which[i]] = 0.0;
2681        }
2682    } else {
2683        int i;
2684        for ( i = 0; i < numberMembers; i++) {
2685            if (weights[i] >= separator_) {
2686                break;
2687            } else {
2688                solver->setColUpper(which[i], 0.0);
2689                upper[which[i]] = 0.0;
2690            }
2691        }
2692        assert (i < numberMembers);
2693    }
2694}
2695// Print what would happen
2696void
2697CbcSOSBranchingObject::print()
2698{
2699    int numberMembers = set_->numberMembers();
2700    const int * which = set_->members();
2701    const double * weights = set_->weights();
2702    OsiSolverInterface * solver = model_->solver();
2703    //const double * lower = solver->getColLower();
2704    const double * upper = solver->getColUpper();
2705    int first = numberMembers;
2706    int last = -1;
2707    int numberFixed = 0;
2708    int numberOther = 0;
2709    int i;
2710    for ( i = 0; i < numberMembers; i++) {
2711        double bound = upper[which[i]];
2712        if (bound) {
2713            first = CoinMin(first, i);
2714            last = CoinMax(last, i);
2715        }
2716    }
2717    // *** for way - up means fix all those in down section
2718    if (way_ < 0) {
2719        printf("SOS Down");
2720        for ( i = 0; i < numberMembers; i++) {
2721            double bound = upper[which[i]];
2722            if (weights[i] > separator_)
2723                break;
2724            else if (bound)
2725                numberOther++;
2726        }
2727        assert (i < numberMembers);
2728        for (; i < numberMembers; i++) {
2729            double bound = upper[which[i]];
2730            if (bound)
2731                numberFixed++;
2732        }
2733    } else {
2734        printf("SOS Up");
2735        for ( i = 0; i < numberMembers; i++) {
2736            double bound = upper[which[i]];
2737            if (weights[i] >= separator_)
2738                break;
2739            else if (bound)
2740                numberFixed++;
2741        }
2742        assert (i < numberMembers);
2743        for (; i < numberMembers; i++) {
2744            double bound = upper[which[i]];
2745            if (bound)
2746                numberOther++;
2747        }
2748    }
2749    printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
2750           separator_, which[first], weights[first], which[last], weights[last], numberFixed, numberOther);
2751}
2752
2753/** Compare the original object of \c this with the original object of \c
2754    brObj. Assumes that there is an ordering of the original objects.
2755    This method should be invoked only if \c this and brObj are of the same
2756    type.
2757    Return negative/0/positive depending on whether \c this is
2758    smaller/same/larger than the argument.
2759*/
2760int
2761CbcSOSBranchingObject::compareOriginalObject
2762(const CbcBranchingObject* brObj) const
2763{
2764    const CbcSOSBranchingObject* br =
2765        dynamic_cast<const CbcSOSBranchingObject*>(brObj);
2766    assert(br);
2767    const CbcSOS* s0 = set_;
2768    const CbcSOS* s1 = br->set_;
2769    if (s0->sosType() != s1->sosType()) {
2770        return s0->sosType() - s1->sosType();
2771    }
2772    if (s0->numberMembers() != s1->numberMembers()) {
2773        return s0->numberMembers() - s1->numberMembers();
2774    }
2775    const int memberCmp = memcmp(s0->members(), s1->members(),
2776                                 s0->numberMembers() * sizeof(int));
2777    if (memberCmp != 0) {
2778        return memberCmp;
2779    }
2780    return memcmp(s0->weights(), s1->weights(),
2781                  s0->numberMembers() * sizeof(double));
2782}
2783
2784/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2785    same type and must have the same original object, but they may have
2786    different feasible regions.
2787    Return the appropriate CbcRangeCompare value (first argument being the
2788    sub/superset if that's the case). In case of overlap (and if \c
2789    replaceIfOverlap is true) replace the current branching object with one
2790    whose feasible region is the overlap.
2791*/
2792CbcRangeCompare
2793CbcSOSBranchingObject::compareBranchingObject
2794(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
2795{
2796    const CbcSOSBranchingObject* br =
2797        dynamic_cast<const CbcSOSBranchingObject*>(brObj);
2798    assert(br);
2799    if (firstNonzero_ < br->firstNonzero_) {
2800        if (lastNonzero_ >= br->lastNonzero_) {
2801            return CbcRangeSuperset;
2802        } else if (lastNonzero_ <= br->firstNonzero_) {
2803            return CbcRangeDisjoint;
2804        } else {
2805            // overlap
2806            if (replaceIfOverlap) {
2807                firstNonzero_ = br->firstNonzero_;
2808            }
2809            return CbcRangeOverlap;
2810        }
2811    } else if (firstNonzero_ > br->firstNonzero_) {
2812        if (lastNonzero_ <= br->lastNonzero_) {
2813            return CbcRangeSubset;
2814        } else if (firstNonzero_ >= br->lastNonzero_) {
2815            return CbcRangeDisjoint;
2816        } else {
2817            // overlap
2818            if (replaceIfOverlap) {
2819                lastNonzero_ = br->lastNonzero_;
2820            }
2821            return CbcRangeOverlap;
2822        }
2823    } else {
2824        if (lastNonzero_ == br->lastNonzero_) {
2825            return CbcRangeSame;
2826        }
2827        return lastNonzero_ < br->lastNonzero_ ? CbcRangeSubset : CbcRangeSuperset;
2828    }
2829    return CbcRangeSame; // fake return
2830}
2831
2832//##############################################################################
2833
2834// Default Constructor
2835CbcBranchDefaultDecision::CbcBranchDefaultDecision()
2836        : CbcBranchDecision()
2837{
2838    bestCriterion_ = 0.0;
2839    bestChangeUp_ = 0.0;
2840    bestNumberUp_ = 0;
2841    bestChangeDown_ = 0.0;
2842    bestObject_ = NULL;
2843    bestNumberDown_ = 0;
2844}
2845
2846// Copy constructor
2847CbcBranchDefaultDecision::CbcBranchDefaultDecision (
2848    const CbcBranchDefaultDecision & rhs)
2849        : CbcBranchDecision(rhs)
2850{
2851    bestCriterion_ = rhs.bestCriterion_;
2852    bestChangeUp_ = rhs.bestChangeUp_;
2853    bestNumberUp_ = rhs.bestNumberUp_;
2854    bestChangeDown_ = rhs.bestChangeDown_;
2855    bestNumberDown_ = rhs.bestNumberDown_;
2856    bestObject_ = rhs.bestObject_;
2857    model_ = rhs.model_;
2858}
2859
2860CbcBranchDefaultDecision::~CbcBranchDefaultDecision()
2861{
2862}
2863
2864// Clone
2865CbcBranchDecision *
2866CbcBranchDefaultDecision::clone() const
2867{
2868    return new CbcBranchDefaultDecision(*this);
2869}
2870
2871// Initialize i.e. before start of choosing at a node
2872void
2873CbcBranchDefaultDecision::initialize(CbcModel * model)
2874{
2875    bestCriterion_ = 0.0;
2876    bestChangeUp_ = 0.0;
2877    bestNumberUp_ = 0;
2878    bestChangeDown_ = 0.0;
2879    bestNumberDown_ = 0;
2880    bestObject_ = NULL;
2881    model_ = model;
2882}
2883
2884
2885/*
2886  Simple default decision algorithm. Compare based on infeasibility (numInfUp,
2887  numInfDn) until a solution is found by search, then switch to change in
2888  objective (changeUp, changeDn). Note that bestSoFar is remembered in
2889  bestObject_, so the parameter bestSoFar is unused.
2890*/
2891
2892int
2893CbcBranchDefaultDecision::betterBranch(CbcBranchingObject * thisOne,
2894                                       CbcBranchingObject * /*bestSoFar*/,
2895                                       double changeUp, int numInfUp,
2896                                       double changeDn, int numInfDn)
2897{
2898    bool beforeSolution = cbcModel()->getSolutionCount() ==
2899                          cbcModel()->getNumberHeuristicSolutions();;
2900    int betterWay = 0;
2901    if (beforeSolution) {
2902        if (!bestObject_) {
2903            bestNumberUp_ = COIN_INT_MAX;
2904            bestNumberDown_ = COIN_INT_MAX;
2905        }
2906        // before solution - choose smallest number
2907        // could add in depth as well
2908        int bestNumber = CoinMin(bestNumberUp_, bestNumberDown_);
2909        if (numInfUp < numInfDn) {
2910            if (numInfUp < bestNumber) {
2911                betterWay = 1;
2912            } else if (numInfUp == bestNumber) {
2913                if (changeUp < bestCriterion_)
2914                    betterWay = 1;
2915            }
2916        } else if (numInfUp > numInfDn) {
2917            if (numInfDn < bestNumber) {
2918                betterWay = -1;
2919            } else if (numInfDn == bestNumber) {
2920                if (changeDn < bestCriterion_)
2921                    betterWay = -1;
2922            }
2923        } else {
2924            // up and down have same number
2925            bool better = false;
2926            if (numInfUp < bestNumber) {
2927                better = true;
2928            } else if (numInfUp == bestNumber) {
2929                if (CoinMin(changeUp, changeDn) < bestCriterion_)
2930                    better = true;;
2931            }
2932            if (better) {
2933                // see which way
2934                if (changeUp <= changeDn)
2935                    betterWay = 1;
2936                else
2937                    betterWay = -1;
2938            }
2939        }
2940    } else {
2941        if (!bestObject_) {
2942            bestCriterion_ = -1.0;
2943        }
2944        // got a solution
2945        if (changeUp <= changeDn) {
2946            if (changeUp > bestCriterion_)
2947                betterWay = 1;
2948        } else {
2949            if (changeDn > bestCriterion_)
2950                betterWay = -1;
2951        }
2952    }
2953    if (betterWay) {
2954        bestCriterion_ = CoinMin(changeUp, changeDn);
2955        bestChangeUp_ = changeUp;
2956        bestNumberUp_ = numInfUp;
2957        bestChangeDown_ = changeDn;
2958        bestNumberDown_ = numInfDn;
2959        bestObject_ = thisOne;
2960        // See if user is overriding way
2961        if (thisOne->object() && thisOne->object()->preferredWay())
2962            betterWay = thisOne->object()->preferredWay();
2963    }
2964    return betterWay;
2965}
2966/* Sets or gets best criterion so far */
2967void
2968CbcBranchDefaultDecision::setBestCriterion(double value)
2969{
2970    bestCriterion_ = value;
2971}
2972double
2973CbcBranchDefaultDecision::getBestCriterion() const
2974{
2975    return bestCriterion_;
2976}
2977
2978/* Compare N branching objects. Return index of best
2979   and sets way of branching in chosen object.
2980
2981   This routine is used only after strong branching.
2982*/
2983
2984int
2985CbcBranchDefaultDecision::bestBranch (CbcBranchingObject ** objects, int numberObjects,
2986                                      int numberUnsatisfied,
2987                                      double * changeUp, int * numberInfeasibilitiesUp,
2988                                      double * changeDown, int * numberInfeasibilitiesDown,
2989                                      double objectiveValue)
2990{
2991
2992    int bestWay = 0;
2993    int whichObject = -1;
2994    if (numberObjects) {
2995        CbcModel * model = cbcModel();
2996        // at continuous
2997        //double continuousObjective = model->getContinuousObjective();
2998        //int continuousInfeasibilities = model->getContinuousInfeasibilities();
2999
3000        // average cost to get rid of infeasibility
3001        //double averageCostPerInfeasibility =
3002        //(objectiveValue-continuousObjective)/
3003        //(double) (abs(continuousInfeasibilities-numberUnsatisfied)+1);
3004        /* beforeSolution is :
3005           0 - before any solution
3006           n - n heuristic solutions but no branched one
3007           -1 - branched solution found
3008        */
3009        int numberSolutions = model->getSolutionCount();
3010        double cutoff = model->getCutoff();
3011        int method = 0;
3012        int i;
3013        if (numberSolutions) {
3014            int numberHeuristic = model->getNumberHeuristicSolutions();
3015            if (numberHeuristic < numberSolutions) {
3016                method = 1;
3017            } else {
3018                method = 2;
3019                // look further
3020                for ( i = 0 ; i < numberObjects ; i++) {
3021                    int numberNext = numberInfeasibilitiesUp[i];
3022
3023                    if (numberNext < numberUnsatisfied) {
3024                        int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
3025                        double perUnsatisfied = changeUp[i] / static_cast<double> (numberUp);
3026                        double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
3027                        if (estimatedObjective < cutoff)
3028                            method = 3;
3029                    }
3030                    numberNext = numberInfeasibilitiesDown[i];
3031                    if (numberNext < numberUnsatisfied) {
3032                        int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
3033                        double perUnsatisfied = changeDown[i] / static_cast<double> (numberDown);
3034                        double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
3035                        if (estimatedObjective < cutoff)
3036                            method = 3;
3037                    }
3038                }
3039            }
3040            method = 2;
3041        } else {
3042            method = 0;
3043        }
3044        // Uncomment next to force method 4
3045        //method=4;
3046        /* Methods :
3047           0 - fewest infeasibilities
3048           1 - largest min change in objective
3049           2 - as 1 but use sum of changes if min close
3050           3 - predicted best solution
3051           4 - take cheapest up branch if infeasibilities same
3052        */
3053        int bestNumber = COIN_INT_MAX;
3054        double bestCriterion = -1.0e50;
3055        double alternativeCriterion = -1.0;
3056        double bestEstimate = 1.0e100;
3057        switch (method) {
3058        case 0:
3059            // could add in depth as well
3060            for ( i = 0 ; i < numberObjects ; i++) {
3061                int thisNumber = CoinMin(numberInfeasibilitiesUp[i], numberInfeasibilitiesDown[i]);
3062                if (thisNumber <= bestNumber) {
3063                    int betterWay = 0;
3064                    if (numberInfeasibilitiesUp[i] < numberInfeasibilitiesDown[i]) {
3065                        if (numberInfeasibilitiesUp[i] < bestNumber) {
3066                            betterWay = 1;
3067                        } else {
3068                            if (changeUp[i] < bestCriterion)
3069                                betterWay = 1;
3070                        }
3071                    } else if (numberInfeasibilitiesUp[i] > numberInfeasibilitiesDown[i]) {
3072                        if (numberInfeasibilitiesDown[i] < bestNumber) {
3073                            betterWay = -1;
3074                        } else {
3075                            if (changeDown[i] < bestCriterion)
3076                                betterWay = -1;
3077                        }
3078                    } else {
3079                        // up and down have same number
3080                        bool better = false;
3081                        if (numberInfeasibilitiesUp[i] < bestNumber) {
3082                            better = true;
3083                        } else if (numberInfeasibilitiesUp[i] == bestNumber) {
3084                            if (CoinMin(changeUp[i], changeDown[i]) < bestCriterion)
3085                                better = true;;
3086                        }
3087                        if (better) {
3088                            // see which way
3089                            if (changeUp[i] <= changeDown[i])
3090                                betterWay = 1;
3091                            else
3092                                betterWay = -1;
3093                        }
3094                    }
3095                    if (betterWay) {
3096                        bestCriterion = CoinMin(changeUp[i], changeDown[i]);
3097                        bestNumber = thisNumber;
3098                        whichObject = i;
3099                        bestWay = betterWay;
3100                    }
3101                }
3102            }
3103            break;
3104        case 1:
3105            for ( i = 0 ; i < numberObjects ; i++) {
3106                int betterWay = 0;
3107                if (changeUp[i] <= changeDown[i]) {
3108                    if (changeUp[i] > bestCriterion)
3109                        betterWay = 1;
3110                } else {
3111                    if (changeDown[i] > bestCriterion)
3112                        betterWay = -1;
3113                }
3114                if (betterWay) {
3115                    bestCriterion = CoinMin(changeUp[i], changeDown[i]);
3116                    whichObject = i;
3117                    bestWay = betterWay;
3118                }
3119            }
3120            break;
3121        case 2:
3122            for ( i = 0 ; i < numberObjects ; i++) {
3123                double change = CoinMin(changeUp[i], changeDown[i]);
3124                double sum = changeUp[i] + changeDown[i];
3125                bool take = false;
3126                if (change > 1.1*bestCriterion)
3127                    take = true;
3128                else if (change > 0.9*bestCriterion && sum + change > bestCriterion + alternativeCriterion)
3129                    take = true;
3130                if (take) {
3131                    if (changeUp[i] <= changeDown[i]) {
3132                        if (changeUp[i] > bestCriterion)
3133                            bestWay = 1;
3134                    } else {
3135                        if (changeDown[i] > bestCriterion)
3136                            bestWay = -1;
3137                    }
3138                    bestCriterion = change;
3139                    alternativeCriterion = sum;
3140                    whichObject = i;
3141                }
3142            }
3143            break;
3144        case 3:
3145            for ( i = 0 ; i < numberObjects ; i++) {
3146                int numberNext = numberInfeasibilitiesUp[i];
3147
3148                if (numberNext < numberUnsatisfied) {
3149                    int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
3150                    double perUnsatisfied = changeUp[i] / static_cast<double> (numberUp);
3151                    double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
3152                    if (estimatedObjective < bestEstimate) {
3153                        bestEstimate = estimatedObjective;
3154                        bestWay = 1;
3155                        whichObject = i;
3156                    }
3157                }
3158                numberNext = numberInfeasibilitiesDown[i];
3159                if (numberNext < numberUnsatisfied) {
3160                    int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
3161                    double perUnsatisfied = changeDown[i] / static_cast<double> (numberDown);
3162                    double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
3163                    if (estimatedObjective < bestEstimate) {
3164                        bestEstimate = estimatedObjective;
3165                        bestWay = -1;
3166                        whichObject = i;
3167                    }
3168                }
3169            }
3170            break;
3171        case 4:
3172            // if number infeas same then cheapest up
3173            // first get best number or when going down
3174            // now choose smallest change up amongst equal number infeas
3175            for ( i = 0 ; i < numberObjects ; i++) {
3176                int thisNumber = CoinMin(numberInfeasibilitiesUp[i], numberInfeasibilitiesDown[i]);
3177                if (thisNumber <= bestNumber) {
3178                    int betterWay = 0;
3179                    if (numberInfeasibilitiesUp[i] < numberInfeasibilitiesDown[i]) {
3180                        if (numberInfeasibilitiesUp[i] < bestNumber) {
3181                            betterWay = 1;
3182                        } else {
3183                            if (changeUp[i] < bestCriterion)
3184                                betterWay = 1;
3185                        }
3186                    } else if (numberInfeasibilitiesUp[i] > numberInfeasibilitiesDown[i]) {
3187                        if (numberInfeasibilitiesDown[i] < bestNumber) {
3188                            betterWay = -1;
3189                        } else {
3190                            if (changeDown[i] < bestCriterion)
3191                                betterWay = -1;
3192                        }
3193                    } else {
3194                        // up and down have same number
3195                        bool better = false;
3196                        if (numberInfeasibilitiesUp[i] < bestNumber) {
3197                            better = true;
3198                        } else if (numberInfeasibilitiesUp[i] == bestNumber) {
3199                            if (CoinMin(changeUp[i], changeDown[i]) < bestCriterion)
3200                                better = true;;
3201                        }
3202                        if (better) {
3203                            // see which way
3204                            if (changeUp[i] <= changeDown[i])
3205                                betterWay = 1;
3206                            else
3207                                betterWay = -1;
3208                        }
3209                    }
3210                    if (betterWay) {
3211                        bestCriterion = CoinMin(changeUp[i], changeDown[i]);
3212                        bestNumber = thisNumber;
3213                        whichObject = i;
3214                        bestWay = betterWay;
3215                    }
3216                }
3217            }
3218            bestCriterion = 1.0e50;
3219            for ( i = 0 ; i < numberObjects ; i++) {
3220                int thisNumber = numberInfeasibilitiesUp[i];
3221                if (thisNumber == bestNumber && changeUp) {
3222                    if (changeUp[i] < bestCriterion) {
3223                        bestCriterion = changeUp[i];
3224                        whichObject = i;
3225                        bestWay = 1;
3226                    }
3227                }
3228            }
3229            break;
3230        }
3231        // set way in best
3232        if (whichObject >= 0) {
3233            CbcBranchingObject * bestObject = objects[whichObject];
3234            if (bestObject->object() && bestObject->object()->preferredWay())
3235                bestWay = bestObject->object()->preferredWay();
3236            bestObject->way(bestWay);
3237        } else {
3238            printf("debug\n");
3239        }
3240    }
3241    return whichObject;
3242}
3243
3244//##############################################################################
3245
3246// Default Constructor
3247CbcFollowOn::CbcFollowOn ()
3248        : CbcObject(),
3249        rhs_(NULL)
3250{
3251}
3252
3253// Useful constructor
3254CbcFollowOn::CbcFollowOn (CbcModel * model)
3255        : CbcObject(model)
3256{
3257    assert (model);
3258    OsiSolverInterface * solver = model_->solver();
3259    matrix_ = *solver->getMatrixByCol();
3260    matrix_.removeGaps();
3261    matrix_.setExtraGap(0.0);
3262    matrixByRow_ = *solver->getMatrixByRow();
3263    int numberRows = matrix_.getNumRows();
3264
3265    rhs_ = new int[numberRows];
3266    int i;
3267    const double * rowLower = solver->getRowLower();
3268    const double * rowUpper = solver->getRowUpper();
3269    // Row copy
3270    const double * elementByRow = matrixByRow_.getElements();
3271    const int * column = matrixByRow_.getIndices();
3272    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
3273    const int * rowLength = matrixByRow_.getVectorLengths();
3274    for (i = 0; i < numberRows; i++) {
3275        rhs_[i] = 0;
3276        double value = rowLower[i];
3277        if (value == rowUpper[i]) {
3278            if (floor(value) == value && value >= 1.0 && value < 10.0) {
3279                // check elements
3280                bool good = true;
3281                for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
3282                    int iColumn = column[j];
3283                    if (!solver->isBinary(iColumn))
3284                        good = false;
3285                    double elValue = elementByRow[j];
3286                    if (floor(elValue) != elValue || value < 1.0)
3287                        good = false;
3288                }
3289                if (good)
3290                    rhs_[i] = static_cast<int> (value);
3291            }
3292        }
3293    }
3294}
3295
3296// Copy constructor
3297CbcFollowOn::CbcFollowOn ( const CbcFollowOn & rhs)
3298        : CbcObject(rhs),
3299        matrix_(rhs.matrix_),
3300        matrixByRow_(rhs.matrixByRow_)
3301{
3302    int numberRows = matrix_.getNumRows();
3303    rhs_ = CoinCopyOfArray(rhs.rhs_, numberRows);
3304}
3305
3306// Clone
3307CbcObject *
3308CbcFollowOn::clone() const
3309{
3310    return new CbcFollowOn(*this);
3311}
3312
3313// Assignment operator
3314CbcFollowOn &
3315CbcFollowOn::operator=( const CbcFollowOn & rhs)
3316{
3317    if (this != &rhs) {
3318        CbcObject::operator=(rhs);
3319        delete [] rhs_;
3320        matrix_ = rhs.matrix_;
3321        matrixByRow_ = rhs.matrixByRow_;
3322        int numberRows = matrix_.getNumRows();
3323        rhs_ = CoinCopyOfArray(rhs.rhs_, numberRows);
3324    }
3325    return *this;
3326}
3327
3328// Destructor
3329CbcFollowOn::~CbcFollowOn ()
3330{
3331    delete [] rhs_;
3332}
3333// As some computation is needed in more than one place - returns row
3334int
3335CbcFollowOn::gutsOfFollowOn(int & otherRow, int & preferredWay) const
3336{
3337    int whichRow = -1;
3338    otherRow = -1;
3339    int numberRows = matrix_.getNumRows();
3340
3341    int i;
3342    // For sorting
3343    int * sort = new int [numberRows];
3344    int * isort = new int [numberRows];
3345    // Column copy
3346    //const double * element = matrix_.getElements();
3347    const int * row = matrix_.getIndices();
3348    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
3349    const int * columnLength = matrix_.getVectorLengths();
3350    // Row copy
3351    const double * elementByRow = matrixByRow_.getElements();
3352    const int * column = matrixByRow_.getIndices();
3353    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
3354    const int * rowLength = matrixByRow_.getVectorLengths();
3355    OsiSolverInterface * solver = model_->solver();
3356    const double * columnLower = solver->getColLower();
3357    const double * columnUpper = solver->getColUpper();
3358    const double * solution = solver->getColSolution();
3359    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
3360    int nSort = 0;
3361    for (i = 0; i < numberRows; i++) {
3362        if (rhs_[i]) {
3363            // check elements
3364            double smallest = 1.0e10;
3365            double largest = 0.0;
3366            int rhsValue = rhs_[i];
3367            int number1 = 0;
3368            int numberUnsatisfied = 0;
3369            for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
3370                int iColumn = column[j];
3371                double value = elementByRow[j];
3372                double solValue = solution[iColumn];
3373                if (columnLower[iColumn] != columnUpper[iColumn]) {
3374                    smallest = CoinMin(smallest, value);
3375                    largest = CoinMax(largest, value);
3376                    if (value == 1.0)
3377                        number1++;
3378                    if (solValue < 1.0 - integerTolerance && solValue > integerTolerance)
3379                        numberUnsatisfied++;
3380                } else {
3381                    rhsValue -= static_cast<int>(value * floor(solValue + 0.5));
3382                }
3383            }
3384            if (numberUnsatisfied > 1) {
3385                if (smallest < largest) {
3386                    // probably no good but check a few things
3387                    assert (largest <= rhsValue);
3388                    if (number1 == 1 && largest == rhsValue)
3389                        printf("could fix\n");
3390                } else if (largest == rhsValue) {
3391                    sort[nSort] = i;
3392                    isort[nSort++] = -numberUnsatisfied;
3393                }
3394            }
3395        }
3396    }
3397    if (nSort > 1) {
3398        CoinSort_2(isort, isort + nSort, sort);
3399        CoinZeroN(isort, numberRows);
3400        double * other = new double[numberRows];
3401        CoinZeroN(other, numberRows);
3402        int * which = new int[numberRows];
3403        //#define COUNT
3404#ifndef COUNT
3405        bool beforeSolution = model_->getSolutionCount() == 0;
3406#endif
3407        for (int k = 0; k < nSort - 1; k++) {
3408            i = sort[k];
3409            int numberUnsatisfied = 0;
3410            int n = 0;
3411            int j;
3412            for (j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
3413                int iColumn = column[j];
3414                if (columnLower[iColumn] != columnUpper[iColumn]) {
3415                    double solValue = solution[iColumn] - columnLower[iColumn];
3416                    if (solValue < 1.0 - integerTolerance && solValue > integerTolerance) {
3417                        numberUnsatisfied++;
3418                        for (int jj = columnStart[iColumn]; jj < columnStart[iColumn] + columnLength[iColumn]; jj++) {
3419                            int iRow = row[jj];
3420                            if (rhs_[iRow]) {
3421                                other[iRow] += solValue;
3422                                if (isort[iRow]) {
3423                                    isort[iRow]++;
3424                                } else {
3425                                    isort[iRow] = 1;
3426                                    which[n++] = iRow;
3427                                }
3428                            }
3429                        }
3430                    }
3431                }
3432            }
3433            double total = 0.0;
3434            // Take out row
3435            double sumThis = other[i];
3436            other[i] = 0.0;
3437            assert (numberUnsatisfied == isort[i]);
3438            // find one nearest half if solution, one if before solution
3439            int iBest = -1;
3440            double dtarget = 0.5 * total;
3441#ifdef COUNT
3442            int target = (numberUnsatisfied + 1) >> 1;
3443            int best = numberUnsatisfied;
3444#else
3445            double best;
3446            if (beforeSolution)
3447                best = dtarget;
3448            else
3449                best = 1.0e30;
3450#endif
3451            for (j = 0; j < n; j++) {
3452                int iRow = which[j];
3453                double dvalue = other[iRow];
3454                other[iRow] = 0.0;
3455#ifdef COUNT
3456                int value = isort[iRow];
3457#endif
3458                isort[iRow] = 0;
3459                if (fabs(dvalue) < 1.0e-8 || fabs(sumThis - dvalue) < 1.0e-8)
3460                    continue;
3461                if (dvalue < integerTolerance || dvalue > 1.0 - integerTolerance)
3462                    continue;
3463#ifdef COUNT
3464                if (abs(value - target) < best && value != numberUnsatisfied) {
3465                    best = abs(value - target);
3466                    iBest = iRow;
3467                    if (dvalue < dtarget)
3468                        preferredWay = 1;
3469                    else
3470                        preferredWay = -1;
3471                }
3472#else
3473                if (beforeSolution) {
3474                    if (fabs(dvalue - dtarget) > best) {
3475                        best = fabs(dvalue - dtarget);
3476                        iBest = iRow;
3477                        if (dvalue < dtarget)
3478                            preferredWay = 1;
3479                        else
3480                            preferredWay = -1;
3481                    }
3482                } else {
3483                    if (fabs(dvalue - dtarget) < best) {
3484                        best = fabs(dvalue - dtarget);
3485                        iBest = iRow;
3486                        if (dvalue < dtarget)
3487                            preferredWay = 1;
3488                        else
3489                            preferredWay = -1;
3490                    }
3491                }
3492#endif
3493            }
3494            if (iBest >= 0) {
3495                whichRow = i;
3496                otherRow = iBest;
3497                break;
3498            }
3499        }
3500        delete [] which;
3501        delete [] other;
3502    }
3503    delete [] sort;
3504    delete [] isort;
3505    return whichRow;
3506}
3507double
3508CbcFollowOn::infeasibility(const OsiBranchingInformation * /*info*/,
3509                           int &preferredWay) const
3510{
3511    int otherRow = 0;
3512    int whichRow = gutsOfFollowOn(otherRow, preferredWay);
3513    if (whichRow < 0)
3514        return 0.0;
3515    else
3516        return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
3517}
3518
3519// This looks at solution and sets bounds to contain solution
3520void
3521CbcFollowOn::feasibleRegion()
3522{
3523}
3524
3525CbcBranchingObject *
3526CbcFollowOn::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way)
3527{
3528    int otherRow = 0;
3529    int preferredWay;
3530    int whichRow = gutsOfFollowOn(otherRow, preferredWay);
3531    assert(way == preferredWay);
3532    assert (whichRow >= 0);
3533    int numberColumns = matrix_.getNumCols();
3534
3535    // Column copy
3536    //const double * element = matrix_.getElements();
3537    const int * row = matrix_.getIndices();
3538    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
3539    const int * columnLength = matrix_.getVectorLengths();
3540    // Row copy
3541    //const double * elementByRow = matrixByRow_.getElements();
3542    const int * column = matrixByRow_.getIndices();
3543    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
3544    const int * rowLength = matrixByRow_.getVectorLengths();
3545    //OsiSolverInterface * solver = model_->solver();
3546    const double * columnLower = solver->getColLower();
3547    const double * columnUpper = solver->getColUpper();
3548    //const double * solution = solver->getColSolution();
3549    int nUp = 0;
3550    int nDown = 0;
3551    int * upList = new int[numberColumns];
3552    int * downList = new int[numberColumns];
3553    int j;
3554    for (j = rowStart[whichRow]; j < rowStart[whichRow] + rowLength[whichRow]; j++) {
3555        int iColumn = column[j];
3556        if (columnLower[iColumn] != columnUpper[iColumn]) {
3557            bool up = true;
3558            for (int jj = columnStart[iColumn]; jj < columnStart[iColumn] + columnLength[iColumn]; jj++) {
3559                int iRow = row[jj];
3560                if (iRow == otherRow) {
3561                    up = false;
3562                    break;
3563                }
3564            }
3565            if (up)
3566                upList[nUp++] = iColumn;
3567            else
3568                downList[nDown++] = iColumn;
3569        }
3570    }
3571    //printf("way %d\n",way);
3572    // create object
3573    //printf("would fix %d down and %d up\n",nDown,nUp);
3574    CbcBranchingObject * branch
3575    = new CbcFixingBranchingObject(model_, way,
3576                                   nDown, downList, nUp, upList);
3577    delete [] upList;
3578    delete [] downList;
3579    return branch;
3580}
3581
3582//##############################################################################
3583
3584// Default Constructor
3585CbcFixingBranchingObject::CbcFixingBranchingObject()
3586        : CbcBranchingObject()
3587{
3588    numberDown_ = 0;
3589    numberUp_ = 0;
3590    downList_ = NULL;
3591    upList_ = NULL;
3592}
3593
3594// Useful constructor
3595CbcFixingBranchingObject::CbcFixingBranchingObject (CbcModel * model,
3596        int way ,
3597        int numberOnDownSide, const int * down,
3598        int numberOnUpSide, const int * up)
3599        : CbcBranchingObject(model, 0, way, 0.5)
3600{
3601    numberDown_ = numberOnDownSide;
3602    numberUp_ = numberOnUpSide;
3603    downList_ = CoinCopyOfArray(down, numberDown_);
3604    upList_ = CoinCopyOfArray(up, numberUp_);
3605}
3606
3607// Copy constructor
3608CbcFixingBranchingObject::CbcFixingBranchingObject ( const CbcFixingBranchingObject & rhs) : CbcBranchingObject(rhs)
3609{
3610    numberDown_ = rhs.numberDown_;
3611    numberUp_ = rhs.numberUp_;
3612    downList_ = CoinCopyOfArray(rhs.downList_, numberDown_);
3613    upList_ = CoinCopyOfArray(rhs.upList_, numberUp_);
3614}
3615
3616// Assignment operator
3617CbcFixingBranchingObject &
3618CbcFixingBranchingObject::operator=( const CbcFixingBranchingObject & rhs)
3619{
3620    if (this != &rhs) {
3621        CbcBranchingObject::operator=(rhs);
3622        delete [] downList_;
3623        delete [] upList_;
3624        numberDown_ = rhs.numberDown_;
3625        numberUp_ = rhs.numberUp_;
3626        downList_ = CoinCopyOfArray(rhs.downList_, numberDown_);
3627        upList_ = CoinCopyOfArray(rhs.upList_, numberUp_);
3628    }
3629    return *this;
3630}
3631CbcBranchingObject *
3632CbcFixingBranchingObject::clone() const
3633{
3634    return (new CbcFixingBranchingObject(*this));
3635}
3636
3637
3638// Destructor
3639CbcFixingBranchingObject::~CbcFixingBranchingObject ()
3640{
3641    delete [] downList_;
3642    delete [] upList_;
3643}
3644double
3645CbcFixingBranchingObject::branch()
3646{
3647    decrementNumberBranchesLeft();
3648    OsiSolverInterface * solver = model_->solver();
3649    const double * columnLower = solver->getColLower();
3650    int i;
3651    // *** for way - up means fix all those in up section
3652    if (way_ < 0) {
3653#ifdef FULL_PRINT
3654        printf("Down Fix ");
3655#endif
3656        //printf("Down Fix %d\n",numberDown_);
3657        for (i = 0; i < numberDown_; i++) {
3658            int iColumn = downList_[i];
3659            model_->solver()->setColUpper(iColumn, columnLower[iColumn]);
3660#ifdef FULL_PRINT
3661            printf("Setting bound on %d to lower bound\n", iColumn);
3662#endif
3663        }
3664        way_ = 1;         // Swap direction
3665    } else {
3666#ifdef FULL_PRINT
3667        printf("Up Fix ");
3668#endif
3669        //printf("Up Fix %d\n",numberUp_);
3670        for (i = 0; i < numberUp_; i++) {
3671            int iColumn = upList_[i];
3672            model_->solver()->setColUpper(iColumn, columnLower[iColumn]);
3673#ifdef FULL_PRINT
3674            printf("Setting bound on %d to lower bound\n", iColumn);
3675#endif
3676        }
3677        way_ = -1;        // Swap direction
3678    }
3679#ifdef FULL_PRINT
3680    printf("\n");
3681#endif
3682    return 0.0;
3683}
3684void
3685CbcFixingBranchingObject::print()
3686{
3687    int i;
3688    // *** for way - up means fix all those in up section
3689    if (way_ < 0) {
3690        printf("Down Fix ");
3691        for (i = 0; i < numberDown_; i++) {
3692            int iColumn = downList_[i];
3693            printf("%d ", iColumn);
3694        }
3695    } else {
3696        printf("Up Fix ");
3697        for (i = 0; i < numberUp_; i++) {
3698            int iColumn = upList_[i];
3699            printf("%d ", iColumn);
3700        }
3701    }
3702    printf("\n");
3703}
3704
3705/** Compare the original object of \c this with the original object of \c
3706    brObj. Assumes that there is an ordering of the original objects.
3707    This method should be invoked only if \c this and brObj are of the same
3708    type.
3709    Return negative/0/positive depending on whether \c this is
3710    smaller/same/larger than the argument.
3711*/
3712int
3713CbcFixingBranchingObject::compareOriginalObject
3714(const CbcBranchingObject* /*brObj*/) const
3715{
3716    throw("must implement");
3717}
3718
3719/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
3720    same type and must have the same original object, but they may have
3721    different feasible regions.
3722    Return the appropriate CbcRangeCompare value (first argument being the
3723    sub/superset if that's the case). In case of overlap (and if \c
3724    replaceIfOverlap is true) replace the current branching object with one
3725    whose feasible region is the overlap.
3726   */
3727CbcRangeCompare
3728CbcFixingBranchingObject::compareBranchingObject
3729(const CbcBranchingObject* /*brObj*/, const bool /*replaceIfOverlap*/)
3730{
3731#if 0 //ndef NDEBUG
3732    const CbcFixingBranchingObject* br =
3733        dynamic_cast<const CbcFixingBranchingObject*>(brObj);
3734    assert(br);
3735#endif
3736    // If two FixingBranchingObject's have the same base object then it's pretty
3737    // much guaranteed
3738    throw("must implement");
3739}
3740
3741//##############################################################################
3742
3743// Default Constructor
3744CbcNWay::CbcNWay ()
3745        : CbcObject(),
3746        numberMembers_(0),
3747        members_(NULL),
3748        consequence_(NULL)
3749{
3750}
3751
3752// Useful constructor (which are integer indices)
3753CbcNWay::CbcNWay (CbcModel * model, int numberMembers,
3754                  const int * which, int identifier)
3755        : CbcObject(model)
3756{
3757    id_ = identifier;
3758    numberMembers_ = numberMembers;
3759    if (numberMembers_) {
3760        members_ = new int[numberMembers_];
3761        memcpy(members_, which, numberMembers_*sizeof(int));
3762    } else {
3763        members_ = NULL;
3764    }
3765    consequence_ = NULL;
3766}
3767
3768// Copy constructor
3769CbcNWay::CbcNWay ( const CbcNWay & rhs)
3770        : CbcObject(rhs)
3771{
3772    numberMembers_ = rhs.numberMembers_;
3773    consequence_ = NULL;
3774    if (numberMembers_) {
3775        members_ = new int[numberMembers_];
3776        memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
3777        if (rhs.consequence_) {
3778            consequence_ = new CbcConsequence * [numberMembers_];
3779            for (int i = 0; i < numberMembers_; i++) {
3780                if (rhs.consequence_[i])
3781                    consequence_[i] = rhs.consequence_[i]->clone();
3782                else
3783                    consequence_[i] = NULL;
3784            }
3785        }
3786    } else {
3787        members_ = NULL;
3788    }
3789}
3790
3791// Clone
3792CbcObject *
3793CbcNWay::clone() const
3794{
3795    return new CbcNWay(*this);
3796}
3797
3798// Assignment operator
3799CbcNWay &
3800CbcNWay::operator=( const CbcNWay & rhs)
3801{
3802    if (this != &rhs) {
3803        CbcObject::operator=(rhs);
3804        delete [] members_;
3805        numberMembers_ = rhs.numberMembers_;
3806        if (consequence_) {
3807            for (int i = 0; i < numberMembers_; i++)
3808                delete consequence_[i];
3809            delete [] consequence_;
3810            consequence_ = NULL;
3811        }
3812        if (numberMembers_) {
3813            members_ = new int[numberMembers_];
3814            memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
3815        } else {
3816            members_ = NULL;
3817        }
3818        if (rhs.consequence_) {
3819            consequence_ = new CbcConsequence * [numberMembers_];
3820            for (int i = 0; i < numberMembers_; i++) {
3821                if (rhs.consequence_[i])
3822                    consequence_[i] = rhs.consequence_[i]->clone();
3823                else
3824                    consequence_[i] = NULL;
3825            }
3826        }
3827    }
3828    return *this;
3829}
3830
3831// Destructor
3832CbcNWay::~CbcNWay ()
3833{
3834    delete [] members_;
3835    if (consequence_) {
3836        for (int i = 0; i < numberMembers_; i++)
3837            delete consequence_[i];
3838        delete [] consequence_;
3839    }
3840}
3841// Set up a consequence for a single member
3842void
3843CbcNWay::setConsequence(int iColumn, const CbcConsequence & consequence)
3844{
3845    if (!consequence_) {
3846        consequence_ = new CbcConsequence * [numberMembers_];
3847        for (int i = 0; i < numberMembers_; i++)
3848            consequence_[i] = NULL;
3849    }
3850    for (int i = 0; i < numberMembers_; i++) {
3851        if (members_[i] == iColumn) {
3852            consequence_[i] = consequence.clone();
3853            break;
3854        }
3855    }
3856}
3857
3858// Applies a consequence for a single member
3859void
3860CbcNWay::applyConsequence(int iSequence, int state) const
3861{
3862    assert (state == -9999 || state == 9999);
3863    if (consequence_) {
3864        CbcConsequence * consequence = consequence_[iSequence];
3865        if (consequence)
3866            consequence->applyToSolver(model_->solver(), state);
3867    }
3868}
3869double
3870CbcNWay::infeasibility(const OsiBranchingInformation * /*info*/,
3871                       int &preferredWay) const
3872{
3873    int numberUnsatis = 0;
3874    int j;
3875    OsiSolverInterface * solver = model_->solver();
3876    const double * solution = model_->testSolution();
3877    const double * lower = solver->getColLower();
3878    const double * upper = solver->getColUpper();
3879    double largestValue = 0.0;
3880
3881    double integerTolerance =
3882        model_->getDblParam(CbcModel::CbcIntegerTolerance);
3883
3884    for (j = 0; j < numberMembers_; j++) {
3885        int iColumn = members_[j];
3886        double value = solution[iColumn];
3887        value = CoinMax(value, lower[iColumn]);
3888        value = CoinMin(value, upper[iColumn]);
3889        double distance = CoinMin(value - lower[iColumn], upper[iColumn] - value);
3890        if (distance > integerTolerance) {
3891            numberUnsatis++;
3892            largestValue = CoinMax(distance, largestValue);
3893        }
3894    }
3895    preferredWay = 1;
3896    if (numberUnsatis) {
3897        return largestValue;
3898    } else {
3899        return 0.0; // satisfied
3900    }
3901}
3902
3903// This looks at solution and sets bounds to contain solution
3904void
3905CbcNWay::feasibleRegion()
3906{
3907    int j;
3908    OsiSolverInterface * solver = model_->solver();
3909    const double * solution = model_->testSolution();
3910    const double * lower = solver->getColLower();
3911    const double * upper = solver->getColUpper();
3912    double integerTolerance =
3913        model_->getDblParam(CbcModel::CbcIntegerTolerance);
3914    for (j = 0; j < numberMembers_; j++) {
3915        int iColumn = members_[j];
3916        double value = solution[iColumn];
3917        value = CoinMax(value, lower[iColumn]);
3918        value = CoinMin(value, upper[iColumn]);
3919        if (value >= upper[iColumn] - integerTolerance) {
3920            solver->setColLower(iColumn, upper[iColumn]);
3921        } else {
3922            assert (value <= lower[iColumn] + integerTolerance);
3923            solver->setColUpper(iColumn, lower[iColumn]);
3924        }
3925    }
3926}
3927// Redoes data when sequence numbers change
3928void
3929CbcNWay::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
3930{
3931    model_ = model;
3932    int n2 = 0;
3933    for (int j = 0; j < numberMembers_; j++) {
3934        int iColumn = members_[j];
3935        int i;
3936        for (i = 0; i < numberColumns; i++) {
3937            if (originalColumns[i] == iColumn)
3938                break;
3939        }
3940        if (i < numberColumns) {
3941            members_[n2] = i;
3942            consequence_[n2++] = consequence_[j];
3943        } else {
3944            delete consequence_[j];
3945        }
3946    }
3947    if (n2 < numberMembers_) {
3948        printf("** NWay number of members reduced from %d to %d!\n", numberMembers_, n2);
3949        numberMembers_ = n2;
3950    }
3951}
3952CbcBranchingObject *
3953CbcNWay::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int /*way*/)
3954{
3955    int numberFree = 0;
3956    int j;
3957
3958    //OsiSolverInterface * solver = model_->solver();
3959    const double * solution = model_->testSolution();
3960    const double * lower = solver->getColLower();
3961    const double * upper = solver->getColUpper();
3962    int * list = new int[numberMembers_];
3963    double * sort = new double[numberMembers_];
3964
3965    for (j = 0; j < numberMembers_; j++) {
3966        int iColumn = members_[j];
3967        double value = solution[iColumn];
3968        value = CoinMax(value, lower[iColumn]);
3969        value = CoinMin(value, upper[iColumn]);
3970        if (upper[iColumn] > lower[iColumn]) {
3971            double distance = upper[iColumn] - value;
3972            list[numberFree] = j;
3973            sort[numberFree++] = distance;
3974        }
3975    }
3976    assert (numberFree);
3977    // sort
3978    CoinSort_2(sort, sort + numberFree, list);
3979    // create object
3980    CbcBranchingObject * branch;
3981    branch = new CbcNWayBranchingObject(model_, this, numberFree, list);
3982    branch->setOriginalObject(this);
3983    delete [] list;
3984    delete [] sort;
3985    return branch;
3986}
3987
3988//##############################################################################
3989
3990// Default Constructor
3991CbcNWayBranchingObject::CbcNWayBranchingObject()
3992        : CbcBranchingObject()
3993{
3994    order_ = NULL;
3995    object_ = NULL;
3996    numberInSet_ = 0;
3997    way_ = 0;
3998}
3999
4000// Useful constructor
4001CbcNWayBranchingObject::CbcNWayBranchingObject (CbcModel * model,
4002        const CbcNWay * nway,
4003        int number, const int * order)
4004        : CbcBranchingObject(model, nway->id(), -1, 0.5)
4005{
4006    numberBranches_ = number;
4007    order_ = new int [number];
4008    object_ = nway;
4009    numberInSet_ = number;
4010    memcpy(order_, order, number*sizeof(int));
4011}
4012
4013// Copy constructor
4014CbcNWayBranchingObject::CbcNWayBranchingObject ( const CbcNWayBranchingObject & rhs) : CbcBranchingObject(rhs)
4015{
4016    numberInSet_ = rhs.numberInSet_;
4017    object_ = rhs.object_;
4018    if (numberInSet_) {
4019        order_ = new int [numberInSet_];
4020        memcpy(order_, rhs.order_, numberInSet_*sizeof(int));
4021    } else {
4022        order_ = NULL;
4023    }
4024}
4025
4026// Assignment operator
4027CbcNWayBranchingObject &
4028CbcNWayBranchingObject::operator=( const CbcNWayBranchingObject & rhs)
4029{
4030    if (this != &rhs) {
4031        CbcBranchingObject::operator=(rhs);
4032        object_ = rhs.object_;
4033        delete [] order_;
4034        numberInSet_ = rhs.numberInSet_;
4035        if (numberInSet_) {
4036            order_ = new int [numberInSet_];
4037            memcpy(order_, rhs.order_, numberInSet_*sizeof(int));
4038        } else {
4039            order_ = NULL;
4040        }
4041    }
4042    return *this;
4043}
4044CbcBranchingObject *
4045CbcNWayBranchingObject::clone() const
4046{
4047    return (new CbcNWayBranchingObject(*this));
4048}
4049
4050
4051// Destructor
4052CbcNWayBranchingObject::~CbcNWayBranchingObject ()
4053{
4054    delete [] order_;
4055}
4056double
4057CbcNWayBranchingObject::branch()
4058{
4059    int which = branchIndex_;
4060    branchIndex_++;
4061    assert (numberBranchesLeft() >= 0);
4062    if (which == 0) {
4063        // first branch so way_ may mean something
4064        assert (way_ == -1 || way_ == 1);
4065        if (way_ == -1)
4066            which++;
4067    } else if (which == 1) {
4068        // second branch so way_ may mean something
4069        assert (way_ == -1 || way_ == 1);
4070        if (way_ == -1)
4071            which--;
4072        // switch way off
4073        way_ = 0;
4074    }
4075    const double * lower = model_->solver()->getColLower();
4076    const double * upper = model_->solver()->getColUpper();
4077    const int * members = object_->members();
4078    for (int j = 0; j < numberInSet_; j++) {
4079        int iSequence = order_[j];
4080        int iColumn = members[iSequence];
4081        if (j != which) {
4082            model_->solver()->setColUpper(iColumn, lower[iColumn]);
4083            //model_->solver()->setColLower(iColumn,lower[iColumn]);
4084            assert (lower[iColumn] > -1.0e20);
4085            // apply any consequences
4086            object_->applyConsequence(iSequence, -9999);
4087        } else {
4088            model_->solver()->setColLower(iColumn, upper[iColumn]);
4089            //model_->solver()->setColUpper(iColumn,upper[iColumn]);
4090#ifdef FULL_PRINT
4091            printf("Up Fix %d to %g\n", iColumn, upper[iColumn]);
4092#endif
4093            assert (upper[iColumn] < 1.0e20);
4094            // apply any consequences
4095            object_->applyConsequence(iSequence, 9999);
4096        }
4097    }
4098    return 0.0;
4099}
4100void
4101CbcNWayBranchingObject::print()
4102{
4103    printf("NWay - Up Fix ");
4104    const int * members = object_->members();
4105    for (int j = 0; j < way_; j++) {
4106        int iColumn = members[order_[j]];
4107        printf("%d ", iColumn);
4108    }
4109    printf("\n");
4110}
4111
4112/** Compare the original object of \c this with the original object of \c
4113    brObj. Assumes that there is an ordering of the original objects.
4114    This method should be invoked only if \c this and brObj are of the same
4115    type.
4116    Return negative/0/positive depending on whether \c this is
4117    smaller/same/larger than the argument.
4118*/
4119int
4120CbcNWayBranchingObject::compareOriginalObject
4121(const CbcBranchingObject* /*brObj*/) const
4122{
4123    throw("must implement");
4124}
4125
4126/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
4127    same type and must have the same original object, but they may have
4128    different feasible regions.
4129    Return the appropriate CbcRangeCompare value (first argument being the
4130    sub/superset if that's the case). In case of overlap (and if \c
4131    replaceIfOverlap is true) replace the current branching object with one
4132    whose feasible region is the overlap.
4133*/
4134CbcRangeCompare
4135CbcNWayBranchingObject::compareBranchingObject
4136(const CbcBranchingObject* /*brObj*/, const bool /*replaceIfOverlap*/)
4137{
4138    throw("must implement");
4139}
4140
4141//##############################################################################
4142
4143// Default Constructor
4144CbcFixVariable::CbcFixVariable ()
4145        : CbcConsequence(),
4146        numberStates_(0),
4147        states_(NULL),
4148        startLower_(NULL),
4149        startUpper_(NULL),
4150        newBound_(NULL),
4151        variable_(NULL)
4152{
4153}
4154
4155// One useful Constructor
4156CbcFixVariable::CbcFixVariable (int numberStates, const int * states, const int * numberNewLower,
4157                                const int ** newLowerValue,
4158                                const int ** lowerColumn,
4159                                const int * numberNewUpper, const int ** newUpperValue,
4160                                const int ** upperColumn)
4161        : CbcConsequence(),
4162        states_(NULL),
4163        startLower_(NULL),
4164        startUpper_(NULL),
4165        newBound_(NULL),
4166        variable_(NULL)
4167{
4168    // How much space
4169    numberStates_ = numberStates;
4170    if (numberStates_) {
4171        states_ = new int[numberStates_];
4172        memcpy(states_, states, numberStates_*sizeof(int));
4173        int i;
4174        int n = 0;
4175        startLower_ = new int[numberStates_+1];
4176        startUpper_ = new int[numberStates_+1];
4177        startLower_[0] = 0;
4178        //count
4179        for (i = 0; i < numberStates_; i++) {
4180            n += numberNewLower[i];
4181            startUpper_[i] = n;
4182            n += numberNewUpper[i];
4183            startLower_[i+1] = n;
4184        }
4185        newBound_ = new double [n];
4186        variable_ = new int [n];
4187        n = 0;
4188        for (i = 0; i < numberStates_; i++) {
4189            int j;
4190            int k;
4191            const int * bound;
4192            const int * variable;
4193            k = numberNewLower[i];
4194            bound = newLowerValue[i];
4195            variable = lowerColumn[i];
4196            for (j = 0; j < k; j++) {
4197                newBound_[n] = bound[j];
4198                variable_[n++] = variable[j];
4199            }
4200            k = numberNewUpper[i];
4201            bound = newUpperValue[i];
4202            variable = upperColumn[i];
4203            for (j = 0; j < k; j++) {
4204                newBound_[n] = bound[j];
4205                variable_[n++] = variable[j];
4206            }
4207        }
4208    }
4209}
4210
4211// Copy constructor
4212CbcFixVariable::CbcFixVariable ( const CbcFixVariable & rhs)
4213        : CbcConsequence(rhs)
4214{
4215    numberStates_ = rhs.numberStates_;
4216    states_ = NULL;
4217    startLower_ = NULL;
4218    startUpper_ = NULL;
4219    newBound_ = NULL;
4220    variable_ = NULL;
4221    if (numberStates_) {
4222        states_ = CoinCopyOfArray(rhs.states_, numberStates_);
4223        startLower_ = CoinCopyOfArray(rhs.startLower_, numberStates_ + 1);
4224        startUpper_ = CoinCopyOfArray(rhs.startUpper_, numberStates_ + 1);
4225        int n = startLower_[numberStates_];
4226        newBound_ = CoinCopyOfArray(rhs.newBound_, n);
4227        variable_ = CoinCopyOfArray(rhs.variable_, n);
4228    }
4229}
4230
4231// Clone
4232CbcConsequence *
4233CbcFixVariable::clone() const
4234{
4235    return new CbcFixVariable(*this);
4236}
4237
4238// Assignment operator
4239CbcFixVariable &
4240CbcFixVariable::operator=( const CbcFixVariable & rhs)
4241{
4242    if (this != &rhs) {
4243        CbcConsequence::operator=(rhs);
4244        delete [] states_;
4245        delete [] startLower_;
4246        delete [] startUpper_;
4247        delete [] newBound_;
4248        delete [] variable_;
4249        states_ = NULL;
4250        startLower_ = NULL;
4251        startUpper_ = NULL;
4252        newBound_ = NULL;
4253        variable_ = NULL;
4254        numberStates_ = rhs.numberStates_;
4255        if (numberStates_) {
4256            states_ = CoinCopyOfArray(rhs.states_, numberStates_);
4257            startLower_ = CoinCopyOfArray(rhs.startLower_, numberStates_ + 1);
4258            startUpper_ = CoinCopyOfArray(rhs.startUpper_, numberStates_ + 1);
4259            int n = startLower_[numberStates_];
4260            newBound_ = CoinCopyOfArray(rhs.newBound_, n);
4261            variable_ = CoinCopyOfArray(rhs.variable_, n);
4262        }
4263    }
4264    return *this;
4265}
4266
4267// Destructor
4268CbcFixVariable::~CbcFixVariable ()
4269{
4270    delete [] states_;
4271    delete [] startLower_;
4272    delete [] startUpper_;
4273    delete [] newBound_;
4274    delete [] variable_;
4275}
4276// Set up a startLower for a single member
4277void
4278CbcFixVariable::applyToSolver(OsiSolverInterface * solver, int state) const
4279{
4280    assert (state == -9999 || state == 9999);
4281    // Find state
4282    int find;
4283    for (find = 0; find < numberStates_; find++)
4284        if (states_[find] == state)
4285            break;
4286    if (find == numberStates_)
4287        return;
4288    int i;
4289    // Set new lower bounds
4290    for (i = startLower_[find]; i < startUpper_[find]; i++) {
4291        int iColumn = variable_[i];
4292        double value = newBound_[i];
4293        double oldValue = solver->getColLower()[iColumn];
4294        //printf("for %d old lower bound %g, new %g",iColumn,oldValue,value);
4295        solver->setColLower(iColumn, CoinMax(value, oldValue));
4296        //printf(" => %g\n",solver->getColLower()[iColumn]);
4297    }
4298    // Set new upper bounds
4299    for (i = startUpper_[find]; i < startLower_[find+1]; i++) {
4300        int iColumn = variable_[i];
4301        double value = newBound_[i];
4302        double oldValue = solver->getColUpper()[iColumn];
4303        //printf("for %d old upper bound %g, new %g",iColumn,oldValue,value);
4304        solver->setColUpper(iColumn, CoinMin(value, oldValue));
4305        //printf(" => %g\n",solver->getColUpper()[iColumn]);
4306    }
4307}
4308
4309//##############################################################################
4310
4311// Default Constructor
4312CbcDummyBranchingObject::CbcDummyBranchingObject(CbcModel * model)
4313        : CbcBranchingObject(model, 0, 0, 0.5)
4314{
4315    setNumberBranchesLeft(1);
4316}
4317
4318
4319// Copy constructor
4320CbcDummyBranchingObject::CbcDummyBranchingObject ( const CbcDummyBranchingObject & rhs) : CbcBranchingObject(rhs)
4321{
4322}
4323
4324// Assignment operator
4325CbcDummyBranchingObject &
4326CbcDummyBranchingObject::operator=( const CbcDummyBranchingObject & rhs)
4327{
4328    if (this != &rhs) {
4329        CbcBranchingObject::operator=(rhs);
4330    }
4331    return *this;
4332}
4333CbcBranchingObject *
4334CbcDummyBranchingObject::clone() const
4335{
4336    return (new CbcDummyBranchingObject(*this));
4337}
4338
4339
4340// Destructor
4341CbcDummyBranchingObject::~CbcDummyBranchingObject ()
4342{
4343}
4344
4345/*
4346  Perform a dummy branch
4347*/
4348double
4349CbcDummyBranchingObject::branch()
4350{
4351    decrementNumberBranchesLeft();
4352    return 0.0;
4353}
4354// Print what would happen
4355void
4356CbcDummyBranchingObject::print()
4357{
4358    printf("Dummy branch\n");
4359}
4360
4361// Default Constructor
4362CbcGeneral::CbcGeneral()
4363        : CbcObject()
4364{
4365}
4366
4367// Constructor from model
4368CbcGeneral::CbcGeneral(CbcModel * model)
4369        : CbcObject(model)
4370{
4371}
4372
4373
4374// Destructor
4375CbcGeneral::~CbcGeneral ()
4376{
4377}
4378
4379// Copy constructor
4380CbcGeneral::CbcGeneral ( const CbcGeneral & rhs)
4381        : CbcObject(rhs)
4382{
4383}
4384#ifdef COIN_HAS_CLP
4385#include "OsiClpSolverInterface.hpp"
4386#include "CoinWarmStartBasis.hpp"
4387#include "ClpNode.hpp"
4388#include "CbcBranchDynamic.hpp"
4389// Assignment operator
4390CbcGeneral &
4391CbcGeneral::operator=( const CbcGeneral & rhs)
4392{
4393    if (this != &rhs) {
4394        CbcObject::operator=(rhs);
4395    }
4396    return *this;
4397}
4398// Infeasibility - large is 0.5
4399double
4400CbcGeneral::infeasibility(const OsiBranchingInformation * /*info*/,
4401                          int &/*preferredWay*/) const
4402{
4403    abort();
4404    return 0.0;
4405}
4406CbcBranchingObject *
4407CbcGeneral::createCbcBranch(OsiSolverInterface * /*solver*/, const OsiBranchingInformation * /*info*/, int /*way*/)
4408{
4409    abort();
4410    return NULL;
4411}
4412
4413// Default Constructor
4414CbcGeneralDepth::CbcGeneralDepth ()
4415        : CbcGeneral(),
4416        maximumDepth_(0),
4417        maximumNodes_(0),
4418        whichSolution_(-1),
4419        numberNodes_(0),
4420        nodeInfo_(NULL)
4421{
4422}
4423
4424// Useful constructor (which are integer indices)
4425CbcGeneralDepth::CbcGeneralDepth (CbcModel * model, int maximumDepth)
4426        : CbcGeneral(model),
4427        maximumDepth_(maximumDepth),
4428        maximumNodes_(0),
4429        whichSolution_(-1),
4430        numberNodes_(0),
4431        nodeInfo_(NULL)
4432{
4433    assert(maximumDepth_ < 1000000);
4434    if (maximumDepth_ > 0)
4435        maximumNodes_ = (1 << maximumDepth_) + 1 + maximumDepth_;
4436    else if (maximumDepth_ < 0)
4437        maximumNodes_ = 1 + 1 - maximumDepth_;
4438    else
4439        maximumNodes_ = 0;
4440#define MAX_NODES 100
4441    maximumNodes_ = CoinMin(maximumNodes_, 1 + maximumDepth_ + MAX_NODES);
4442    if (maximumNodes_) {
4443        nodeInfo_ = new ClpNodeStuff();
4444        nodeInfo_->maximumNodes_ = maximumNodes_;
4445        ClpNodeStuff * info = nodeInfo_;
4446        // for reduced costs and duals
4447        info->solverOptions_ |= 7;
4448        if (maximumDepth_ > 0) {
4449            info->nDepth_ = maximumDepth_;
4450        } else {
4451            info->nDepth_ = - maximumDepth_;
4452            info->solverOptions_ |= 32;
4453        }
4454        ClpNode ** nodeInfo = new ClpNode * [maximumNodes_];
4455        for (int i = 0; i < maximumNodes_; i++)
4456            nodeInfo[i] = NULL;
4457        info->nodeInfo_ = nodeInfo;
4458    } else {
4459        nodeInfo_ = NULL;
4460    }
4461}
4462
4463// Copy constructor
4464CbcGeneralDepth::CbcGeneralDepth ( const CbcGeneralDepth & rhs)
4465        : CbcGeneral(rhs)
4466{
4467    maximumDepth_ = rhs.maximumDepth_;
4468    maximumNodes_ = rhs.maximumNodes_;
4469    whichSolution_ = -1;
4470    numberNodes_ = 0;
4471    if (maximumNodes_) {
4472        assert (rhs.nodeInfo_);
4473        nodeInfo_ = new ClpNodeStuff(*rhs.nodeInfo_);
4474        nodeInfo_->maximumNodes_ = maximumNodes_;
4475        ClpNodeStuff * info = nodeInfo_;
4476        if (maximumDepth_ > 0) {
4477            info->nDepth_ = maximumDepth_;
4478        } else {
4479            info->nDepth_ = - maximumDepth_;
4480            info->solverOptions_ |= 32;
4481        }
4482        if (!info->nodeInfo_) {
4483            ClpNode ** nodeInfo = new ClpNode * [maximumNodes_];
4484            for (int i = 0; i < maximumNodes_; i++)
4485                nodeInfo[i] = NULL;
4486            info->nodeInfo_ = nodeInfo;
4487        }
4488    } else {
4489        nodeInfo_ = NULL;
4490    }
4491}
4492
4493// Clone
4494CbcObject *
4495CbcGeneralDepth::clone() const
4496{
4497    return new CbcGeneralDepth(*this);
4498}
4499
4500// Assignment operator
4501CbcGeneralDepth &
4502CbcGeneralDepth::operator=( const CbcGeneralDepth & rhs)
4503{
4504    if (this != &rhs) {
4505        CbcGeneral::operator=(rhs);
4506        delete nodeInfo_;
4507        maximumDepth_ = rhs.maximumDepth_;
4508        maximumNodes_ = rhs.maximumNodes_;
4509        whichSolution_ = -1;
4510        numberNodes_ = 0;
4511        if (maximumDepth_) {
4512            assert (rhs.nodeInfo_);
4513            nodeInfo_ = new ClpNodeStuff(*rhs.nodeInfo_);
4514            nodeInfo_->maximumNodes_ = maximumNodes_;
4515        } else {
4516            nodeInfo_ = NULL;
4517        }
4518    }
4519    return *this;
4520}
4521
4522// Destructor
4523CbcGeneralDepth::~CbcGeneralDepth ()
4524{
4525    delete nodeInfo_;
4526}
4527// Infeasibility - large is 0.5
4528double
4529CbcGeneralDepth::infeasibility(const OsiBranchingInformation * /*info*/,
4530                               int &/*preferredWay*/) const
4531{
4532    whichSolution_ = -1;
4533    // should use genuine OsiBranchingInformation usefulInfo = model_->usefulInformation();
4534    // for now assume only called when correct
4535    //if (usefulInfo.depth_>=4&&!model_->parentModel()
4536    //     &&(usefulInfo.depth_%2)==0) {
4537    if (true) {
4538        OsiSolverInterface * solver = model_->solver();
4539        OsiClpSolverInterface * clpSolver
4540        = dynamic_cast<OsiClpSolverInterface *> (solver);
4541        if (clpSolver) {
4542            ClpNodeStuff * info = nodeInfo_;
4543            info->integerTolerance_ = model_->getIntegerTolerance();
4544            info->integerIncrement_ = model_->getCutoffIncrement();
4545            info->numberBeforeTrust_ = model_->numberBeforeTrust();
4546            info->stateOfSearch_ = model_->stateOfSearch();
4547            // Compute "small" change in branch
4548            int nBranches = model_->getIntParam(CbcModel::CbcNumberBranches);
4549            if (nBranches) {
4550                double average = model_->getDblParam(CbcModel::CbcSumChange) / static_cast<double>(nBranches);
4551                info->smallChange_ =
4552                    CoinMax(average * 1.0e-5, model_->getDblParam(CbcModel::CbcSmallestChange));
4553                info->smallChange_ = CoinMax(info->smallChange_, 1.0e-8);
4554            } else {
4555                info->smallChange_ = 1.0e-8;
4556            }
4557            int numberIntegers = model_->numberIntegers();
4558            double * down = new double[numberIntegers];
4559            double * up = new double[numberIntegers];
4560            int * priority = new int[numberIntegers];
4561            int * numberDown = new int[numberIntegers];
4562            int * numberUp = new int[numberIntegers];
4563            int * numberDownInfeasible = new int[numberIntegers];
4564            int * numberUpInfeasible = new int[numberIntegers];
4565            model_->fillPseudoCosts(down, up, priority, numberDown, numberUp,
4566                                    numberDownInfeasible, numberUpInfeasible);
4567            info->fillPseudoCosts(down, up, priority, numberDown, numberUp,
4568                                  numberDownInfeasible,
4569                                  numberUpInfeasible, numberIntegers);
4570            info->presolveType_ = 1;
4571            delete [] down;
4572            delete [] up;
4573            delete [] numberDown;
4574            delete [] numberUp;
4575            delete [] numberDownInfeasible;
4576            delete [] numberUpInfeasible;
4577            bool takeHint;
4578            OsiHintStrength strength;
4579            solver->getHintParam(OsiDoReducePrint, takeHint, strength);
4580            ClpSimplex * simplex = clpSolver->getModelPtr();
4581            int saveLevel = simplex->logLevel();
4582            if (strength != OsiHintIgnore && takeHint && saveLevel == 1)
4583                simplex->setLogLevel(0);
4584            clpSolver->setBasis();
4585            whichSolution_ = simplex->fathomMany(info);
4586            //printf("FAT %d nodes, %d iterations\n",
4587            //info->numberNodesExplored_,info->numberIterations_);
4588            //printf("CbcBranch %d rows, %d columns\n",clpSolver->getNumRows(),
4589            //     clpSolver->getNumCols());
4590            model_->incrementExtra(info->numberNodesExplored_,
4591                                   info->numberIterations_);
4592            // update pseudo costs
4593            double smallest = 1.0e50;
4594            double largest = -1.0;
4595            OsiObject ** objects = model_->objects();
4596#ifndef NDEBUG
4597            const int * integerVariable = model_->integerVariable();
4598#endif
4599            for (int i = 0; i < numberIntegers; i++) {
4600#ifndef NDEBUG
4601                CbcSimpleIntegerDynamicPseudoCost * obj =
4602                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(objects[i]) ;
4603                assert (obj && obj->columnNumber() == integerVariable[i]);
4604#else
4605                CbcSimpleIntegerDynamicPseudoCost * obj =
4606                    static_cast <CbcSimpleIntegerDynamicPseudoCost *>(objects[i]) ;
4607#endif
4608                if (info->numberUp_[i] > 0) {
4609                    if (info->downPseudo_[i] > largest)
4610                        largest = info->downPseudo_[i];
4611                    if (info->downPseudo_[i] < smallest)
4612                        smallest = info->downPseudo_[i];
4613                    if (info->upPseudo_[i] > largest)
4614                        largest = info->upPseudo_[i];
4615                    if (info->upPseudo_[i] < smallest)
4616                        smallest = info->upPseudo_[i];
4617                    obj->updateAfterMini(info->numberDown_[i],
4618                                         info->numberDownInfeasible_[i],
4619                                         info->downPseudo_[i],
4620                                         info->numberUp_[i],
4621                                         info->numberUpInfeasible_[i],
4622                                         info->upPseudo_[i]);
4623                }
4624            }
4625            //printf("range of costs %g to %g\n",smallest,largest);
4626            simplex->setLogLevel(saveLevel);
4627            numberNodes_ = info->nNodes_;
4628            int numberDo = numberNodes_;
4629            if (whichSolution_ >= 0)
4630                numberDo--;
4631            if (numberDo > 0) {
4632                return 0.5;
4633            } else {
4634                // no solution
4635                return COIN_DBL_MAX; // say infeasible
4636            }
4637        } else {
4638            return -1.0;
4639        }
4640    } else {
4641        return -1.0;
4642    }
4643}
4644
4645// This looks at solution and sets bounds to contain solution
4646void
4647CbcGeneralDepth::feasibleRegion()
4648{
4649    // Other stuff should have done this
4650}
4651// Redoes data when sequence numbers change
4652void
4653CbcGeneralDepth::redoSequenceEtc(CbcModel * /*model*/,
4654                                 int /*numberColumns*/,
4655                                 const int * /*originalColumns*/)
4656{
4657}
4658
4659//#define CHECK_PATH
4660#ifdef CHECK_PATH
4661extern const double * debuggerSolution_Z;
4662extern int numberColumns_Z;
4663extern int gotGoodNode_Z;
4664#endif
4665CbcBranchingObject *
4666CbcGeneralDepth::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int /*way*/)
4667{
4668    int numberDo = numberNodes_;
4669    if (whichSolution_ >= 0)
4670        numberDo--;
4671    assert (numberDo > 0);
4672    // create object
4673    CbcGeneralBranchingObject * branch = new CbcGeneralBranchingObject(model_);
4674    // skip solution
4675    branch->numberSubProblems_ = numberDo;
4676    // If parentBranch_ back in then will have to be 2*
4677    branch->numberSubLeft_ = numberDo;
4678    branch->setNumberBranches(numberDo);
4679    CbcSubProblem * sub = new CbcSubProblem[numberDo];
4680    int iProb = 0;
4681    branch->subProblems_ = sub;
4682    branch->numberRows_ = model_->solver()->getNumRows();
4683    int iNode;
4684    //OsiSolverInterface * solver = model_->solver();
4685    OsiClpSolverInterface * clpSolver
4686    = dynamic_cast<OsiClpSolverInterface *> (solver);
4687    assert (clpSolver);
4688    ClpSimplex * simplex = clpSolver->getModelPtr();
4689    int numberColumns = simplex->numberColumns();
4690    double * lowerBefore = CoinCopyOfArray(simplex->getColLower(),
4691                                           numberColumns);
4692    double * upperBefore = CoinCopyOfArray(simplex->getColUpper(),
4693                                           numberColumns);
4694    ClpNodeStuff * info = nodeInfo_;
4695    double * weight = new double[numberNodes_];
4696    int * whichNode = new int [numberNodes_];
4697    // Sort
4698    for (iNode = 0; iNode < numberNodes_; iNode++) {
4699        if (iNode != whichSolution_) {
4700            double objectiveValue = info->nodeInfo_[iNode]->objectiveValue();
4701            double sumInfeasibilities = info->nodeInfo_[iNode]->sumInfeasibilities();
4702            int numberInfeasibilities = info->nodeInfo_[iNode]->numberInfeasibilities();
4703            double thisWeight = 0.0;
4704#if 1
4705            // just closest
4706            thisWeight = 1.0e9 * numberInfeasibilities;
4707            thisWeight += sumInfeasibilities;
4708            thisWeight += 1.0e-7 * objectiveValue;
4709            // Try estimate
4710            thisWeight = info->nodeInfo_[iNode]->estimatedSolution();
4711#else
4712            thisWeight = 1.0e-3 * numberInfeasibilities;
4713            thisWeight += 1.0e-5 * sumInfeasibilities;
4714            thisWeight += objectiveValue;
4715#endif
4716            whichNode[iProb] = iNode;
4717            weight[iProb++] = thisWeight;
4718        }
4719    }
4720    assert (iProb == numberDo);
4721    CoinSort_2(weight, weight + numberDo, whichNode);
4722    for (iProb = 0; iProb < numberDo; iProb++) {
4723        iNode = whichNode[iProb];
4724        ClpNode * node = info->nodeInfo_[iNode];
4725        // move bounds
4726        node->applyNode(simplex, 3);
4727        // create subproblem
4728        sub[iProb] = CbcSubProblem(clpSolver, lowerBefore, upperBefore,
4729                                   node->statusArray(), node->depth());
4730        sub[iProb].objectiveValue_ = node->objectiveValue();
4731        sub[iProb].sumInfeasibilities_ = node->sumInfeasibilities();
4732        sub[iProb].numberInfeasibilities_ = node->numberInfeasibilities();
4733#ifdef CHECK_PATH
4734        if (simplex->numberColumns() == numberColumns_Z) {
4735            bool onOptimal = true;
4736            const double * columnLower = simplex->columnLower();
4737            const double * columnUpper = simplex->columnUpper();
4738            for (int i = 0; i < numberColumns_Z; i++) {
4739                if (iNode == gotGoodNode_Z)
4740                    printf("good %d %d %g %g\n", iNode, i, columnLower[i], columnUpper[i]);
4741                if (columnUpper[i] < debuggerSolution_Z[i] || columnLower[i] > debuggerSolution_Z[i] && simplex->isInteger(i)) {
4742                    onOptimal = false;
4743                    break;
4744                }
4745            }
4746            if (onOptimal) {
4747                printf("adding to node %x as %d - objs\n", this, iProb);
4748                for (int j = 0; j <= iProb; j++)
4749                    printf("%d %g\n", j, sub[j].objectiveValue_);
4750            }
4751        }
4752#endif
4753    }
4754    delete [] weight;
4755    delete [] whichNode;
4756    const double * lower = solver->getColLower();
4757    const double * upper = solver->getColUpper();
4758    // restore bounds
4759    for ( int j = 0; j < numberColumns; j++) {
4760        if (lowerBefore[j] != lower[j])
4761            solver->setColLower(j, lowerBefore[j]);
4762        if (upperBefore[j] != upper[j])
4763            solver->setColUpper(j, upperBefore[j]);
4764    }
4765    delete [] upperBefore;
4766    delete [] lowerBefore;
4767    return branch;
4768}
4769
4770// Default Constructor
4771CbcGeneralBranchingObject::CbcGeneralBranchingObject()
4772        : CbcBranchingObject(),
4773        subProblems_(NULL),
4774        node_(NULL),
4775        numberSubProblems_(0),
4776        numberSubLeft_(0),
4777        whichNode_(-1),
4778        numberRows_(0)
4779{
4780    //  printf("CbcGeneral %x default constructor\n",this);
4781}
4782
4783// Useful constructor
4784CbcGeneralBranchingObject::CbcGeneralBranchingObject (CbcModel * model)
4785        : CbcBranchingObject(model, -1, -1, 0.5),
4786        subProblems_(NULL),
4787        node_(NULL),
4788        numberSubProblems_(0),
4789        numberSubLeft_(0),
4790        whichNode_(-1),
4791        numberRows_(0)
4792{
4793    //printf("CbcGeneral %x useful constructor\n",this);
4794}
4795
4796// Copy constructor
4797CbcGeneralBranchingObject::CbcGeneralBranchingObject ( const CbcGeneralBranchingObject & rhs)
4798        : CbcBranchingObject(rhs),
4799        subProblems_(NULL),
4800        node_(rhs.node_),
4801        numberSubProblems_(rhs.numberSubProblems_),
4802        numberSubLeft_(rhs.numberSubLeft_),
4803        whichNode_(rhs.whichNode_),
4804        numberRows_(rhs.numberRows_)
4805{
4806    abort();
4807    if (numberSubProblems_) {
4808        subProblems_ = new CbcSubProblem[numberSubProblems_];
4809        for (int i = 0; i < numberSubProblems_; i++)
4810            subProblems_[i] = rhs.subProblems_[i];
4811    }
4812}
4813
4814// Assignment operator
4815CbcGeneralBranchingObject &
4816CbcGeneralBranchingObject::operator=( const CbcGeneralBranchingObject & rhs)
4817{
4818    if (this != &rhs) {
4819        abort();
4820        CbcBranchingObject::operator=(rhs);
4821        delete [] subProblems_;
4822        numberSubProblems_ = rhs.numberSubProblems_;
4823        numberSubLeft_ = rhs.numberSubLeft_;
4824        whichNode_ = rhs.whichNode_;
4825        numberRows_ = rhs.numberRows_;
4826        if (numberSubProblems_) {
4827            subProblems_ = new CbcSubProblem[numberSubProblems_];
4828            for (int i = 0; i < numberSubProblems_; i++)
4829                subProblems_[i] = rhs.subProblems_[i];
4830        } else {
4831            subProblems_ = NULL;
4832        }
4833        node_ = rhs.node_;
4834    }
4835    return *this;
4836}
4837CbcBranchingObject *
4838CbcGeneralBranchingObject::clone() const
4839{
4840    return (new CbcGeneralBranchingObject(*this));
4841}
4842
4843
4844// Destructor
4845CbcGeneralBranchingObject::~CbcGeneralBranchingObject ()
4846{
4847    //printf("CbcGeneral %x destructor\n",this);
4848    delete [] subProblems_;
4849}
4850bool doingDoneBranch = false;
4851double
4852CbcGeneralBranchingObject::branch()
4853{
4854    double cutoff = model_->getCutoff();
4855    //printf("GenB %x whichNode %d numberLeft %d which %d\n",
4856    // this,whichNode_,numberBranchesLeft(),branchIndex());
4857    if (whichNode_ < 0) {
4858        assert (node_);
4859        bool applied = false;
4860        while (numberBranchesLeft()) {
4861            int which = branchIndex();
4862            decrementNumberBranchesLeft();
4863            CbcSubProblem * thisProb = subProblems_ + which;
4864            if (thisProb->objectiveValue_ < cutoff) {
4865                //printf("branch %x (sub %x) which now %d\n",this,
4866                //     subProblems_,which);
4867                OsiSolverInterface * solver = model_->solver();
4868                thisProb->apply(solver);
4869                OsiClpSolverInterface * clpSolver
4870                = dynamic_cast<OsiClpSolverInterface *> (solver);
4871                assert (clpSolver);
4872                // Move status to basis
4873                clpSolver->setWarmStart(NULL);
4874                //ClpSimplex * simplex = clpSolver->getModelPtr();
4875                node_->setObjectiveValue(thisProb->objectiveValue_);
4876                node_->setSumInfeasibilities(thisProb->sumInfeasibilities_);
4877                node_->setNumberUnsatisfied(thisProb->numberInfeasibilities_);
4878                applied = true;
4879                doingDoneBranch = true;
4880                break;
4881            } else if (numberBranchesLeft()) {
4882                node_->nodeInfo()->branchedOn() ;
4883            }
4884        }
4885        if (!applied) {
4886            // no good one
4887            node_->setObjectiveValue(cutoff + 1.0e20);
4888            node_->setSumInfeasibilities(1.0);
4889            node_->setNumberUnsatisfied(1);
4890            assert (whichNode_ < 0);
4891        }
4892    } else {
4893        decrementNumberBranchesLeft();
4894        CbcSubProblem * thisProb = subProblems_ + whichNode_;
4895        assert (thisProb->objectiveValue_ < cutoff);
4896        OsiSolverInterface * solver = model_->solver();
4897        thisProb->apply(solver);
4898        //OsiClpSolverInterface * clpSolver
4899        //= dynamic_cast<OsiClpSolverInterface *> (solver);
4900        //assert (clpSolver);
4901        // Move status to basis
4902        //clpSolver->setWarmStart(NULL);
4903    }
4904    return 0.0;
4905}
4906/* Double checks in case node can change its mind!
4907   Can change objective etc */
4908void
4909CbcGeneralBranchingObject::checkIsCutoff(double cutoff)
4910{
4911    assert (node_);
4912    int first = branchIndex();
4913    int last = first + numberBranchesLeft();
4914    for (int which = first; which < last; which++) {
4915        CbcSubProblem * thisProb = subProblems_ + which;
4916        if (thisProb->objectiveValue_ < cutoff) {
4917            node_->setObjectiveValue(thisProb->objectiveValue_);
4918            node_->setSumInfeasibilities(thisProb->sumInfeasibilities_);
4919            node_->setNumberUnsatisfied(thisProb->numberInfeasibilities_);
4920            break;
4921        }
4922    }
4923}
4924// Print what would happen
4925void
4926CbcGeneralBranchingObject::print()
4927{
4928    //printf("CbcGeneralObject has %d subproblems\n",numberSubProblems_);
4929}
4930// Fill in current objective etc
4931void
4932CbcGeneralBranchingObject::state(double & objectiveValue,
4933                                 double & sumInfeasibilities,
4934                                 int & numberUnsatisfied, int which) const
4935{
4936    assert (which >= 0 && which < numberSubProblems_);
4937    const CbcSubProblem * thisProb = subProblems_ + which;
4938    objectiveValue = thisProb->objectiveValue_;
4939    sumInfeasibilities = thisProb->sumInfeasibilities_;
4940    numberUnsatisfied = thisProb->numberInfeasibilities_;
4941}
4942/** Compare the original object of \c this with the original object of \c
4943    brObj. Assumes that there is an ordering of the original objects.
4944    This method should be invoked only if \c this and brObj are of the same
4945    type.
4946    Return negative/0/positive depending on whether \c this is
4947    smaller/same/larger than the argument.
4948*/
4949int
4950CbcGeneralBranchingObject::compareOriginalObject
4951(const CbcBranchingObject* /*brObj*/) const
4952{
4953    throw("must implement");
4954}
4955
4956/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
4957    same type and must have the same original object, but they may have
4958    different feasible regions.
4959    Return the appropriate CbcRangeCompare value (first argument being the
4960    sub/superset if that's the case). In case of overlap (and if \c
4961    replaceIfOverlap is true) replace the current branching object with one
4962    whose feasible region is the overlap.
4963*/
4964CbcRangeCompare
4965CbcGeneralBranchingObject::compareBranchingObject
4966(const CbcBranchingObject* /*brObj*/, const bool /*replaceIfOverlap*/)
4967{
4968    throw("must implement");
4969}
4970
4971// Default Constructor
4972CbcOneGeneralBranchingObject::CbcOneGeneralBranchingObject()
4973        : CbcBranchingObject(),
4974        object_(NULL),
4975        whichOne_(-1)
4976{
4977    //printf("CbcOneGeneral %x default constructor\n",this);
4978}
4979
4980// Useful constructor
4981CbcOneGeneralBranchingObject::CbcOneGeneralBranchingObject (CbcModel * model,
4982        CbcGeneralBranchingObject * object,
4983        int whichOne)
4984        : CbcBranchingObject(model, -1, -1, 0.5),
4985        object_(object),
4986        whichOne_(whichOne)
4987{
4988    //printf("CbcOneGeneral %x useful constructor object %x %d left\n",this,
4989    //   object_,object_->numberSubLeft_);
4990    numberBranches_ = 1;
4991}
4992
4993// Copy constructor
4994CbcOneGeneralBranchingObject::CbcOneGeneralBranchingObject ( const CbcOneGeneralBranchingObject & rhs)
4995        : CbcBranchingObject(rhs),
4996        object_(rhs.object_),
4997        whichOne_(rhs.whichOne_)
4998{
4999}
5000
5001// Assignment operator
5002CbcOneGeneralBranchingObject &
5003CbcOneGeneralBranchingObject::operator=( const CbcOneGeneralBranchingObject & rhs)
5004{
5005    if (this != &rhs) {
5006        CbcBranchingObject::operator=(rhs);
5007        object_ = rhs.object_;
5008        whichOne_ = rhs.whichOne_;
5009    }
5010    return *this;
5011}
5012CbcBranchingObject *
5013CbcOneGeneralBranchingObject::clone() const
5014{
5015    return (new CbcOneGeneralBranchingObject(*this));
5016}
5017
5018
5019// Destructor
5020CbcOneGeneralBranchingObject::~CbcOneGeneralBranchingObject ()
5021{
5022    //printf("CbcOneGeneral %x destructor object %x %d left\n",this,
5023    // object_,object_->numberSubLeft_);
5024    assert (object_->numberSubLeft_ > 0 &&
5025            object_->numberSubLeft_ < 1000000);
5026    if (!object_->decrementNumberLeft()) {
5027        // printf("CbcGeneral %x yy destructor\n",object_);
5028        delete object_;
5029    }
5030}
5031double
5032CbcOneGeneralBranchingObject::branch()
5033{
5034    assert (numberBranchesLeft());
5035    decrementNumberBranchesLeft();
5036    assert (!numberBranchesLeft());
5037    object_->setWhichNode(whichOne_);
5038    object_->branch();
5039    return 0.0;
5040}
5041/* Double checks in case node can change its mind!
5042   Can change objective etc */
5043void
5044CbcOneGeneralBranchingObject::checkIsCutoff(double /*cutoff*/)
5045{
5046    assert (numberBranchesLeft());
5047}
5048// Print what would happen
5049void
5050CbcOneGeneralBranchingObject::print()
5051{
5052    //printf("CbcOneGeneralObject has 1 subproblem\n");
5053}
5054/** Compare the original object of \c this with the original object of \c
5055    brObj. Assumes that there is an ordering of the original objects.
5056    This method should be invoked only if \c this and brObj are of the same
5057    type.
5058    Return negative/0/positive depending on whether \c this is
5059    smaller/same/larger than the argument.
5060*/
5061int
5062CbcOneGeneralBranchingObject::compareOriginalObject
5063(const CbcBranchingObject* /*brObj*/) const
5064{
5065    throw("must implement");
5066}
5067
5068/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
5069    same type and must have the same original object, but they may have
5070    different feasible regions.
5071    Return the appropriate CbcRangeCompare value (first argument being the
5072    sub/superset if that's the case). In case of overlap (and if \c
5073    replaceIfOverlap is true) replace the current branching object with one
5074    whose feasible region is the overlap.
5075*/
5076CbcRangeCompare
5077CbcOneGeneralBranchingObject::compareBranchingObject
5078(const CbcBranchingObject* /*brObj*/, const bool /*replaceIfOverlap*/)
5079{
5080    throw("must implement");
5081}
5082// Default Constructor
5083CbcSubProblem::CbcSubProblem()
5084        : objectiveValue_(0.0),
5085        sumInfeasibilities_(0.0),
5086        variables_(NULL),
5087        newBounds_(NULL),
5088        status_(NULL),
5089        depth_(0),
5090        numberChangedBounds_(0),
5091        numberInfeasibilities_(0)
5092{
5093}
5094
5095// Useful constructor
5096CbcSubProblem::CbcSubProblem (const OsiSolverInterface * solver,
5097                              const double * lastLower,
5098                              const double * lastUpper,
5099                              const unsigned char * status,
5100                              int depth)
5101        : objectiveValue_(0.0),
5102        sumInfeasibilities_(0.0),
5103        variables_(NULL),
5104        newBounds_(NULL),
5105        status_(NULL),
5106        depth_(depth),
5107        numberChangedBounds_(0),
5108        numberInfeasibilities_(0)
5109{
5110    const double * lower = solver->getColLower();
5111    const double * upper = solver->getColUpper();
5112
5113    numberChangedBounds_ = 0;
5114    int numberColumns = solver->getNumCols();
5115    int i;
5116    for (i = 0; i < numberColumns; i++) {
5117        if (lower[i] != lastLower[i])
5118            numberChangedBounds_++;
5119        if (upper[i] != lastUpper[i])
5120            numberChangedBounds_++;
5121    }
5122    if (numberChangedBounds_) {
5123        newBounds_ = new double [numberChangedBounds_] ;
5124        variables_ = new int [numberChangedBounds_] ;
5125        numberChangedBounds_ = 0;
5126        for (i = 0; i < numberColumns; i++) {
5127            if (lower[i] != lastLower[i]) {
5128                variables_[numberChangedBounds_] = i;
5129                newBounds_[numberChangedBounds_++] = lower[i];
5130            }
5131            if (upper[i] != lastUpper[i]) {
5132                variables_[numberChangedBounds_] = i | 0x80000000;
5133                newBounds_[numberChangedBounds_++] = upper[i];
5134            }
5135#ifdef CBC_DEBUG
5136            if (lower[i] != lastLower[i]) {
5137                std::cout
5138                    << "lower on " << i << " changed from "
5139                    << lastLower[i] << " to " << lower[i] << std::endl ;
5140            }
5141            if (upper[i] != lastUpper[i]) {
5142                std::cout
5143                    << "upper on " << i << " changed from "
5144                    << lastUpper[i] << " to " << upper[i] << std::endl ;
5145            }
5146#endif
5147        }
5148#ifdef CBC_DEBUG
5149        std::cout << numberChangedBounds_ << " changed bounds." << std::endl ;
5150#endif
5151    }
5152    const OsiClpSolverInterface * clpSolver
5153    = dynamic_cast<const OsiClpSolverInterface *> (solver);
5154    assert (clpSolver);
5155    // Do difference
5156    // Current basis
5157    status_ = clpSolver->getBasis(status);
5158    assert (status_->fullBasis());
5159    //status_->print();
5160}
5161
5162// Copy constructor
5163CbcSubProblem::CbcSubProblem ( const CbcSubProblem & rhs)
5164        : objectiveValue_(rhs.objectiveValue_),
5165        sumInfeasibilities_(rhs.sumInfeasibilities_),
5166        variables_(NULL),
5167        newBounds_(NULL),
5168        status_(NULL),
5169        depth_(rhs.depth_),
5170        numberChangedBounds_(rhs.numberChangedBounds_),
5171        numberInfeasibilities_(rhs.numberInfeasibilities_)
5172{
5173    if (numberChangedBounds_) {
5174        variables_ = CoinCopyOfArray(rhs.variables_, numberChangedBounds_);
5175        newBounds_ = CoinCopyOfArray(rhs.newBounds_, numberChangedBounds_);
5176    }
5177    if (rhs.status_) {
5178        status_ = new CoinWarmStartBasis(*rhs.status_);
5179    }
5180}
5181
5182// Assignment operator
5183CbcSubProblem &
5184CbcSubProblem::operator=( const CbcSubProblem & rhs)
5185{
5186    if (this != &rhs) {
5187        delete [] variables_;
5188        delete [] newBounds_;
5189        delete status_;
5190        objectiveValue_ = rhs.objectiveValue_;
5191        sumInfeasibilities_ = rhs.sumInfeasibilities_;
5192        depth_ = rhs.depth_;
5193        numberChangedBounds_ = rhs.numberChangedBounds_;
5194        numberInfeasibilities_ = rhs.numberInfeasibilities_;
5195        if (numberChangedBounds_) {
5196            variables_ = CoinCopyOfArray(rhs.variables_, numberChangedBounds_);
5197            newBounds_ = CoinCopyOfArray(rhs.newBounds_, numberChangedBounds_);
5198        } else {
5199            variables_ = NULL;
5200            newBounds_ = NULL;
5201        }
5202        if (rhs.status_) {
5203            status_ = new CoinWarmStartBasis(*rhs.status_);
5204        } else {
5205            status_ = NULL;
5206        }
5207    }
5208    return *this;
5209}
5210
5211// Destructor
5212CbcSubProblem::~CbcSubProblem ()
5213{
5214    delete [] variables_;
5215    delete [] newBounds_;
5216    delete status_;
5217}
5218// Apply subproblem
5219void
5220CbcSubProblem::apply(OsiSolverInterface * solver, int what) const
5221{
5222    int i;
5223    if ((what&1) != 0) {
5224        int nSame = 0;
5225        for (i = 0; i < numberChangedBounds_; i++) {
5226            int variable = variables_[i];
5227            int k = variable & 0x3fffffff;
5228            if ((variable&0x80000000) == 0) {
5229                // lower bound changing
5230                //#define CBC_PRINT2
5231#ifdef CBC_PRINT2
5232                if (solver->getColLower()[k] != newBounds_[i])
5233                    printf("lower change for column %d - from %g to %g\n", k, solver->getColLower()[k], newBounds_[i]);
5234#endif
5235#ifndef NDEBUG
5236                if ((variable&0x40000000) == 0 && true) {
5237                    double oldValue = solver->getColLower()[k];
5238                    assert (newBounds_[i] > oldValue - 1.0e-8);
5239                    if (newBounds_[i] < oldValue + 1.0e-8) {
5240#ifdef CBC_PRINT2
5241                        printf("bad null lower change for column %d - bound %g\n", k, oldValue);
5242#endif
5243                        if (newBounds_[i] == oldValue)
5244                            nSame++;
5245                    }
5246                }
5247#endif
5248                solver->setColLower(k, newBounds_[i]);
5249            } else {
5250                // upper bound changing
5251#ifdef CBC_PRINT2
5252                if (solver->getColUpper()[k] != newBounds_[i])
5253                    printf("upper change for column %d - from %g to %g\n", k, solver->getColUpper()[k], newBounds_[i]);
5254#endif
5255#ifndef NDEBUG
5256                if ((variable&0x40000000) == 0 && true) {
5257                    double oldValue = solver->getColUpper()[k];
5258                    assert (newBounds_[i] < oldValue + 1.0e-8);
5259                    if (newBounds_[i] > oldValue - 1.0e-8) {
5260#ifdef CBC_PRINT2
5261                        printf("bad null upper change for column %d - bound %g\n", k, oldValue);
5262#endif
5263                        if (newBounds_[i] == oldValue)
5264                            nSame++;
5265                    }
5266                }
5267#endif
5268                solver->setColUpper(k, newBounds_[i]);
5269            }
5270        }
5271        if (nSame && (nSame < numberChangedBounds_ || (what&3) != 3))
5272            printf("%d changes out of %d redundant %d\n",
5273                   nSame, numberChangedBounds_, what);
5274        else if (numberChangedBounds_ && what == 7 && !nSame)
5275            printf("%d good changes %d\n",
5276                   numberChangedBounds_, what);
5277    }
5278#if 0
5279    if ((what&2) != 0) {
5280        OsiClpSolverInterface * clpSolver
5281        = dynamic_cast<OsiClpSolverInterface *> (solver);
5282        assert (clpSolver);
5283        //assert (clpSolver->getNumRows()==numberRows_);
5284        //clpSolver->setBasis(*status_);
5285        // Current basis
5286        CoinWarmStartBasis * basis = clpSolver->getPointerToWarmStart();
5287        printf("BBBB\n");
5288        basis->print();
5289        assert (basis->fullBasis());
5290        basis->applyDiff(status_);
5291        printf("diff applied %x\n", status_);
5292        printf("CCCC\n");
5293        basis->print();
5294        assert (basis->fullBasis());
5295#ifndef NDEBUG
5296        if (!basis->fullBasis())
5297            printf("Debug this basis!!\n");
5298#endif
5299        clpSolver->setBasis(*basis);
5300    }
5301#endif
5302    if ((what&8) != 0) {
5303        OsiClpSolverInterface * clpSolver
5304        = dynamic_cast<OsiClpSolverInterface *> (solver);
5305        assert (clpSolver);
5306        clpSolver->setBasis(*status_);
5307        delete status_;
5308        status_ = NULL;
5309    }
5310}
5311#endif
5312
5313/** Compare the original object of \c this with the original object of \c
5314    brObj. Assumes that there is an ordering of the original objects.
5315    This method should be invoked only if \c this and brObj are of the same
5316    type.
5317    Return negative/0/positive depending on whether \c this is
5318    smaller/same/larger than the argument.
5319*/
5320int
5321CbcDummyBranchingObject::compareOriginalObject
5322(const CbcBranchingObject* /*brObj*/) const
5323{
5324    throw("must implement");
5325}
5326
5327/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
5328    same type and must have the same original object, but they may have
5329    different feasible regions.
5330    Return the appropriate CbcRangeCompare value (first argument being the
5331    sub/superset if that's the case). In case of overlap (and if \c
5332    replaceIfOverlap is true) replace the current branching object with one
5333    whose feasible region is the overlap.
5334*/
5335CbcRangeCompare
5336CbcDummyBranchingObject::compareBranchingObject
5337(const CbcBranchingObject* /*brObj*/, const bool /*replaceIfOverlap*/)
5338{
5339    throw("must implement");
5340}
Note: See TracBrowser for help on using the repository browser.