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

Last change on this file since 2367 was 2367, checked in by forrest, 3 years ago

take out assert

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