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

Last change on this file since 2364 was 2364, checked in by forrest, 16 months ago

should not be tolerance in SOS

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.6 KB
Line 
1// $Id: CbcSOS.cpp 2364 2018-02-22 12:01:40Z 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    assert (iWhere<lastNonZero);
589    if (sosType_ == 1) {
590        // SOS 1
591        separator = 0.5 * (weights_[iWhere] + weights_[iWhere+1]);
592    } else {
593        // SOS 2
594        if (iWhere == firstNonFixed)
595            iWhere++;;
596        if (iWhere == lastNonFixed - 1)
597            iWhere = lastNonFixed - 2;
598        separator = weights_[iWhere+1];
599    }
600    // create object
601    CbcBranchingObject * branch;
602    branch = new CbcSOSBranchingObject(model_, this, way, separator);
603    branch->setOriginalObject(this);
604    return branch;
605}
606/* Pass in information on branch just done and create CbcObjectUpdateData instance.
607   If object does not need data then backward pointer will be NULL.
608   Assumes can get information from solver */
609CbcObjectUpdateData
610CbcSOS::createUpdateInformation(const OsiSolverInterface * solver,
611                                const CbcNode * node,
612                                const CbcBranchingObject * branchingObject)
613{
614    double originalValue = node->objectiveValue();
615    int originalUnsatisfied = node->numberUnsatisfied();
616    double objectiveValue = solver->getObjValue() * solver->getObjSense();
617    int unsatisfied = 0;
618    int i;
619    //might be base model - doesn't matter
620    int numberIntegers = model_->numberIntegers();;
621    const double * solution = solver->getColSolution();
622    //const double * lower = solver->getColLower();
623    //const double * upper = solver->getColUpper();
624    double change = CoinMax(0.0, objectiveValue - originalValue);
625    int iStatus;
626    if (solver->isProvenOptimal())
627        iStatus = 0; // optimal
628    else if (solver->isIterationLimitReached()
629             && !solver->isDualObjectiveLimitReached())
630        iStatus = 2; // unknown
631    else
632        iStatus = 1; // infeasible
633
634    bool feasible = iStatus != 1;
635    if (feasible) {
636#ifndef ZERO_SOS_TOLERANCE
637        double integerTolerance =
638          model_->getDblParam(CbcModel::CbcIntegerTolerance);
639#else
640        double integerTolerance = ZERO_SOS_TOLERANCE;
641#endif
642        const int * integerVariable = model_->integerVariable();
643        for (i = 0; i < numberIntegers; i++) {
644            int j = integerVariable[i];
645            double value = solution[j];
646            double nearest = floor(value + 0.5);
647            if (fabs(value - nearest) > integerTolerance)
648                unsatisfied++;
649        }
650    }
651    int way = branchingObject->way();
652    way = - way; // because after branch so moved on
653    double value = branchingObject->value();
654    CbcObjectUpdateData newData (this, way,
655                                 change, iStatus,
656                                 originalUnsatisfied - unsatisfied, value);
657    newData.originalObjective_ = originalValue;
658    // Solvers know about direction
659    double direction = solver->getObjSense();
660    solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
661    newData.cutoff_ *= direction;
662    return newData;
663}
664// Update object by CbcObjectUpdateData
665void
666CbcSOS::updateInformation(const CbcObjectUpdateData & data)
667{
668    bool feasible = data.status_ != 1;
669    int way = data.way_;
670    //double value = data.branchingValue_;
671    double originalValue = data.originalObjective_;
672    double change = data.change_;
673    if (way < 0) {
674        // down
675        if (!feasible) {
676            double distanceToCutoff = 0.0;
677            //double objectiveValue = model_->getCurrentMinimizationObjValue();
678            distanceToCutoff =  model_->getCutoff()  - originalValue;
679            if (distanceToCutoff < 1.0e20)
680                change = distanceToCutoff * 2.0;
681            else
682                change = (downDynamicPseudoRatio_ * shadowEstimateDown_ + 1.0e-3) * 10.0;
683        }
684        change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
685#ifdef PRINT_SHADOW
686        if (numberTimesDown_)
687            printf("Updating id %d - down change %g (true %g) - ndown %d estimated change %g - raw shadow estimate %g\n",
688                   id_, change, data.change_, numberTimesDown_, shadowEstimateDown_*
689                   (downDynamicPseudoRatio_ / ((double) numberTimesDown_)),
690                   shadowEstimateDown_);
691        else
692            printf("Updating id %d - down change %g (true %g) - shadow estimate %g\n",
693                   id_, change, data.change_, shadowEstimateDown_);
694#endif
695        numberTimesDown_++;
696        downDynamicPseudoRatio_ += change / shadowEstimateDown_;
697    } else {
698        // up
699        if (!feasible) {
700            double distanceToCutoff = 0.0;
701            //double objectiveValue = model_->getCurrentMinimizationObjValue();
702            distanceToCutoff =  model_->getCutoff()  - originalValue;
703            if (distanceToCutoff < 1.0e20)
704                change = distanceToCutoff * 2.0;
705            else
706                change = (upDynamicPseudoRatio_ * shadowEstimateUp_ + 1.0e-3) * 10.0;
707        }
708        change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
709#ifdef PRINT_SHADOW
710        if (numberTimesUp_)
711            printf("Updating id %d - up change %g (true %g) - nup %d estimated change %g - raw shadow estimate %g\n",
712                   id_, change, data.change_, numberTimesUp_, shadowEstimateUp_*
713                   (upDynamicPseudoRatio_ / ((double) numberTimesUp_)),
714                   shadowEstimateUp_);
715        else
716            printf("Updating id %d - up change %g (true %g) - shadow estimate %g\n",
717                   id_, change, data.change_, shadowEstimateUp_);
718#endif
719        numberTimesUp_++;
720        upDynamicPseudoRatio_ += change / shadowEstimateUp_;
721    }
722}
723
724/* Create an OsiSolverBranch object
725
726This returns NULL if branch not represented by bound changes
727*/
728OsiSolverBranch *
729CbcSOS::solverBranch() const
730{
731    int j;
732    const double * solution = model_->testSolution();
733#ifndef ZERO_SOS_TOLERANCE
734    double integerTolerance =
735        model_->getDblParam(CbcModel::CbcIntegerTolerance);
736#else
737    double integerTolerance = ZERO_SOS_TOLERANCE;
738#endif
739    OsiSolverInterface * solver = model_->solver();
740    const double * upper = solver->getColUpper();
741    int firstNonFixed = -1;
742    int lastNonFixed = -1;
743    int firstNonZero = -1;
744    int lastNonZero = -1;
745    double weight = 0.0;
746    double sum = 0.0;
747    double * fix = new double[numberMembers_];
748    int * which = new int[numberMembers_];
749    for (j = 0; j < numberMembers_; j++) {
750        int iColumn = members_[j];
751        // fix all on one side or other (even if fixed)
752        fix[j] = 0.0;
753        which[j] = iColumn;
754        if (upper[iColumn] || oddValues_) {
755            double value = CoinMax(0.0, solution[iColumn]);
756            sum += value;
757            if (firstNonFixed < 0)
758                firstNonFixed = j;
759            lastNonFixed = j;
760            if (value > integerTolerance) {
761                weight += weights_[j] * value;
762                if (firstNonZero < 0)
763                    firstNonZero = j;
764                lastNonZero = j;
765            }
766        }
767    }
768    assert (lastNonZero - firstNonZero >= sosType_) ;
769    // find where to branch
770    if (!oddValues_)
771      weight /= sum;
772    else
773      weight = 0.5*(weights_[firstNonZero]+weights_[lastNonZero]);
774    // down branch fixes ones above weight to 0
775    int iWhere;
776    int iDownStart = 0;
777    int iUpEnd = 0;
778    for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
779        if (weight < weights_[iWhere+1])
780            break;
781    if (sosType_ == 1) {
782        // SOS 1
783        iUpEnd = iWhere + 1;
784        iDownStart = iUpEnd;
785    } else {
786        // SOS 2
787        if (iWhere == firstNonFixed)
788            iWhere++;;
789        if (iWhere == lastNonFixed - 1)
790            iWhere = lastNonFixed - 2;
791        iUpEnd = iWhere + 1;
792        iDownStart = iUpEnd + 1;
793    }
794    //
795    OsiSolverBranch * branch = new OsiSolverBranch();
796    branch->addBranch(-1, 0, NULL, NULL, numberMembers_ - iDownStart, which + iDownStart, fix);
797    branch->addBranch(1, 0, NULL, NULL, iUpEnd, which, fix);
798    delete [] fix;
799    delete [] which;
800    return branch;
801}
802// Construct an OsiSOS object
803OsiSOS *
804CbcSOS::osiObject(const OsiSolverInterface * solver) const
805{
806    OsiSOS * obj = new OsiSOS(solver, numberMembers_, members_, weights_, sosType_);
807    obj->setPriority(priority());
808    return obj;
809}
810
811// Default Constructor
812CbcSOSBranchingObject::CbcSOSBranchingObject()
813        : CbcBranchingObject(),
814        firstNonzero_(-1),
815        lastNonzero_(-1)
816{
817    set_ = NULL;
818    separator_ = 0.0;
819}
820
821// Useful constructor
822CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
823        const CbcSOS * set,
824        int way ,
825        double separator)
826        : CbcBranchingObject(model, set->id(), way, 0.5)
827{
828    set_ = set;
829    separator_ = separator;
830    computeNonzeroRange();
831}
832
833// Copy constructor
834CbcSOSBranchingObject::CbcSOSBranchingObject (const CbcSOSBranchingObject & rhs)
835        : CbcBranchingObject(rhs),
836        firstNonzero_(rhs.firstNonzero_),
837        lastNonzero_(rhs.lastNonzero_)
838{
839    set_ = rhs.set_;
840    separator_ = rhs.separator_;
841}
842
843// Assignment operator
844CbcSOSBranchingObject &
845CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject & rhs)
846{
847    if (this != &rhs) {
848        CbcBranchingObject::operator=(rhs);
849        set_ = rhs.set_;
850        separator_ = rhs.separator_;
851        firstNonzero_ = rhs.firstNonzero_;
852        lastNonzero_ = rhs.lastNonzero_;
853    }
854    return *this;
855}
856CbcBranchingObject *
857CbcSOSBranchingObject::clone() const
858{
859    return (new CbcSOSBranchingObject(*this));
860}
861
862
863// Destructor
864CbcSOSBranchingObject::~CbcSOSBranchingObject ()
865{
866}
867
868void
869CbcSOSBranchingObject::computeNonzeroRange()
870{
871    const int numberMembers = set_->numberMembers();
872    const double * weights = set_->weights();
873    int i = 0;
874    if (way_ < 0) {
875        for ( i = 0; i < numberMembers; i++) {
876            if (weights[i] > separator_)
877                break;
878        }
879        assert (i < numberMembers);
880        firstNonzero_ = 0;
881        lastNonzero_ = i;
882    } else {
883        for ( i = 0; i < numberMembers; i++) {
884            if (weights[i] >= separator_)
885                break;
886        }
887        assert (i < numberMembers);
888        firstNonzero_ = i;
889        lastNonzero_ = numberMembers;
890    }
891}
892
893double
894CbcSOSBranchingObject::branch()
895{
896    decrementNumberBranchesLeft();
897    int numberMembers = set_->numberMembers();
898    const int * which = set_->members();
899    const double * weights = set_->weights();
900    OsiSolverInterface * solver = model_->solver();
901    const double * lower = solver->getColLower();
902    const double * upper = solver->getColUpper();
903    // *** for way - up means fix all those in down section
904    if (way_ < 0) {
905        int i;
906        for ( i = 0; i < numberMembers; i++) {
907            if (weights[i] > separator_)
908                break;
909        }
910        assert (i < numberMembers);
911        for (; i < numberMembers; i++) {
912            solver->setColLower(which[i], 0.0);
913            solver->setColUpper(which[i], 0.0);
914        }
915        way_ = 1;         // Swap direction
916    } else {
917        int i;
918        for ( i = 0; i < numberMembers; i++) {
919          if (weights[i] >= separator_) {
920                break;
921          } else {
922                solver->setColLower(which[i], 0.0);
923                solver->setColUpper(which[i], 0.0);
924          }
925        }
926        assert (i < numberMembers);
927        way_ = -1;        // Swap direction
928    }
929    computeNonzeroRange();
930    double predictedChange=0.0;
931    for (int i = 0; i < numberMembers; i++) {
932      int iColumn=which[i];
933      if (lower[iColumn]>upper[iColumn])
934        predictedChange=COIN_DBL_MAX;
935    }
936    return predictedChange;
937}
938/* Update bounds in solver as in 'branch' and update given bounds.
939   branchState is -1 for 'down' +1 for 'up' */
940void
941CbcSOSBranchingObject::fix(OsiSolverInterface * solver,
942                           double * lower, double * upper,
943                           int branchState) const
944{
945    int numberMembers = set_->numberMembers();
946    const int * which = set_->members();
947    const double * weights = set_->weights();
948    //const double * lower = solver->getColLower();
949    //const double * upper = solver->getColUpper();
950    // *** for way - up means fix all those in down section
951    if (branchState < 0) {
952        int i;
953        for ( i = 0; i < numberMembers; i++) {
954            if (weights[i] > separator_)
955                break;
956        }
957        assert (i < numberMembers);
958        for (; i < numberMembers; i++) {
959            solver->setColLower(which[i], 0.0);
960            lower[which[i]] = 0.0;
961            solver->setColUpper(which[i], 0.0);
962            upper[which[i]] = 0.0;
963        }
964    } else {
965        int i;
966        for ( i = 0; i < numberMembers; i++) {
967            if (weights[i] >= separator_) {
968                break;
969            } else {
970                solver->setColLower(which[i], 0.0);
971                lower[which[i]] = 0.0;
972                solver->setColUpper(which[i], 0.0);
973                upper[which[i]] = 0.0;
974            }
975        }
976        assert (i < numberMembers);
977    }
978}
979// Print what would happen
980void
981CbcSOSBranchingObject::print()
982{
983    int numberMembers = set_->numberMembers();
984    const int * which = set_->members();
985    const double * weights = set_->weights();
986    OsiSolverInterface * solver = model_->solver();
987    //const double * lower = solver->getColLower();
988    const double * upper = solver->getColUpper();
989    int first = numberMembers;
990    int last = -1;
991    int numberFixed = 0;
992    int numberOther = 0;
993    int i;
994    for ( i = 0; i < numberMembers; i++) {
995        double bound = upper[which[i]];
996        if (bound) {
997            first = CoinMin(first, i);
998            last = CoinMax(last, i);
999        }
1000    }
1001    // *** for way - up means fix all those in down section
1002    if (way_ < 0) {
1003        printf("SOS Down");
1004        for ( i = 0; i < numberMembers; i++) {
1005            double bound = upper[which[i]];
1006            if (weights[i] > separator_)
1007                break;
1008            else if (bound)
1009                numberOther++;
1010        }
1011        assert (i < numberMembers);
1012        for (; i < numberMembers; i++) {
1013            double bound = upper[which[i]];
1014            if (bound)
1015                numberFixed++;
1016        }
1017    } else {
1018        printf("SOS Up");
1019        for ( i = 0; i < numberMembers; i++) {
1020            double bound = upper[which[i]];
1021            if (weights[i] >= separator_)
1022                break;
1023            else if (bound)
1024                numberFixed++;
1025        }
1026        assert (i < numberMembers);
1027        for (; i < numberMembers; i++) {
1028            double bound = upper[which[i]];
1029            if (bound)
1030                numberOther++;
1031        }
1032    }
1033    printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
1034           separator_, which[first], weights[first], which[last], weights[last], numberFixed, numberOther);
1035}
1036
1037/** Compare the original object of \c this with the original object of \c
1038    brObj. Assumes that there is an ordering of the original objects.
1039    This method should be invoked only if \c this and brObj are of the same
1040    type.
1041    Return negative/0/positive depending on whether \c this is
1042    smaller/same/larger than the argument.
1043*/
1044int
1045CbcSOSBranchingObject::compareOriginalObject
1046(const CbcBranchingObject* brObj) const
1047{
1048    const CbcSOSBranchingObject* br =
1049        dynamic_cast<const CbcSOSBranchingObject*>(brObj);
1050    assert(br);
1051    const CbcSOS* s0 = set_;
1052    const CbcSOS* s1 = br->set_;
1053    if (s0->sosType() != s1->sosType()) {
1054        return s0->sosType() - s1->sosType();
1055    }
1056    if (s0->numberMembers() != s1->numberMembers()) {
1057        return s0->numberMembers() - s1->numberMembers();
1058    }
1059    const int memberCmp = memcmp(s0->members(), s1->members(),
1060                                 s0->numberMembers() * sizeof(int));
1061    if (memberCmp != 0) {
1062        return memberCmp;
1063    }
1064    return memcmp(s0->weights(), s1->weights(),
1065                  s0->numberMembers() * sizeof(double));
1066}
1067
1068/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1069    same type and must have the same original object, but they may have
1070    different feasible regions.
1071    Return the appropriate CbcRangeCompare value (first argument being the
1072    sub/superset if that's the case). In case of overlap (and if \c
1073    replaceIfOverlap is true) replace the current branching object with one
1074    whose feasible region is the overlap.
1075*/
1076CbcRangeCompare
1077CbcSOSBranchingObject::compareBranchingObject
1078(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1079{
1080    const CbcSOSBranchingObject* br =
1081        dynamic_cast<const CbcSOSBranchingObject*>(brObj);
1082    assert(br);
1083    if (firstNonzero_ < br->firstNonzero_) {
1084        if (lastNonzero_ >= br->lastNonzero_) {
1085            return CbcRangeSuperset;
1086        } else if (lastNonzero_ <= br->firstNonzero_) {
1087            return CbcRangeDisjoint;
1088        } else {
1089            // overlap
1090            if (replaceIfOverlap) {
1091                firstNonzero_ = br->firstNonzero_;
1092            }
1093            return CbcRangeOverlap;
1094        }
1095    } else if (firstNonzero_ > br->firstNonzero_) {
1096        if (lastNonzero_ <= br->lastNonzero_) {
1097            return CbcRangeSubset;
1098        } else if (firstNonzero_ >= br->lastNonzero_) {
1099            return CbcRangeDisjoint;
1100        } else {
1101            // overlap
1102            if (replaceIfOverlap) {
1103                lastNonzero_ = br->lastNonzero_;
1104            }
1105            return CbcRangeOverlap;
1106        }
1107    } else {
1108        if (lastNonzero_ == br->lastNonzero_) {
1109            return CbcRangeSame;
1110        }
1111        return lastNonzero_ < br->lastNonzero_ ? CbcRangeSubset : CbcRangeSuperset;
1112    }
1113    return CbcRangeSame; // fake return
1114}
1115
Note: See TracBrowser for help on using the repository browser.