source: stable/2.8/Cbc/src/CbcSOS.cpp @ 1942

Last change on this file since 1942 was 1902, checked in by stefan, 6 years ago

sync with trunk rev 1901

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