source: trunk/Cbc/src/CbcSOS.cpp @ 1433

Last change on this file since 1433 was 1432, checked in by bjarni, 10 years ago

Added extra return at end of each source file where needed, to remove possible linefeed conflicts (NightlyBuild? errors)

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