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

Last change on this file since 2070 was 2070, checked in by forrest, 4 years ago

allow sos in lp files and fix odd SOS branching

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