source: trunk/Cbc/src/CbcSOS.cpp

Last change on this file was 2491, checked in by forrest, 3 months ago

better SOS in mipstart, ctrl-c back, improve symmetric

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