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

Last change on this file since 2055 was 2055, checked in by forrest, 5 years ago

Fix SOS problem and also bad bestSolution_

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 35.7 KB
Line 
1// $Id: CbcSOS.cpp 2055 2014-08-09 16:05:41Z 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{
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    // Might get here in odd situation if so fix all
464    if (lastNonZero - firstNonZero < sosType_) {
465      for (j = 0; j < firstNonZero; j++) {
466        int iColumn = members_[j];
467        solver->setColUpper(iColumn, 0.0);
468      }
469      for (j = lastNonZero + 1; j < numberMembers_; j++) {
470        int iColumn = members_[j];
471        solver->setColUpper(iColumn, 0.0);
472      }
473    } else {
474      for (j = 0; j < numberMembers_; j++) {
475        int iColumn = members_[j];
476        solver->setColUpper(iColumn, 0.0);
477        solver->setColLower(iColumn, 1.0);
478      }
479    }
480}
481// Redoes data when sequence numbers change
482void
483CbcSOS::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
484{
485    model_ = model;
486    int n2 = 0;
487    for (int j = 0; j < numberMembers_; j++) {
488        int iColumn = members_[j];
489        int i;
490        for (i = 0; i < numberColumns; i++) {
491            if (originalColumns[i] == iColumn)
492                break;
493        }
494        if (i < numberColumns) {
495            members_[n2] = i;
496            weights_[n2++] = weights_[j];
497        }
498    }
499    if (n2 < numberMembers_) {
500      COIN_DETAIL_PRINT(printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2));
501        numberMembers_ = n2;
502    }
503}
504CbcBranchingObject *
505CbcSOS::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way)
506{
507    int j;
508    const double * solution = model_->testSolution();
509    double integerTolerance =
510        model_->getDblParam(CbcModel::CbcIntegerTolerance);
511    //OsiSolverInterface * solver = model_->solver();
512    const double * upper = solver->getColUpper();
513    int firstNonFixed = -1;
514    int lastNonFixed = -1;
515    int firstNonZero = -1;
516    int lastNonZero = -1;
517    double weight = 0.0;
518    double sum = 0.0;
519    for (j = 0; j < numberMembers_; j++) {
520        int iColumn = members_[j];
521        if (upper[iColumn]) {
522            double value = CoinMax(0.0, solution[iColumn]);
523            sum += value;
524            if (firstNonFixed < 0)
525                firstNonFixed = j;
526            lastNonFixed = j;
527            if (value > integerTolerance) {
528                weight += weights_[j] * value;
529                if (firstNonZero < 0)
530                    firstNonZero = j;
531                lastNonZero = j;
532            }
533        }
534    }
535    assert (lastNonZero - firstNonZero >= sosType_) ;
536    // find where to branch
537    assert (sum > 0.0);
538    weight /= sum;
539    int iWhere;
540    double separator = 0.0;
541    for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
542        if (weight < weights_[iWhere+1])
543            break;
544    if (sosType_ == 1) {
545        // SOS 1
546        separator = 0.5 * (weights_[iWhere] + weights_[iWhere+1]);
547    } else {
548        // SOS 2
549        if (iWhere == firstNonFixed)
550            iWhere++;;
551        if (iWhere == lastNonFixed - 1)
552            iWhere = lastNonFixed - 2;
553        separator = weights_[iWhere+1];
554    }
555    // create object
556    CbcBranchingObject * branch;
557    branch = new CbcSOSBranchingObject(model_, this, way, separator);
558    branch->setOriginalObject(this);
559    return branch;
560}
561/* Pass in information on branch just done and create CbcObjectUpdateData instance.
562   If object does not need data then backward pointer will be NULL.
563   Assumes can get information from solver */
564CbcObjectUpdateData
565CbcSOS::createUpdateInformation(const OsiSolverInterface * solver,
566                                const CbcNode * node,
567                                const CbcBranchingObject * branchingObject)
568{
569    double originalValue = node->objectiveValue();
570    int originalUnsatisfied = node->numberUnsatisfied();
571    double objectiveValue = solver->getObjValue() * solver->getObjSense();
572    int unsatisfied = 0;
573    int i;
574    //might be base model - doesn't matter
575    int numberIntegers = model_->numberIntegers();;
576    const double * solution = solver->getColSolution();
577    //const double * lower = solver->getColLower();
578    //const double * upper = solver->getColUpper();
579    double change = CoinMax(0.0, objectiveValue - originalValue);
580    int iStatus;
581    if (solver->isProvenOptimal())
582        iStatus = 0; // optimal
583    else if (solver->isIterationLimitReached()
584             && !solver->isDualObjectiveLimitReached())
585        iStatus = 2; // unknown
586    else
587        iStatus = 1; // infeasible
588
589    bool feasible = iStatus != 1;
590    if (feasible) {
591        double integerTolerance =
592            model_->getDblParam(CbcModel::CbcIntegerTolerance);
593        const int * integerVariable = model_->integerVariable();
594        for (i = 0; i < numberIntegers; i++) {
595            int j = integerVariable[i];
596            double value = solution[j];
597            double nearest = floor(value + 0.5);
598            if (fabs(value - nearest) > integerTolerance)
599                unsatisfied++;
600        }
601    }
602    int way = branchingObject->way();
603    way = - way; // because after branch so moved on
604    double value = branchingObject->value();
605    CbcObjectUpdateData newData (this, way,
606                                 change, iStatus,
607                                 originalUnsatisfied - unsatisfied, value);
608    newData.originalObjective_ = originalValue;
609    // Solvers know about direction
610    double direction = solver->getObjSense();
611    solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
612    newData.cutoff_ *= direction;
613    return newData;
614}
615// Update object by CbcObjectUpdateData
616void
617CbcSOS::updateInformation(const CbcObjectUpdateData & data)
618{
619    bool feasible = data.status_ != 1;
620    int way = data.way_;
621    //double value = data.branchingValue_;
622    double originalValue = data.originalObjective_;
623    double change = data.change_;
624    if (way < 0) {
625        // down
626        if (!feasible) {
627            double distanceToCutoff = 0.0;
628            //double objectiveValue = model_->getCurrentMinimizationObjValue();
629            distanceToCutoff =  model_->getCutoff()  - originalValue;
630            if (distanceToCutoff < 1.0e20)
631                change = distanceToCutoff * 2.0;
632            else
633                change = (downDynamicPseudoRatio_ * shadowEstimateDown_ + 1.0e-3) * 10.0;
634        }
635        change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
636#ifdef PRINT_SHADOW
637        if (numberTimesDown_)
638            printf("Updating id %d - down change %g (true %g) - ndown %d estimated change %g - raw shadow estimate %g\n",
639                   id_, change, data.change_, numberTimesDown_, shadowEstimateDown_*
640                   (downDynamicPseudoRatio_ / ((double) numberTimesDown_)),
641                   shadowEstimateDown_);
642        else
643            printf("Updating id %d - down change %g (true %g) - shadow estimate %g\n",
644                   id_, change, data.change_, shadowEstimateDown_);
645#endif
646        numberTimesDown_++;
647        downDynamicPseudoRatio_ += change / shadowEstimateDown_;
648    } else {
649        // up
650        if (!feasible) {
651            double distanceToCutoff = 0.0;
652            //double objectiveValue = model_->getCurrentMinimizationObjValue();
653            distanceToCutoff =  model_->getCutoff()  - originalValue;
654            if (distanceToCutoff < 1.0e20)
655                change = distanceToCutoff * 2.0;
656            else
657                change = (upDynamicPseudoRatio_ * shadowEstimateUp_ + 1.0e-3) * 10.0;
658        }
659        change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
660#ifdef PRINT_SHADOW
661        if (numberTimesUp_)
662            printf("Updating id %d - up change %g (true %g) - nup %d estimated change %g - raw shadow estimate %g\n",
663                   id_, change, data.change_, numberTimesUp_, shadowEstimateUp_*
664                   (upDynamicPseudoRatio_ / ((double) numberTimesUp_)),
665                   shadowEstimateUp_);
666        else
667            printf("Updating id %d - up change %g (true %g) - shadow estimate %g\n",
668                   id_, change, data.change_, shadowEstimateUp_);
669#endif
670        numberTimesUp_++;
671        upDynamicPseudoRatio_ += change / shadowEstimateUp_;
672    }
673}
674
675/* Create an OsiSolverBranch object
676
677This returns NULL if branch not represented by bound changes
678*/
679OsiSolverBranch *
680CbcSOS::solverBranch() const
681{
682    int j;
683    const double * solution = model_->testSolution();
684    double integerTolerance =
685        model_->getDblParam(CbcModel::CbcIntegerTolerance);
686    OsiSolverInterface * solver = model_->solver();
687    const double * upper = solver->getColUpper();
688    int firstNonFixed = -1;
689    int lastNonFixed = -1;
690    int firstNonZero = -1;
691    int lastNonZero = -1;
692    double weight = 0.0;
693    double sum = 0.0;
694    double * fix = new double[numberMembers_];
695    int * which = new int[numberMembers_];
696    for (j = 0; j < numberMembers_; j++) {
697        int iColumn = members_[j];
698        // fix all on one side or other (even if fixed)
699        fix[j] = 0.0;
700        which[j] = iColumn;
701        if (upper[iColumn]) {
702            double value = CoinMax(0.0, solution[iColumn]);
703            sum += value;
704            if (firstNonFixed < 0)
705                firstNonFixed = j;
706            lastNonFixed = j;
707            if (value > integerTolerance) {
708                weight += weights_[j] * value;
709                if (firstNonZero < 0)
710                    firstNonZero = j;
711                lastNonZero = j;
712            }
713        }
714    }
715    assert (lastNonZero - firstNonZero >= sosType_) ;
716    // find where to branch
717    assert (sum > 0.0);
718    weight /= sum;
719    // down branch fixes ones above weight to 0
720    int iWhere;
721    int iDownStart = 0;
722    int iUpEnd = 0;
723    for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
724        if (weight < weights_[iWhere+1])
725            break;
726    if (sosType_ == 1) {
727        // SOS 1
728        iUpEnd = iWhere + 1;
729        iDownStart = iUpEnd;
730    } else {
731        // SOS 2
732        if (iWhere == firstNonFixed)
733            iWhere++;;
734        if (iWhere == lastNonFixed - 1)
735            iWhere = lastNonFixed - 2;
736        iUpEnd = iWhere + 1;
737        iDownStart = iUpEnd + 1;
738    }
739    //
740    OsiSolverBranch * branch = new OsiSolverBranch();
741    branch->addBranch(-1, 0, NULL, NULL, numberMembers_ - iDownStart, which + iDownStart, fix);
742    branch->addBranch(1, 0, NULL, NULL, iUpEnd, which, fix);
743    delete [] fix;
744    delete [] which;
745    return branch;
746}
747// Construct an OsiSOS object
748OsiSOS *
749CbcSOS::osiObject(const OsiSolverInterface * solver) const
750{
751    OsiSOS * obj = new OsiSOS(solver, numberMembers_, members_, weights_, sosType_);
752    obj->setPriority(priority());
753    return obj;
754}
755
756// Default Constructor
757CbcSOSBranchingObject::CbcSOSBranchingObject()
758        : CbcBranchingObject(),
759        firstNonzero_(-1),
760        lastNonzero_(-1)
761{
762    set_ = NULL;
763    separator_ = 0.0;
764}
765
766// Useful constructor
767CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
768        const CbcSOS * set,
769        int way ,
770        double separator)
771        : CbcBranchingObject(model, set->id(), way, 0.5)
772{
773    set_ = set;
774    separator_ = separator;
775    computeNonzeroRange();
776}
777
778// Copy constructor
779CbcSOSBranchingObject::CbcSOSBranchingObject (const CbcSOSBranchingObject & rhs)
780        : CbcBranchingObject(rhs),
781        firstNonzero_(rhs.firstNonzero_),
782        lastNonzero_(rhs.lastNonzero_)
783{
784    set_ = rhs.set_;
785    separator_ = rhs.separator_;
786}
787
788// Assignment operator
789CbcSOSBranchingObject &
790CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject & rhs)
791{
792    if (this != &rhs) {
793        CbcBranchingObject::operator=(rhs);
794        set_ = rhs.set_;
795        separator_ = rhs.separator_;
796        firstNonzero_ = rhs.firstNonzero_;
797        lastNonzero_ = rhs.lastNonzero_;
798    }
799    return *this;
800}
801CbcBranchingObject *
802CbcSOSBranchingObject::clone() const
803{
804    return (new CbcSOSBranchingObject(*this));
805}
806
807
808// Destructor
809CbcSOSBranchingObject::~CbcSOSBranchingObject ()
810{
811}
812
813void
814CbcSOSBranchingObject::computeNonzeroRange()
815{
816    const int numberMembers = set_->numberMembers();
817    const double * weights = set_->weights();
818    int i = 0;
819    if (way_ < 0) {
820        for ( i = 0; i < numberMembers; i++) {
821            if (weights[i] > separator_)
822                break;
823        }
824        assert (i < numberMembers);
825        firstNonzero_ = 0;
826        lastNonzero_ = i;
827    } else {
828        for ( i = 0; i < numberMembers; i++) {
829            if (weights[i] >= separator_)
830                break;
831        }
832        assert (i < numberMembers);
833        firstNonzero_ = i;
834        lastNonzero_ = numberMembers;
835    }
836}
837
838double
839CbcSOSBranchingObject::branch()
840{
841    decrementNumberBranchesLeft();
842    int numberMembers = set_->numberMembers();
843    const int * which = set_->members();
844    const double * weights = set_->weights();
845    OsiSolverInterface * solver = model_->solver();
846    const double * lower = solver->getColLower();
847    const double * upper = solver->getColUpper();
848    // *** for way - up means fix all those in down section
849    if (way_ < 0) {
850        int i;
851        for ( i = 0; i < numberMembers; i++) {
852            if (weights[i] > separator_)
853                break;
854        }
855        assert (i < numberMembers);
856        for (; i < numberMembers; i++)
857            solver->setColUpper(which[i], 0.0);
858        way_ = 1;         // Swap direction
859    } else {
860        int i;
861        for ( i = 0; i < numberMembers; i++) {
862            if (weights[i] >= separator_)
863                break;
864            else
865                solver->setColUpper(which[i], 0.0);
866        }
867        assert (i < numberMembers);
868        way_ = -1;        // Swap direction
869    }
870    computeNonzeroRange();
871    double predictedChange=0.0;
872    for (int i = 0; i < numberMembers; i++) {
873      int iColumn=which[i];
874      if (lower[iColumn]>upper[iColumn])
875        predictedChange=COIN_DBL_MAX;
876    }
877    return predictedChange;
878}
879/* Update bounds in solver as in 'branch' and update given bounds.
880   branchState is -1 for 'down' +1 for 'up' */
881void
882CbcSOSBranchingObject::fix(OsiSolverInterface * solver,
883                           double * /*lower*/, double * upper,
884                           int branchState) const
885{
886    int numberMembers = set_->numberMembers();
887    const int * which = set_->members();
888    const double * weights = set_->weights();
889    //const double * lower = solver->getColLower();
890    //const double * upper = solver->getColUpper();
891    // *** for way - up means fix all those in down section
892    if (branchState < 0) {
893        int i;
894        for ( i = 0; i < numberMembers; i++) {
895            if (weights[i] > separator_)
896                break;
897        }
898        assert (i < numberMembers);
899        for (; i < numberMembers; i++) {
900            solver->setColUpper(which[i], 0.0);
901            upper[which[i]] = 0.0;
902        }
903    } else {
904        int i;
905        for ( i = 0; i < numberMembers; i++) {
906            if (weights[i] >= separator_) {
907                break;
908            } else {
909                solver->setColUpper(which[i], 0.0);
910                upper[which[i]] = 0.0;
911            }
912        }
913        assert (i < numberMembers);
914    }
915}
916// Print what would happen
917void
918CbcSOSBranchingObject::print()
919{
920    int numberMembers = set_->numberMembers();
921    const int * which = set_->members();
922    const double * weights = set_->weights();
923    OsiSolverInterface * solver = model_->solver();
924    //const double * lower = solver->getColLower();
925    const double * upper = solver->getColUpper();
926    int first = numberMembers;
927    int last = -1;
928    int numberFixed = 0;
929    int numberOther = 0;
930    int i;
931    for ( i = 0; i < numberMembers; i++) {
932        double bound = upper[which[i]];
933        if (bound) {
934            first = CoinMin(first, i);
935            last = CoinMax(last, i);
936        }
937    }
938    // *** for way - up means fix all those in down section
939    if (way_ < 0) {
940        printf("SOS Down");
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                numberOther++;
947        }
948        assert (i < numberMembers);
949        for (; i < numberMembers; i++) {
950            double bound = upper[which[i]];
951            if (bound)
952                numberFixed++;
953        }
954    } else {
955        printf("SOS Up");
956        for ( i = 0; i < numberMembers; i++) {
957            double bound = upper[which[i]];
958            if (weights[i] >= separator_)
959                break;
960            else if (bound)
961                numberFixed++;
962        }
963        assert (i < numberMembers);
964        for (; i < numberMembers; i++) {
965            double bound = upper[which[i]];
966            if (bound)
967                numberOther++;
968        }
969    }
970    printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
971           separator_, which[first], weights[first], which[last], weights[last], numberFixed, numberOther);
972}
973
974/** Compare the original object of \c this with the original object of \c
975    brObj. Assumes that there is an ordering of the original objects.
976    This method should be invoked only if \c this and brObj are of the same
977    type.
978    Return negative/0/positive depending on whether \c this is
979    smaller/same/larger than the argument.
980*/
981int
982CbcSOSBranchingObject::compareOriginalObject
983(const CbcBranchingObject* brObj) const
984{
985    const CbcSOSBranchingObject* br =
986        dynamic_cast<const CbcSOSBranchingObject*>(brObj);
987    assert(br);
988    const CbcSOS* s0 = set_;
989    const CbcSOS* s1 = br->set_;
990    if (s0->sosType() != s1->sosType()) {
991        return s0->sosType() - s1->sosType();
992    }
993    if (s0->numberMembers() != s1->numberMembers()) {
994        return s0->numberMembers() - s1->numberMembers();
995    }
996    const int memberCmp = memcmp(s0->members(), s1->members(),
997                                 s0->numberMembers() * sizeof(int));
998    if (memberCmp != 0) {
999        return memberCmp;
1000    }
1001    return memcmp(s0->weights(), s1->weights(),
1002                  s0->numberMembers() * sizeof(double));
1003}
1004
1005/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1006    same type and must have the same original object, but they may have
1007    different feasible regions.
1008    Return the appropriate CbcRangeCompare value (first argument being the
1009    sub/superset if that's the case). In case of overlap (and if \c
1010    replaceIfOverlap is true) replace the current branching object with one
1011    whose feasible region is the overlap.
1012*/
1013CbcRangeCompare
1014CbcSOSBranchingObject::compareBranchingObject
1015(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1016{
1017    const CbcSOSBranchingObject* br =
1018        dynamic_cast<const CbcSOSBranchingObject*>(brObj);
1019    assert(br);
1020    if (firstNonzero_ < br->firstNonzero_) {
1021        if (lastNonzero_ >= br->lastNonzero_) {
1022            return CbcRangeSuperset;
1023        } else if (lastNonzero_ <= br->firstNonzero_) {
1024            return CbcRangeDisjoint;
1025        } else {
1026            // overlap
1027            if (replaceIfOverlap) {
1028                firstNonzero_ = br->firstNonzero_;
1029            }
1030            return CbcRangeOverlap;
1031        }
1032    } else if (firstNonzero_ > br->firstNonzero_) {
1033        if (lastNonzero_ <= br->lastNonzero_) {
1034            return CbcRangeSubset;
1035        } else if (firstNonzero_ >= br->lastNonzero_) {
1036            return CbcRangeDisjoint;
1037        } else {
1038            // overlap
1039            if (replaceIfOverlap) {
1040                lastNonzero_ = br->lastNonzero_;
1041            }
1042            return CbcRangeOverlap;
1043        }
1044    } else {
1045        if (lastNonzero_ == br->lastNonzero_) {
1046            return CbcRangeSame;
1047        }
1048        return lastNonzero_ < br->lastNonzero_ ? CbcRangeSubset : CbcRangeSuperset;
1049    }
1050    return CbcRangeSame; // fake return
1051}
1052
Note: See TracBrowser for help on using the repository browser.