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

Last change on this file since 2420 was 2420, checked in by forrest, 2 years ago

try and fix problem with tiny values

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