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

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

I am a stupid idiot id_ was in so don't need setNumber_

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