source: branches/sandbox/Cbc/src/CbcSOS.cpp @ 1350

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

Combined CbcSOS with CbcSOSBranchingObject
Combined CbcNWay with CbcNWayBranchingObject
Combined CbcFollowOn? with CbcFixingBranchingObject?

File size: 34.5 KB
Line 
1// Edwin 11/9/2009-- carved out of CbcBranchActual
2
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cstdlib>
9#include <cmath>
10#include <cfloat>
11//#define CBC_DEBUG
12
13#include "CoinTypes.hpp"
14#include "OsiSolverInterface.hpp"
15#include "OsiSolverBranch.hpp"
16#include "CbcModel.hpp"
17#include "CbcMessage.hpp"
18#include "CbcSOS.hpp"
19#include "CbcBranchActual.hpp"
20#include "CoinSort.hpp"
21#include "CoinError.hpp"
22
23//##############################################################################
24
25// Default Constructor
26CbcSOS::CbcSOS ()
27        : CbcObject(),
28        members_(NULL),
29        weights_(NULL),
30        shadowEstimateDown_(1.0),
31        shadowEstimateUp_(1.0),
32        downDynamicPseudoRatio_(0.0),
33        upDynamicPseudoRatio_(0.0),
34        numberTimesDown_(0),
35        numberTimesUp_(0),
36        numberMembers_(0),
37        sosType_(-1),
38        integerValued_(false)
39{
40}
41
42// Useful constructor (which are indices)
43CbcSOS::CbcSOS (CbcModel * model,  int numberMembers,
44                const int * which, const double * weights, int identifier, int type)
45        : CbcObject(model),
46        shadowEstimateDown_(1.0),
47        shadowEstimateUp_(1.0),
48        downDynamicPseudoRatio_(0.0),
49        upDynamicPseudoRatio_(0.0),
50        numberTimesDown_(0),
51        numberTimesUp_(0),
52        numberMembers_(numberMembers),
53        sosType_(type)
54{
55    id_ = identifier;
56    integerValued_ = type == 1;
57    if (integerValued_) {
58        // check all members integer
59        OsiSolverInterface * solver = model->solver();
60        if (solver) {
61            for (int i = 0; i < numberMembers_; i++) {
62                if (!solver->isInteger(which[i]))
63                    integerValued_ = false;
64            }
65        } else {
66            // can't tell
67            integerValued_ = false;
68        }
69    }
70    if (numberMembers_) {
71        members_ = new int[numberMembers_];
72        weights_ = new double[numberMembers_];
73        memcpy(members_, which, numberMembers_*sizeof(int));
74        if (weights) {
75            memcpy(weights_, weights, numberMembers_*sizeof(double));
76        } else {
77            for (int i = 0; i < numberMembers_; i++)
78                weights_[i] = i;
79        }
80        // sort so weights increasing
81        CoinSort_2(weights_, weights_ + numberMembers_, members_);
82        double last = -COIN_DBL_MAX;
83        int i;
84        for (i = 0; i < numberMembers_; i++) {
85            double possible = CoinMax(last + 1.0e-10, weights_[i]);
86            weights_[i] = possible;
87            last = possible;
88        }
89    } else {
90        members_ = NULL;
91        weights_ = NULL;
92    }
93    assert (sosType_ > 0 && sosType_ < 3);
94}
95
96// Copy constructor
97CbcSOS::CbcSOS ( const CbcSOS & rhs)
98        : CbcObject(rhs)
99{
100    shadowEstimateDown_ = rhs.shadowEstimateDown_;
101    shadowEstimateUp_ = rhs.shadowEstimateUp_;
102    downDynamicPseudoRatio_ = rhs.downDynamicPseudoRatio_;
103    upDynamicPseudoRatio_ = rhs.upDynamicPseudoRatio_;
104    numberTimesDown_ = rhs.numberTimesDown_;
105    numberTimesUp_ = rhs.numberTimesUp_;
106    numberMembers_ = rhs.numberMembers_;
107    sosType_ = rhs.sosType_;
108    integerValued_ = rhs.integerValued_;
109    if (numberMembers_) {
110        members_ = new int[numberMembers_];
111        weights_ = new double[numberMembers_];
112        memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
113        memcpy(weights_, rhs.weights_, numberMembers_*sizeof(double));
114    } else {
115        members_ = NULL;
116        weights_ = NULL;
117    }
118}
119
120// Clone
121CbcObject *
122CbcSOS::clone() const
123{
124    return new CbcSOS(*this);
125}
126
127// Assignment operator
128CbcSOS &
129CbcSOS::operator=( const CbcSOS & rhs)
130{
131    if (this != &rhs) {
132        CbcObject::operator=(rhs);
133        delete [] members_;
134        delete [] weights_;
135        shadowEstimateDown_ = rhs.shadowEstimateDown_;
136        shadowEstimateUp_ = rhs.shadowEstimateUp_;
137        downDynamicPseudoRatio_ = rhs.downDynamicPseudoRatio_;
138        upDynamicPseudoRatio_ = rhs.upDynamicPseudoRatio_;
139        numberTimesDown_ = rhs.numberTimesDown_;
140        numberTimesUp_ = rhs.numberTimesUp_;
141        numberMembers_ = rhs.numberMembers_;
142        sosType_ = rhs.sosType_;
143        integerValued_ = rhs.integerValued_;
144        if (numberMembers_) {
145            members_ = new int[numberMembers_];
146            weights_ = new double[numberMembers_];
147            memcpy(members_, rhs.members_, numberMembers_*sizeof(int));
148            memcpy(weights_, rhs.weights_, numberMembers_*sizeof(double));
149        } else {
150            members_ = NULL;
151            weights_ = NULL;
152        }
153    }
154    return *this;
155}
156
157// Destructor
158CbcSOS::~CbcSOS ()
159{
160    delete [] members_;
161    delete [] weights_;
162}
163double
164CbcSOS::infeasibility(const OsiBranchingInformation * info,
165                      int &preferredWay) const
166{
167    int j;
168    int firstNonZero = -1;
169    int lastNonZero = -1;
170    OsiSolverInterface * solver = model_->solver();
171    const double * solution = model_->testSolution();
172    //const double * lower = solver->getColLower();
173    const double * upper = solver->getColUpper();
174    //double largestValue=0.0;
175    double integerTolerance =
176        model_->getDblParam(CbcModel::CbcIntegerTolerance);
177    double weight = 0.0;
178    double sum = 0.0;
179
180    // check bounds etc
181    double lastWeight = -1.0e100;
182    for (j = 0; j < numberMembers_; j++) {
183        int iColumn = members_[j];
184        if (lastWeight >= weights_[j] - 1.0e-7)
185            throw CoinError("Weights too close together in SOS", "infeasibility", "CbcSOS");
186        double value = CoinMax(0.0, solution[iColumn]);
187        sum += value;
188        if (value > integerTolerance && upper[iColumn]) {
189            // Possibly due to scaling a fixed variable might slip through
190            if (value > upper[iColumn]) {
191                value = upper[iColumn];
192                // Could change to #ifdef CBC_DEBUG
193#ifndef NDEBUG
194                if (model_->messageHandler()->logLevel() > 2)
195                    printf("** Variable %d (%d) has value %g and upper bound of %g\n",
196                           iColumn, j, value, upper[iColumn]);
197#endif
198            }
199            weight += weights_[j] * value;
200            if (firstNonZero < 0)
201                firstNonZero = j;
202            lastNonZero = j;
203        }
204    }
205    preferredWay = 1;
206    if (lastNonZero - firstNonZero >= sosType_) {
207        // find where to branch
208        assert (sum > 0.0);
209        weight /= sum;
210        if (info->defaultDual_ >= 0.0 && info->usefulRegion_ && info->columnStart_) {
211            assert (sosType_ == 1);
212            int iWhere;
213            for (iWhere = firstNonZero; iWhere < lastNonZero - 1; iWhere++) {
214                if (weight < weights_[iWhere+1]) {
215                    break;
216                }
217            }
218            int jColumnDown = members_[iWhere];
219            int jColumnUp = members_[iWhere+1];
220            int n = 0;
221            CoinBigIndex j;
222            double objMove = info->objective_[jColumnDown];
223            for (j = info->columnStart_[jColumnDown];
224                    j < info->columnStart_[jColumnDown] + info->columnLength_[jColumnDown]; j++) {
225                double value = info->elementByColumn_[j];
226                int iRow = info->row_[j];
227                info->indexRegion_[n++] = iRow;
228                info->usefulRegion_[iRow] = value;
229            }
230            for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++) {
231                int jColumn = members_[iWhere];
232                double solValue = info->solution_[jColumn];
233                if (!solValue)
234                    continue;
235                objMove -= info->objective_[jColumn] * solValue;
236                for (j = info->columnStart_[jColumn];
237                        j < info->columnStart_[jColumn] + info->columnLength_[jColumn]; j++) {
238                    double value = -info->elementByColumn_[j] * solValue;
239                    int iRow = info->row_[j];
240                    double oldValue = info->usefulRegion_[iRow];
241                    if (!oldValue) {
242                        info->indexRegion_[n++] = iRow;
243                    } else {
244                        value += oldValue;
245                        if (!value)
246                            value = 1.0e-100;
247                    }
248                    info->usefulRegion_[iRow] = value;
249                }
250            }
251            const double * pi = info->pi_;
252            const double * activity = info->rowActivity_;
253            const double * lower = info->rowLower_;
254            const double * upper = info->rowUpper_;
255            double tolerance = info->primalTolerance_;
256            double direction = info->direction_;
257            shadowEstimateDown_ = objMove * direction;
258            bool infeasible = false;
259            for (int k = 0; k < n; k++) {
260                int iRow = info->indexRegion_[k];
261                double movement = info->usefulRegion_[iRow];
262                // not this time info->usefulRegion_[iRow]=0.0;
263                if (lower[iRow] < -1.0e20)
264                    assert (pi[iRow] <= 1.0e-3);
265                if (upper[iRow] > 1.0e20)
266                    assert (pi[iRow] >= -1.0e-3);
267                double valueP = pi[iRow] * direction;
268                // if move makes infeasible then make at least default
269                double newValue = activity[iRow] + movement;
270                if (newValue > upper[iRow] + tolerance || newValue < lower[iRow] - tolerance) {
271                    shadowEstimateDown_ += fabs(movement) * CoinMax(fabs(valueP), info->defaultDual_);
272                    infeasible = true;
273                }
274            }
275            if (shadowEstimateDown_ < info->integerTolerance_) {
276                if (!infeasible) {
277                    shadowEstimateDown_ = 1.0e-10;
278#ifdef COIN_DEVELOP
279                    printf("zero pseudoShadowPrice\n");
280#endif
281                } else
282                    shadowEstimateDown_ = info->integerTolerance_;
283            }
284            // And other way
285            // take off
286            objMove -= info->objective_[jColumnDown];
287            for (j = info->columnStart_[jColumnDown];
288                    j < info->columnStart_[jColumnDown] + info->columnLength_[jColumnDown]; j++) {
289                double value = -info->elementByColumn_[j];
290                int iRow = info->row_[j];
291                double oldValue = info->usefulRegion_[iRow];
292                if (!oldValue) {
293                    info->indexRegion_[n++] = iRow;
294                } else {
295                    value += oldValue;
296                    if (!value)
297                        value = 1.0e-100;
298                }
299                info->usefulRegion_[iRow] = value;
300            }
301            // add on
302            objMove += info->objective_[jColumnUp];
303            for (j = info->columnStart_[jColumnUp];
304                    j < info->columnStart_[jColumnUp] + info->columnLength_[jColumnUp]; j++) {
305                double value = info->elementByColumn_[j];
306                int iRow = info->row_[j];
307                double oldValue = info->usefulRegion_[iRow];
308                if (!oldValue) {
309                    info->indexRegion_[n++] = iRow;
310                } else {
311                    value += oldValue;
312                    if (!value)
313                        value = 1.0e-100;
314                }
315                info->usefulRegion_[iRow] = value;
316            }
317            shadowEstimateUp_ = objMove * direction;
318            infeasible = false;
319            for (int k = 0; k < n; k++) {
320                int iRow = info->indexRegion_[k];
321                double movement = info->usefulRegion_[iRow];
322                info->usefulRegion_[iRow] = 0.0;
323                if (lower[iRow] < -1.0e20)
324                    assert (pi[iRow] <= 1.0e-3);
325                if (upper[iRow] > 1.0e20)
326                    assert (pi[iRow] >= -1.0e-3);
327                double valueP = pi[iRow] * direction;
328                // if move makes infeasible then make at least default
329                double newValue = activity[iRow] + movement;
330                if (newValue > upper[iRow] + tolerance || newValue < lower[iRow] - tolerance) {
331                    shadowEstimateUp_ += fabs(movement) * CoinMax(fabs(valueP), info->defaultDual_);
332                    infeasible = true;
333                }
334            }
335            if (shadowEstimateUp_ < info->integerTolerance_) {
336                if (!infeasible) {
337                    shadowEstimateUp_ = 1.0e-10;
338#ifdef COIN_DEVELOP
339                    printf("zero pseudoShadowPrice\n");
340#endif
341                } else
342                    shadowEstimateUp_ = info->integerTolerance_;
343            }
344            // adjust
345            double downCost = shadowEstimateDown_;
346            double upCost = shadowEstimateUp_;
347            if (numberTimesDown_)
348                downCost *= downDynamicPseudoRatio_ /
349                            static_cast<double> (numberTimesDown_);
350            if (numberTimesUp_)
351                upCost *= upDynamicPseudoRatio_ /
352                          static_cast<double> (numberTimesUp_);
353#define WEIGHT_AFTER 0.7
354#define WEIGHT_BEFORE 0.1
355            int stateOfSearch = model_->stateOfSearch() % 10;
356            double returnValue = 0.0;
357            double minValue = CoinMin(downCost, upCost);
358            double maxValue = CoinMax(downCost, upCost);
359            if (stateOfSearch <= 2) {
360                // no branching solution
361                returnValue = WEIGHT_BEFORE * minValue + (1.0 - WEIGHT_BEFORE) * maxValue;
362            } else {
363                returnValue = WEIGHT_AFTER * minValue + (1.0 - WEIGHT_AFTER) * maxValue;
364            }
365#ifdef PRINT_SHADOW
366            printf("%d id - down %d %g up %d %g shadow %g, %g returned %g\n",
367                   id_, numberTimesDown_, downDynamicPseudoRatio_,
368                   numberTimesUp_, upDynamicPseudoRatio_, shadowEstimateDown_,
369                   shadowEstimateUp_, returnValue);
370#endif
371            return returnValue;
372        } else {
373            double value = lastNonZero - firstNonZero + 1;
374            value *= 0.5 / static_cast<double> (numberMembers_);
375            return value;
376        }
377    } else {
378        return 0.0; // satisfied
379    }
380}
381
382// This looks at solution and sets bounds to contain solution
383void
384CbcSOS::feasibleRegion()
385{
386    int j;
387    int firstNonZero = -1;
388    int lastNonZero = -1;
389    OsiSolverInterface * solver = model_->solver();
390    const double * solution = model_->testSolution();
391    //const double * lower = solver->getColLower();
392    const double * upper = solver->getColUpper();
393    double integerTolerance =
394        model_->getDblParam(CbcModel::CbcIntegerTolerance);
395    double weight = 0.0;
396    double sum = 0.0;
397
398    for (j = 0; j < numberMembers_; j++) {
399        int iColumn = members_[j];
400        double value = CoinMax(0.0, solution[iColumn]);
401        sum += value;
402        if (value > integerTolerance && upper[iColumn]) {
403            weight += weights_[j] * value;
404            if (firstNonZero < 0)
405                firstNonZero = j;
406            lastNonZero = j;
407        }
408    }
409    assert (lastNonZero - firstNonZero < sosType_) ;
410    for (j = 0; j < firstNonZero; j++) {
411        int iColumn = members_[j];
412        solver->setColUpper(iColumn, 0.0);
413    }
414    for (j = lastNonZero + 1; j < numberMembers_; j++) {
415        int iColumn = members_[j];
416        solver->setColUpper(iColumn, 0.0);
417    }
418}
419// Redoes data when sequence numbers change
420void
421CbcSOS::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
422{
423    model_ = model;
424    int n2 = 0;
425    for (int j = 0; j < numberMembers_; j++) {
426        int iColumn = members_[j];
427        int i;
428        for (i = 0; i < numberColumns; i++) {
429            if (originalColumns[i] == iColumn)
430                break;
431        }
432        if (i < numberColumns) {
433            members_[n2] = i;
434            weights_[n2++] = weights_[j];
435        }
436    }
437    if (n2 < numberMembers_) {
438        //printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
439        numberMembers_ = n2;
440    }
441}
442CbcBranchingObject *
443CbcSOS::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way)
444{
445    int j;
446    const double * solution = model_->testSolution();
447    double integerTolerance =
448        model_->getDblParam(CbcModel::CbcIntegerTolerance);
449    //OsiSolverInterface * solver = model_->solver();
450    const double * upper = solver->getColUpper();
451    int firstNonFixed = -1;
452    int lastNonFixed = -1;
453    int firstNonZero = -1;
454    int lastNonZero = -1;
455    double weight = 0.0;
456    double sum = 0.0;
457    for (j = 0; j < numberMembers_; j++) {
458        int iColumn = members_[j];
459        if (upper[iColumn]) {
460            double value = CoinMax(0.0, solution[iColumn]);
461            sum += value;
462            if (firstNonFixed < 0)
463                firstNonFixed = j;
464            lastNonFixed = j;
465            if (value > integerTolerance) {
466                weight += weights_[j] * value;
467                if (firstNonZero < 0)
468                    firstNonZero = j;
469                lastNonZero = j;
470            }
471        }
472    }
473    assert (lastNonZero - firstNonZero >= sosType_) ;
474    // find where to branch
475    assert (sum > 0.0);
476    weight /= sum;
477    int iWhere;
478    double separator = 0.0;
479    for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
480        if (weight < weights_[iWhere+1])
481            break;
482    if (sosType_ == 1) {
483        // SOS 1
484        separator = 0.5 * (weights_[iWhere] + weights_[iWhere+1]);
485    } else {
486        // SOS 2
487        if (iWhere == firstNonFixed)
488            iWhere++;;
489        if (iWhere == lastNonFixed - 1)
490            iWhere = lastNonFixed - 2;
491        separator = weights_[iWhere+1];
492    }
493    // create object
494    CbcBranchingObject * branch;
495    branch = new CbcSOSBranchingObject(model_, this, way, separator);
496    branch->setOriginalObject(this);
497    return branch;
498}
499/* Pass in information on branch just done and create CbcObjectUpdateData instance.
500   If object does not need data then backward pointer will be NULL.
501   Assumes can get information from solver */
502CbcObjectUpdateData
503CbcSOS::createUpdateInformation(const OsiSolverInterface * solver,
504                                const CbcNode * node,
505                                const CbcBranchingObject * branchingObject)
506{
507    double originalValue = node->objectiveValue();
508    int originalUnsatisfied = node->numberUnsatisfied();
509    double objectiveValue = solver->getObjValue() * solver->getObjSense();
510    int unsatisfied = 0;
511    int i;
512    //might be base model - doesn't matter
513    int numberIntegers = model_->numberIntegers();;
514    const double * solution = solver->getColSolution();
515    //const double * lower = solver->getColLower();
516    //const double * upper = solver->getColUpper();
517    double change = CoinMax(0.0, objectiveValue - originalValue);
518    int iStatus;
519    if (solver->isProvenOptimal())
520        iStatus = 0; // optimal
521    else if (solver->isIterationLimitReached()
522             && !solver->isDualObjectiveLimitReached())
523        iStatus = 2; // unknown
524    else
525        iStatus = 1; // infeasible
526
527    bool feasible = iStatus != 1;
528    if (feasible) {
529        double integerTolerance =
530            model_->getDblParam(CbcModel::CbcIntegerTolerance);
531        const int * integerVariable = model_->integerVariable();
532        for (i = 0; i < numberIntegers; i++) {
533            int j = integerVariable[i];
534            double value = solution[j];
535            double nearest = floor(value + 0.5);
536            if (fabs(value - nearest) > integerTolerance)
537                unsatisfied++;
538        }
539    }
540    int way = branchingObject->way();
541    way = - way; // because after branch so moved on
542    double value = branchingObject->value();
543    CbcObjectUpdateData newData (this, way,
544                                 change, iStatus,
545                                 originalUnsatisfied - unsatisfied, value);
546    newData.originalObjective_ = originalValue;
547    // Solvers know about direction
548    double direction = solver->getObjSense();
549    solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
550    newData.cutoff_ *= direction;
551    return newData;
552}
553// Update object by CbcObjectUpdateData
554void
555CbcSOS::updateInformation(const CbcObjectUpdateData & data)
556{
557    bool feasible = data.status_ != 1;
558    int way = data.way_;
559    //double value = data.branchingValue_;
560    double originalValue = data.originalObjective_;
561    double change = data.change_;
562    if (way < 0) {
563        // down
564        if (!feasible) {
565            double distanceToCutoff = 0.0;
566            //double objectiveValue = model_->getCurrentMinimizationObjValue();
567            distanceToCutoff =  model_->getCutoff()  - originalValue;
568            if (distanceToCutoff < 1.0e20)
569                change = distanceToCutoff * 2.0;
570            else
571                change = (downDynamicPseudoRatio_ * shadowEstimateDown_ + 1.0e-3) * 10.0;
572        }
573        change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
574#ifdef PRINT_SHADOW
575        if (numberTimesDown_)
576            printf("Updating id %d - down change %g (true %g) - ndown %d estimated change %g - raw shadow estimate %g\n",
577                   id_, change, data.change_, numberTimesDown_, shadowEstimateDown_*
578                   (downDynamicPseudoRatio_ / ((double) numberTimesDown_)),
579                   shadowEstimateDown_);
580        else
581            printf("Updating id %d - down change %g (true %g) - shadow estimate %g\n",
582                   id_, change, data.change_, shadowEstimateDown_);
583#endif
584        numberTimesDown_++;
585        downDynamicPseudoRatio_ += change / shadowEstimateDown_;
586    } else {
587        // up
588        if (!feasible) {
589            double distanceToCutoff = 0.0;
590            //double objectiveValue = model_->getCurrentMinimizationObjValue();
591            distanceToCutoff =  model_->getCutoff()  - originalValue;
592            if (distanceToCutoff < 1.0e20)
593                change = distanceToCutoff * 2.0;
594            else
595                change = (upDynamicPseudoRatio_ * shadowEstimateUp_ + 1.0e-3) * 10.0;
596        }
597        change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
598#ifdef PRINT_SHADOW
599        if (numberTimesUp_)
600            printf("Updating id %d - up change %g (true %g) - nup %d estimated change %g - raw shadow estimate %g\n",
601                   id_, change, data.change_, numberTimesUp_, shadowEstimateUp_*
602                   (upDynamicPseudoRatio_ / ((double) numberTimesUp_)),
603                   shadowEstimateUp_);
604        else
605            printf("Updating id %d - up change %g (true %g) - shadow estimate %g\n",
606                   id_, change, data.change_, shadowEstimateUp_);
607#endif
608        numberTimesUp_++;
609        upDynamicPseudoRatio_ += change / shadowEstimateUp_;
610    }
611}
612
613/* Create an OsiSolverBranch object
614
615This returns NULL if branch not represented by bound changes
616*/
617OsiSolverBranch *
618CbcSOS::solverBranch() const
619{
620    int j;
621    const double * solution = model_->testSolution();
622    double integerTolerance =
623        model_->getDblParam(CbcModel::CbcIntegerTolerance);
624    OsiSolverInterface * solver = model_->solver();
625    const double * upper = solver->getColUpper();
626    int firstNonFixed = -1;
627    int lastNonFixed = -1;
628    int firstNonZero = -1;
629    int lastNonZero = -1;
630    double weight = 0.0;
631    double sum = 0.0;
632    double * fix = new double[numberMembers_];
633    int * which = new int[numberMembers_];
634    for (j = 0; j < numberMembers_; j++) {
635        int iColumn = members_[j];
636        // fix all on one side or other (even if fixed)
637        fix[j] = 0.0;
638        which[j] = iColumn;
639        if (upper[iColumn]) {
640            double value = CoinMax(0.0, solution[iColumn]);
641            sum += value;
642            if (firstNonFixed < 0)
643                firstNonFixed = j;
644            lastNonFixed = j;
645            if (value > integerTolerance) {
646                weight += weights_[j] * value;
647                if (firstNonZero < 0)
648                    firstNonZero = j;
649                lastNonZero = j;
650            }
651        }
652    }
653    assert (lastNonZero - firstNonZero >= sosType_) ;
654    // find where to branch
655    assert (sum > 0.0);
656    weight /= sum;
657    // down branch fixes ones above weight to 0
658    int iWhere;
659    int iDownStart = 0;
660    int iUpEnd = 0;
661    for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
662        if (weight < weights_[iWhere+1])
663            break;
664    if (sosType_ == 1) {
665        // SOS 1
666        iUpEnd = iWhere + 1;
667        iDownStart = iUpEnd;
668    } else {
669        // SOS 2
670        if (iWhere == firstNonFixed)
671            iWhere++;;
672        if (iWhere == lastNonFixed - 1)
673            iWhere = lastNonFixed - 2;
674        iUpEnd = iWhere + 1;
675        iDownStart = iUpEnd + 1;
676    }
677    //
678    OsiSolverBranch * branch = new OsiSolverBranch();
679    branch->addBranch(-1, 0, NULL, NULL, numberMembers_ - iDownStart, which + iDownStart, fix);
680    branch->addBranch(1, 0, NULL, NULL, iUpEnd, which, fix);
681    delete [] fix;
682    delete [] which;
683    return branch;
684}
685// Construct an OsiSOS object
686OsiSOS *
687CbcSOS::osiObject(const OsiSolverInterface * solver) const
688{
689    OsiSOS * obj = new OsiSOS(solver, numberMembers_, members_, weights_, sosType_);
690    obj->setPriority(priority());
691    return obj;
692}
693
694// Default Constructor
695CbcSOSBranchingObject::CbcSOSBranchingObject()
696        : CbcBranchingObject(),
697        firstNonzero_(-1),
698        lastNonzero_(-1)
699{
700    set_ = NULL;
701    separator_ = 0.0;
702}
703
704// Useful constructor
705CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
706        const CbcSOS * set,
707        int way ,
708        double separator)
709        : CbcBranchingObject(model, set->id(), way, 0.5)
710{
711    set_ = set;
712    separator_ = separator;
713    computeNonzeroRange();
714}
715
716// Copy constructor
717CbcSOSBranchingObject::CbcSOSBranchingObject (const CbcSOSBranchingObject & rhs)
718        : CbcBranchingObject(rhs),
719        firstNonzero_(rhs.firstNonzero_),
720        lastNonzero_(rhs.lastNonzero_)
721{
722    set_ = rhs.set_;
723    separator_ = rhs.separator_;
724}
725
726// Assignment operator
727CbcSOSBranchingObject &
728CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject & rhs)
729{
730    if (this != &rhs) {
731        CbcBranchingObject::operator=(rhs);
732        set_ = rhs.set_;
733        separator_ = rhs.separator_;
734        firstNonzero_ = rhs.firstNonzero_;
735        lastNonzero_ = rhs.lastNonzero_;
736    }
737    return *this;
738}
739CbcBranchingObject *
740CbcSOSBranchingObject::clone() const
741{
742    return (new CbcSOSBranchingObject(*this));
743}
744
745
746// Destructor
747CbcSOSBranchingObject::~CbcSOSBranchingObject ()
748{
749}
750
751void
752CbcSOSBranchingObject::computeNonzeroRange()
753{
754    const int numberMembers = set_->numberMembers();
755    const double * weights = set_->weights();
756    int i = 0;
757    if (way_ < 0) {
758        for ( i = 0; i < numberMembers; i++) {
759            if (weights[i] > separator_)
760                break;
761        }
762        assert (i < numberMembers);
763        firstNonzero_ = 0;
764        lastNonzero_ = i;
765    } else {
766        for ( i = 0; i < numberMembers; i++) {
767            if (weights[i] >= separator_)
768                break;
769        }
770        assert (i < numberMembers);
771        firstNonzero_ = i;
772        lastNonzero_ = numberMembers;
773    }
774}
775
776double
777CbcSOSBranchingObject::branch()
778{
779    decrementNumberBranchesLeft();
780    int numberMembers = set_->numberMembers();
781    const int * which = set_->members();
782    const double * weights = set_->weights();
783    OsiSolverInterface * solver = model_->solver();
784    //const double * lower = solver->getColLower();
785    //const double * upper = solver->getColUpper();
786    // *** for way - up means fix all those in down section
787    if (way_ < 0) {
788        int i;
789        for ( i = 0; i < numberMembers; i++) {
790            if (weights[i] > separator_)
791                break;
792        }
793        assert (i < numberMembers);
794        for (; i < numberMembers; i++)
795            solver->setColUpper(which[i], 0.0);
796        way_ = 1;         // Swap direction
797    } else {
798        int i;
799        for ( i = 0; i < numberMembers; i++) {
800            if (weights[i] >= separator_)
801                break;
802            else
803                solver->setColUpper(which[i], 0.0);
804        }
805        assert (i < numberMembers);
806        way_ = -1;        // Swap direction
807    }
808    computeNonzeroRange();
809    return 0.0;
810}
811/* Update bounds in solver as in 'branch' and update given bounds.
812   branchState is -1 for 'down' +1 for 'up' */
813void
814CbcSOSBranchingObject::fix(OsiSolverInterface * solver,
815                           double * /*lower*/, double * upper,
816                           int branchState) const
817{
818    int numberMembers = set_->numberMembers();
819    const int * which = set_->members();
820    const double * weights = set_->weights();
821    //const double * lower = solver->getColLower();
822    //const double * upper = solver->getColUpper();
823    // *** for way - up means fix all those in down section
824    if (branchState < 0) {
825        int i;
826        for ( i = 0; i < numberMembers; i++) {
827            if (weights[i] > separator_)
828                break;
829        }
830        assert (i < numberMembers);
831        for (; i < numberMembers; i++) {
832            solver->setColUpper(which[i], 0.0);
833            upper[which[i]] = 0.0;
834        }
835    } else {
836        int i;
837        for ( i = 0; i < numberMembers; i++) {
838            if (weights[i] >= separator_) {
839                break;
840            } else {
841                solver->setColUpper(which[i], 0.0);
842                upper[which[i]] = 0.0;
843            }
844        }
845        assert (i < numberMembers);
846    }
847}
848// Print what would happen
849void
850CbcSOSBranchingObject::print()
851{
852    int numberMembers = set_->numberMembers();
853    const int * which = set_->members();
854    const double * weights = set_->weights();
855    OsiSolverInterface * solver = model_->solver();
856    //const double * lower = solver->getColLower();
857    const double * upper = solver->getColUpper();
858    int first = numberMembers;
859    int last = -1;
860    int numberFixed = 0;
861    int numberOther = 0;
862    int i;
863    for ( i = 0; i < numberMembers; i++) {
864        double bound = upper[which[i]];
865        if (bound) {
866            first = CoinMin(first, i);
867            last = CoinMax(last, i);
868        }
869    }
870    // *** for way - up means fix all those in down section
871    if (way_ < 0) {
872        printf("SOS Down");
873        for ( i = 0; i < numberMembers; i++) {
874            double bound = upper[which[i]];
875            if (weights[i] > separator_)
876                break;
877            else if (bound)
878                numberOther++;
879        }
880        assert (i < numberMembers);
881        for (; i < numberMembers; i++) {
882            double bound = upper[which[i]];
883            if (bound)
884                numberFixed++;
885        }
886    } else {
887        printf("SOS Up");
888        for ( i = 0; i < numberMembers; i++) {
889            double bound = upper[which[i]];
890            if (weights[i] >= separator_)
891                break;
892            else if (bound)
893                numberFixed++;
894        }
895        assert (i < numberMembers);
896        for (; i < numberMembers; i++) {
897            double bound = upper[which[i]];
898            if (bound)
899                numberOther++;
900        }
901    }
902    printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
903           separator_, which[first], weights[first], which[last], weights[last], numberFixed, numberOther);
904}
905
906/** Compare the original object of \c this with the original object of \c
907    brObj. Assumes that there is an ordering of the original objects.
908    This method should be invoked only if \c this and brObj are of the same
909    type.
910    Return negative/0/positive depending on whether \c this is
911    smaller/same/larger than the argument.
912*/
913int
914CbcSOSBranchingObject::compareOriginalObject
915(const CbcBranchingObject* brObj) const
916{
917    const CbcSOSBranchingObject* br =
918        dynamic_cast<const CbcSOSBranchingObject*>(brObj);
919    assert(br);
920    const CbcSOS* s0 = set_;
921    const CbcSOS* s1 = br->set_;
922    if (s0->sosType() != s1->sosType()) {
923        return s0->sosType() - s1->sosType();
924    }
925    if (s0->numberMembers() != s1->numberMembers()) {
926        return s0->numberMembers() - s1->numberMembers();
927    }
928    const int memberCmp = memcmp(s0->members(), s1->members(),
929                                 s0->numberMembers() * sizeof(int));
930    if (memberCmp != 0) {
931        return memberCmp;
932    }
933    return memcmp(s0->weights(), s1->weights(),
934                  s0->numberMembers() * sizeof(double));
935}
936
937/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
938    same type and must have the same original object, but they may have
939    different feasible regions.
940    Return the appropriate CbcRangeCompare value (first argument being the
941    sub/superset if that's the case). In case of overlap (and if \c
942    replaceIfOverlap is true) replace the current branching object with one
943    whose feasible region is the overlap.
944*/
945CbcRangeCompare
946CbcSOSBranchingObject::compareBranchingObject
947(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
948{
949    const CbcSOSBranchingObject* br =
950        dynamic_cast<const CbcSOSBranchingObject*>(brObj);
951    assert(br);
952    if (firstNonzero_ < br->firstNonzero_) {
953        if (lastNonzero_ >= br->lastNonzero_) {
954            return CbcRangeSuperset;
955        } else if (lastNonzero_ <= br->firstNonzero_) {
956            return CbcRangeDisjoint;
957        } else {
958            // overlap
959            if (replaceIfOverlap) {
960                firstNonzero_ = br->firstNonzero_;
961            }
962            return CbcRangeOverlap;
963        }
964    } else if (firstNonzero_ > br->firstNonzero_) {
965        if (lastNonzero_ <= br->lastNonzero_) {
966            return CbcRangeSubset;
967        } else if (firstNonzero_ >= br->lastNonzero_) {
968            return CbcRangeDisjoint;
969        } else {
970            // overlap
971            if (replaceIfOverlap) {
972                lastNonzero_ = br->lastNonzero_;
973            }
974            return CbcRangeOverlap;
975        }
976    } else {
977        if (lastNonzero_ == br->lastNonzero_) {
978            return CbcRangeSame;
979        }
980        return lastNonzero_ < br->lastNonzero_ ? CbcRangeSubset : CbcRangeSuperset;
981    }
982    return CbcRangeSame; // fake return
983}
Note: See TracBrowser for help on using the repository browser.