source: stable/2.7/Cbc/src/CbcSOS.cpp @ 1675

Last change on this file since 1675 was 1675, checked in by stefan, 8 years ago

sync with trunk rev1674

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