source: trunk/Cbc/examples/driver5.cpp

Last change on this file was 2469, checked in by unxusr, 6 weeks ago

formatting

File size: 16.5 KB
Line 
1// $Id: driver5.cpp 2101 2014-12-03 17:43:20Z forrest $
2// Copyright (C) 2007, 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#include <cassert>
7#include <iomanip>
8
9#include "CoinPragma.hpp"
10#include "CbcModel.hpp"
11#include "OsiClpSolverInterface.hpp"
12#include "CbcSolver.hpp"
13#include "CbcSOS.hpp"
14
15#include "CoinTime.hpp"
16
17//#############################################################################
18
19/************************************************************************
20
21This main program shows how to take advantage of the standalone cbc in your program,
22while still making major modifications.
23First it reads in an integer model from an mps file
24Then it initializes the integer model with cbc defaults
25Then it calls CbcMain1 passing all parameters apart from first
26Then it modifies all SOS objects on first CbcEvent
27Then it solves
28Finally it prints solution
29
30************************************************************************/
31/*
32  This is obviously stupid - just to show what can be done -
33  just branch halfway rather than use weights and solution values
34 */
35class CbcUserSOS : public CbcSOS {
36
37public:
38  // Default Constructor
39  CbcUserSOS();
40
41  /* Constructor from CbcSOS
42     Just adds dummy integer
43  */
44
45  CbcUserSOS(CbcSOS object, int dummy);
46
47  // Copy constructor
48  CbcUserSOS(const CbcUserSOS &);
49
50  /// Clone
51  virtual CbcObject *clone() const;
52
53  // Assignment operator
54  CbcUserSOS &operator=(const CbcUserSOS &rhs);
55
56  // Destructor
57  virtual ~CbcUserSOS();
58
59  /// Infeasibility - large is 0.5
60  virtual double infeasibility(const OsiBranchingInformation *info,
61    int &preferredWay) const;
62
63  /// Creates a branching object
64  virtual CbcBranchingObject *createCbcBranch(OsiSolverInterface *solver, const OsiBranchingInformation *info, int way);
65
66private:
67  /// data
68
69  /// Dummy integer
70  int dummyInteger_;
71};
72//##############################################################################
73
74// Default Constructor
75CbcUserSOS::CbcUserSOS()
76  : CbcSOS()
77  , dummyInteger_(0)
78{
79}
80
81// Constructor from CbcSOS
82CbcUserSOS::CbcUserSOS(CbcSOS object, int dummyInteger)
83  : CbcSOS(object)
84  , dummyInteger_(dummyInteger)
85{
86}
87
88// Copy constructor
89CbcUserSOS::CbcUserSOS(const CbcUserSOS &rhs)
90  : CbcSOS(rhs)
91{
92  dummyInteger_ = rhs.dummyInteger_;
93}
94
95// Clone
96CbcObject *
97CbcUserSOS::clone() const
98{
99  return new CbcUserSOS(*this);
100}
101
102// Assignment operator
103CbcUserSOS &
104CbcUserSOS::operator=(const CbcUserSOS &rhs)
105{
106  if (this != &rhs) {
107    CbcSOS::operator=(rhs);
108    dummyInteger_ = rhs.dummyInteger_;
109  }
110  return *this;
111}
112
113// Destructor
114CbcUserSOS::~CbcUserSOS()
115{
116}
117/*
118  Routine to calculate standard infeasibility of an SOS set and return a
119  preferred branching direction.
120  This is just a copy of CbcSOS
121*/
122double
123CbcUserSOS::infeasibility(const OsiBranchingInformation *info,
124  int &preferredWay) const
125{
126  int j;
127  int firstNonZero = -1;
128  int lastNonZero = -1;
129  OsiSolverInterface *solver = model_->solver();
130  const double *solution = model_->testSolution();
131  const double *lower = solver->getColLower();
132  const double *upper = solver->getColUpper();
133  //double largestValue=0.0;
134#ifndef ZERO_ODD_TOLERANCE
135#define ZERO_SOS_TOLERANCE 1.0e-14
136#else
137#define ZERO_SOS_TOLERANCE ZERO_ODD_TOLERANCE
138#endif
139  //double integerTolerance =
140  //  model_->getDblParam(CbcModel::CbcIntegerTolerance);
141  double integerTolerance = ZERO_SOS_TOLERANCE;
142  double weight = 0.0;
143  double sum = 0.0;
144
145  // check bounds etc
146  double lastWeight = -1.0e100;
147  for (j = 0; j < numberMembers_; j++) {
148    int iColumn = members_[j];
149    /*
150          The value used here (1.0e-7) is larger than the value enforced in the
151          constructor.
152        */
153
154    if (lastWeight >= weights_[j] - 1.0e-7)
155      throw CoinError("Weights too close together in SOS", "infeasibility", "CbcSOS");
156    double value = CoinMax(lower[iColumn], solution[iColumn]);
157    value = CoinMin(upper[iColumn], value);
158    sum += value;
159    if (fabs(value) > integerTolerance && (upper[iColumn] > 0.0 || oddValues_)) {
160      weight += weights_[j] * value;
161      if (firstNonZero < 0)
162        firstNonZero = j;
163      lastNonZero = j;
164    }
165  }
166  /* ?? */
167  preferredWay = 1;
168  /*
169  SOS1 allows one nonzero; SOS2 allows two consecutive nonzeros. Infeasibility
170  is calculated as (.5)(range of nonzero values)/(number of members). So if
171  the first and last elements of the set are nonzero, we have maximum
172  infeasibility.
173*/
174  if (lastNonZero - firstNonZero >= sosType_) {
175    // find where to branch
176    if (!oddValues_)
177      weight /= sum;
178    else
179      weight = 0.5 * (weights_[firstNonZero] + weights_[lastNonZero]);
180    if (info->defaultDual_ >= 0.0 && info->usefulRegion_ && info->columnStart_) {
181      assert(sosType_ == 1);
182      int iWhere;
183      for (iWhere = firstNonZero; iWhere < lastNonZero - 1; iWhere++) {
184        if (weight < weights_[iWhere + 1]) {
185          break;
186        }
187      }
188      int jColumnDown = members_[iWhere];
189      int jColumnUp = members_[iWhere + 1];
190      int n = 0;
191      CoinBigIndex j;
192      double objMove = info->objective_[jColumnDown];
193      for (j = info->columnStart_[jColumnDown];
194           j < info->columnStart_[jColumnDown] + info->columnLength_[jColumnDown]; j++) {
195        double value = info->elementByColumn_[j];
196        int iRow = info->row_[j];
197        info->indexRegion_[n++] = iRow;
198        info->usefulRegion_[iRow] = value;
199      }
200      for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++) {
201        int jColumn = members_[iWhere];
202        double solValue = info->solution_[jColumn];
203        if (!solValue)
204          continue;
205        objMove -= info->objective_[jColumn] * solValue;
206        for (j = info->columnStart_[jColumn];
207             j < info->columnStart_[jColumn] + info->columnLength_[jColumn]; j++) {
208          double value = -info->elementByColumn_[j] * solValue;
209          int iRow = info->row_[j];
210          double oldValue = info->usefulRegion_[iRow];
211          if (!oldValue) {
212            info->indexRegion_[n++] = iRow;
213          } else {
214            value += oldValue;
215            if (!value)
216              value = 1.0e-100;
217          }
218          info->usefulRegion_[iRow] = value;
219        }
220      }
221      const double *pi = info->pi_;
222      const double *activity = info->rowActivity_;
223      const double *lower = info->rowLower_;
224      const double *upper = info->rowUpper_;
225      double tolerance = info->primalTolerance_;
226      double direction = info->direction_;
227      shadowEstimateDown_ = objMove * direction;
228      bool infeasible = false;
229      for (int k = 0; k < n; k++) {
230        int iRow = info->indexRegion_[k];
231        double movement = info->usefulRegion_[iRow];
232        // not this time info->usefulRegion_[iRow]=0.0;
233        double valueP = pi[iRow] * direction;
234        // if move makes infeasible then make at least default
235        double newValue = activity[iRow] + movement;
236        if (newValue > upper[iRow] + tolerance || newValue < lower[iRow] - tolerance) {
237          shadowEstimateDown_ += fabs(movement) * CoinMax(fabs(valueP), info->defaultDual_);
238          infeasible = true;
239        }
240      }
241      if (shadowEstimateDown_ < info->integerTolerance_) {
242        if (!infeasible) {
243          shadowEstimateDown_ = 1.0e-10;
244        } else
245          shadowEstimateDown_ = info->integerTolerance_;
246      }
247      // And other way
248      // take off
249      objMove -= info->objective_[jColumnDown];
250      for (j = info->columnStart_[jColumnDown];
251           j < info->columnStart_[jColumnDown] + info->columnLength_[jColumnDown]; j++) {
252        double value = -info->elementByColumn_[j];
253        int iRow = info->row_[j];
254        double oldValue = info->usefulRegion_[iRow];
255        if (!oldValue) {
256          info->indexRegion_[n++] = iRow;
257        } else {
258          value += oldValue;
259          if (!value)
260            value = 1.0e-100;
261        }
262        info->usefulRegion_[iRow] = value;
263      }
264      // add on
265      objMove += info->objective_[jColumnUp];
266      for (j = info->columnStart_[jColumnUp];
267           j < info->columnStart_[jColumnUp] + info->columnLength_[jColumnUp]; j++) {
268        double value = info->elementByColumn_[j];
269        int iRow = info->row_[j];
270        double oldValue = info->usefulRegion_[iRow];
271        if (!oldValue) {
272          info->indexRegion_[n++] = iRow;
273        } else {
274          value += oldValue;
275          if (!value)
276            value = 1.0e-100;
277        }
278        info->usefulRegion_[iRow] = value;
279      }
280      shadowEstimateUp_ = objMove * direction;
281      infeasible = false;
282      for (int k = 0; k < n; k++) {
283        int iRow = info->indexRegion_[k];
284        double movement = info->usefulRegion_[iRow];
285        info->usefulRegion_[iRow] = 0.0;
286        double valueP = pi[iRow] * direction;
287        // if move makes infeasible then make at least default
288        double newValue = activity[iRow] + movement;
289        if (newValue > upper[iRow] + tolerance || newValue < lower[iRow] - tolerance) {
290          shadowEstimateUp_ += fabs(movement) * CoinMax(fabs(valueP), info->defaultDual_);
291          infeasible = true;
292        }
293      }
294      if (shadowEstimateUp_ < info->integerTolerance_) {
295        if (!infeasible) {
296          shadowEstimateUp_ = 1.0e-10;
297        } else
298          shadowEstimateUp_ = info->integerTolerance_;
299      }
300      // adjust
301      double downCost = shadowEstimateDown_;
302      double upCost = shadowEstimateUp_;
303      if (numberTimesDown_)
304        downCost *= downDynamicPseudoRatio_ / static_cast< double >(numberTimesDown_);
305      if (numberTimesUp_)
306        upCost *= upDynamicPseudoRatio_ / static_cast< double >(numberTimesUp_);
307#define WEIGHT_AFTER 0.7
308#define WEIGHT_BEFORE 0.1
309      int stateOfSearch = model_->stateOfSearch() % 10;
310      double returnValue = 0.0;
311      double minValue = CoinMin(downCost, upCost);
312      double maxValue = CoinMax(downCost, upCost);
313      if (stateOfSearch <= 2) {
314        // no branching solution
315        returnValue = WEIGHT_BEFORE * minValue + (1.0 - WEIGHT_BEFORE) * maxValue;
316      } else {
317        returnValue = WEIGHT_AFTER * minValue + (1.0 - WEIGHT_AFTER) * maxValue;
318      }
319      return returnValue;
320    } else {
321      double value = lastNonZero - firstNonZero + 1;
322      value *= 0.5 / static_cast< double >(numberMembers_);
323      return value;
324    }
325  } else {
326    return 0.0; // satisfied
327  }
328}
329
330CbcBranchingObject *
331CbcUserSOS::createCbcBranch(OsiSolverInterface *solver, const OsiBranchingInformation * /*info*/, int way)
332{
333  int j;
334  const double *solution = model_->testSolution();
335#ifndef ZERO_SOS_TOLERANCE
336  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
337#else
338  double integerTolerance = ZERO_SOS_TOLERANCE;
339#endif
340  const double *lower = solver->getColLower();
341  const double *upper = solver->getColUpper();
342  int firstNonZero = -1;
343  int lastNonZero = -1;
344  for (j = 0; j < numberMembers_; j++) {
345    int iColumn = members_[j];
346    double value = CoinMax(lower[iColumn], solution[iColumn]);
347    value = CoinMin(upper[iColumn], value);
348    if (fabs(value) > integerTolerance) {
349      if (firstNonZero < 0)
350        firstNonZero = j;
351      lastNonZero = j;
352    }
353  }
354  assert(lastNonZero - firstNonZero >= sosType_);
355  // find where to branch
356  // stupid but just to show
357  double separator;
358  int iWhere = (firstNonZero + lastNonZero) / 2;
359  if (sosType_ == 1) {
360    // SOS 1
361    separator = 0.5 * (weights_[iWhere] + weights_[iWhere + 1]);
362  } else {
363    // SOS 2
364    if (iWhere + 1 == lastNonZero)
365      iWhere--;
366    separator = weights_[iWhere + 1];
367  }
368  // create object
369  CbcBranchingObject *branch;
370  branch = new CbcSOSBranchingObject(model_, this, way, separator);
371  branch->setOriginalObject(this);
372  return branch;
373}
374#include "CbcEventHandler.hpp"
375/** This is so user can trap events and do useful stuff. 
376
377    CbcModel model_ is available as well as anything else you care
378    to pass in
379*/
380
381class MyEventHandler3 : public CbcEventHandler {
382
383public:
384  /**@name Overrides */
385  //@{
386  virtual CbcAction event(CbcEvent whichEvent);
387  //@}
388
389  /**@name Constructors, destructor etc*/
390  //@{
391  /** Default constructor. */
392  MyEventHandler3();
393  /// Constructor with pointer to model (redundant as setEventHandler does)
394  MyEventHandler3(CbcModel *model);
395  /** Destructor */
396  virtual ~MyEventHandler3();
397  /** The copy constructor. */
398  MyEventHandler3(const MyEventHandler3 &rhs);
399  /// Assignment
400  MyEventHandler3 &operator=(const MyEventHandler3 &rhs);
401  /// Clone
402  virtual CbcEventHandler *clone() const;
403  //@}
404
405protected:
406  // data goes here
407};
408//-------------------------------------------------------------------
409// Default Constructor
410//-------------------------------------------------------------------
411MyEventHandler3::MyEventHandler3()
412  : CbcEventHandler()
413{
414}
415
416//-------------------------------------------------------------------
417// Copy constructor
418//-------------------------------------------------------------------
419MyEventHandler3::MyEventHandler3(const MyEventHandler3 &rhs)
420  : CbcEventHandler(rhs)
421{
422}
423
424// Constructor with pointer to model
425MyEventHandler3::MyEventHandler3(CbcModel *model)
426  : CbcEventHandler(model)
427{
428}
429
430//-------------------------------------------------------------------
431// Destructor
432//-------------------------------------------------------------------
433MyEventHandler3::~MyEventHandler3()
434{
435}
436
437//----------------------------------------------------------------
438// Assignment operator
439//-------------------------------------------------------------------
440MyEventHandler3 &
441MyEventHandler3::operator=(const MyEventHandler3 &rhs)
442{
443  if (this != &rhs) {
444    CbcEventHandler::operator=(rhs);
445  }
446  return *this;
447}
448//-------------------------------------------------------------------
449// Clone
450//-------------------------------------------------------------------
451CbcEventHandler *MyEventHandler3::clone() const
452{
453  return new MyEventHandler3(*this);
454}
455static bool firstTime = true;
456CbcEventHandler::CbcAction
457MyEventHandler3::event(CbcEvent whichEvent)
458{
459  // If in sub tree carry on
460  if (!model_->parentModel()) {
461    if (firstTime) {
462      // Replace all SOS with home grown version
463      // Could of course do same for integers
464      firstTime = false;
465      int numberObjects = model_->numberObjects();
466      OsiObject **objects = model_->objects();
467      for (int i = 0; i < numberObjects; i++) {
468        CbcSOS *sosObj = dynamic_cast< CbcSOS * >(objects[i]);
469        if (sosObj) {
470          objects[i] = new CbcUserSOS(*sosObj, i);
471          delete sosObj;
472        }
473      }
474    }
475  }
476  return noAction; // carry on
477}
478
479int main(int argc, const char *argv[])
480{
481
482  OsiClpSolverInterface solver1;
483  // Read in model using argv[1]
484  // and assert that it is a clean model
485  if (argc < 2) {
486    fprintf(stderr, "Need input file name.\n");
487    exit(1);
488  }
489  int numReadErrors;
490  if (!strstr(argv[1], ".lp"))
491    numReadErrors = solver1.readMps(argv[1], "");
492  else
493    numReadErrors = solver1.readLp(argv[1], 1.0e-10);
494  if (numReadErrors != 0) {
495    printf("%d errors reading file\n", numReadErrors);
496    return numReadErrors;
497  }
498
499  // Messy code below copied from CbcSolver.cpp
500  // Pass to Cbc initialize defaults
501  CbcModel modelA(solver1);
502  CbcModel *model = &modelA;
503  CbcMain0(modelA);
504  // Event handler
505  MyEventHandler3 eventHandler;
506  model->passInEventHandler(&eventHandler);
507  /* Now go into code for standalone solver
508     Could copy arguments and add -quit at end to be safe
509     but this will do
510  */
511  if (argc > 2) {
512    CbcMain1(argc - 1, argv + 1, modelA);
513  } else {
514    const char *argv2[] = { "driver5", "-solve", "-quit" };
515    CbcMain1(3, argv2, modelA);
516  }
517  // Solver was cloned so get current copy
518  OsiSolverInterface *solver = model->solver();
519  // Print solution if finished (could get from model->bestSolution() as well
520
521  if (model->bestSolution()) {
522
523    const double *solution = solver->getColSolution();
524
525    int iColumn;
526    int numberColumns = solver->getNumCols();
527    std::cout << std::setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14);
528
529    std::cout << "--------------------------------------" << std::endl;
530    // names may not be in current solver - use original
531
532    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
533      double value = solution[iColumn];
534      if (fabs(value) > 1.0e-7)
535        std::cout << std::setw(6) << iColumn << " " << std::setw(8) << setiosflags(std::ios::left) << solver1.getModelPtr()->columnName(iColumn)
536                  << resetiosflags(std::ios::adjustfield) << std::setw(14) << " " << value << std::endl;
537    }
538    std::cout << "--------------------------------------" << std::endl;
539
540    std::cout << std::resetiosflags(std::ios::fixed | std::ios::showpoint | std::ios::scientific);
541  } else {
542    std::cout << " No solution!" << std::endl;
543  }
544  return 0;
545}
Note: See TracBrowser for help on using the repository browser.