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

Last change on this file since 1826 was 1799, checked in by forrest, 7 years ago

more logical comparison and out assert on pi for sos

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