source: trunk/Samples/CbcBranchLink.cpp @ 134

Last change on this file since 134 was 134, checked in by forrest, 15 years ago

to match CbcNode?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.0 KB
Line 
1// Copyright (C) 2005, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11
12#include "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcMessage.hpp"
15#include "CbcBranchLink.hpp"
16#include "CoinError.hpp"
17
18// Default Constructor
19CbcLink::CbcLink ()
20  : CbcObject(),
21    weights_(NULL),
22    numberMembers_(0),
23    numberLinks_(0),
24    first_(-1)
25{
26}
27
28// Useful constructor (which are indices)
29CbcLink::CbcLink (CbcModel * model,  int numberMembers,
30           int numberLinks, int first , const double * weights, int identifier)
31  : CbcObject(model),
32    numberMembers_(numberMembers),
33    numberLinks_(numberLinks),
34    first_(first)
35{
36  id_=identifier;
37  if (numberMembers_) {
38    weights_ = new double[numberMembers_];
39    if (weights) {
40      memcpy(weights_,weights,numberMembers_*sizeof(double));
41    } else {
42      for (int i=0;i<numberMembers_;i++)
43        weights_[i]=i;
44    }
45    // weights must be increasing
46    int i;
47    double last=-COIN_DBL_MAX;
48    for (i=0;i<numberMembers_;i++) {
49      assert (weights_[i]>last+1.0e-12);
50      last=weights_[i];
51    }
52  } else {
53    weights_ = NULL;
54  }
55}
56
57// Copy constructor
58CbcLink::CbcLink ( const CbcLink & rhs)
59  :CbcObject(rhs)
60{
61  numberMembers_ = rhs.numberMembers_;
62  numberLinks_ = rhs.numberLinks_;
63  first_ = rhs.first_;
64  if (numberMembers_) {
65    weights_ = new double[numberMembers_];
66    memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
67  } else {
68    weights_ = NULL;
69  }
70}
71
72// Clone
73CbcObject *
74CbcLink::clone() const
75{
76  return new CbcLink(*this);
77}
78
79// Assignment operator
80CbcLink & 
81CbcLink::operator=( const CbcLink& rhs)
82{
83  if (this!=&rhs) {
84    CbcObject::operator=(rhs);
85    delete [] weights_;
86    numberMembers_ = rhs.numberMembers_;
87    numberLinks_ = rhs.numberLinks_;
88    first_ = rhs.first_;
89    if (numberMembers_) {
90      weights_ = new double[numberMembers_];
91      memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
92    } else {
93      weights_ = NULL;
94    }
95  }
96  return *this;
97}
98
99// Destructor
100CbcLink::~CbcLink ()
101{
102  delete [] weights_;
103}
104
105// Infeasibility - large is 0.5
106double 
107CbcLink::infeasibility(int & preferredWay) const
108{
109  int j;
110  int firstNonZero=-1;
111  int lastNonZero = -1;
112  OsiSolverInterface * solver = model_->solver();
113  const double * solution = model_->testSolution();
114  const double * lower = solver->getColLower();
115  const double * upper = solver->getColUpper();
116  double integerTolerance = 
117    model_->getDblParam(CbcModel::CbcIntegerTolerance);
118  double weight = 0.0;
119  double sum =0.0;
120
121  // check bounds etc
122  double lastWeight=-1.0e100;
123  int base=first_;
124  for (j=0;j<numberMembers_;j++) {
125    for (int k=0;k<numberLinks_;k++) {
126      int iColumn = base+k;
127      if (lower[iColumn])
128        throw CoinError("Non zero lower bound in CBCLink","infeasibility","CbcLink");
129      if (lastWeight>=weights_[j]-1.0e-7)
130        throw CoinError("Weights too close together in CBCLink","infeasibility","CbcLink");
131      double value = CoinMax(0.0,solution[iColumn]);
132      sum += value;
133      if (value>integerTolerance&&upper[iColumn]) {
134        // Possibly due to scaling a fixed variable might slip through
135        if (value>upper[iColumn]+1.0e-8) {
136          // Could change to #ifdef CBC_DEBUG
137#ifndef NDEBUG
138          if (model_->messageHandler()->logLevel()>1)
139            printf("** Variable %d (%d) has value %g and upper bound of %g\n",
140                   iColumn,j,value,upper[iColumn]);
141#endif
142        } 
143        weight += weights_[j]*CoinMin(value,upper[iColumn]);
144        if (firstNonZero<0)
145          firstNonZero=j;
146        lastNonZero=j;
147      }
148    }
149    base += numberLinks_;
150  }
151  preferredWay=1;
152  if (lastNonZero-firstNonZero>=1) {
153    // find where to branch
154    assert (sum>0.0);
155    weight /= sum;
156    double value = lastNonZero-firstNonZero+1;
157    value *= 0.5/((double) numberMembers_);
158    return value;
159  } else {
160    return 0.0; // satisfied
161  }
162}
163
164// This looks at solution and sets bounds to contain solution
165void 
166CbcLink::feasibleRegion()
167{
168  int j;
169  int firstNonZero=-1;
170  int lastNonZero = -1;
171  OsiSolverInterface * solver = model_->solver();
172  const double * solution = model_->testSolution();
173  const double * upper = solver->getColUpper();
174  double integerTolerance = 
175    model_->getDblParam(CbcModel::CbcIntegerTolerance);
176  double weight = 0.0;
177  double sum =0.0;
178
179  int base=first_;
180  for (j=0;j<numberMembers_;j++) {
181    for (int k=0;k<numberLinks_;k++) {
182      int iColumn = base+k;
183      double value = CoinMax(0.0,solution[iColumn]);
184      sum += value;
185      if (value>integerTolerance&&upper[iColumn]) {
186        weight += weights_[j]*value;
187        if (firstNonZero<0)
188          firstNonZero=j;
189        lastNonZero=j;
190      }
191    }
192    base += numberLinks_;
193  }
194  assert (lastNonZero-firstNonZero==0) ;
195  base=first_;
196  for (j=0;j<firstNonZero;j++) {
197    for (int k=0;k<numberLinks_;k++) {
198      int iColumn = base+k;
199      solver->setColUpper(iColumn,0.0);
200    }
201    base += numberLinks_;
202  }
203  // skip
204  base += numberLinks_;
205  for (j=lastNonZero+1;j<numberMembers_;j++) {
206    for (int k=0;k<numberLinks_;k++) {
207      int iColumn = base+k;
208      solver->setColUpper(iColumn,0.0);
209    }
210    base += numberLinks_;
211  }
212}
213
214
215// Creates a branching object
216CbcBranchingObject * 
217CbcLink::createBranch(int way) 
218{
219  int j;
220  const double * solution = model_->testSolution();
221  double integerTolerance = 
222      model_->getDblParam(CbcModel::CbcIntegerTolerance);
223  OsiSolverInterface * solver = model_->solver();
224  const double * upper = solver->getColUpper();
225  int firstNonFixed=-1;
226  int lastNonFixed=-1;
227  int firstNonZero=-1;
228  int lastNonZero = -1;
229  double weight = 0.0;
230  double sum =0.0;
231  int base=first_;
232  for (j=0;j<numberMembers_;j++) {
233    for (int k=0;k<numberLinks_;k++) {
234      int iColumn = base+k;
235      if (upper[iColumn]) {
236        double value = CoinMax(0.0,solution[iColumn]);
237        sum += value;
238        if (firstNonFixed<0)
239          firstNonFixed=j;
240        lastNonFixed=j;
241        if (value>integerTolerance) {
242          weight += weights_[j]*value;
243          if (firstNonZero<0)
244            firstNonZero=j;
245          lastNonZero=j;
246        }
247      }
248    }
249    base += numberLinks_;
250  }
251  assert (lastNonZero-firstNonZero>=1) ;
252  // find where to branch
253  assert (sum>0.0);
254  weight /= sum;
255  int iWhere;
256  double separator=0.0;
257  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
258    if (weight<weights_[iWhere+1])
259      break;
260  separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
261  // create object
262  CbcBranchingObject * branch;
263  branch = new CbcLinkBranchingObject(model_,this,way,separator);
264  return branch;
265}
266// Useful constructor
267CbcLinkBranchingObject::CbcLinkBranchingObject (CbcModel * model,
268                                              const CbcLink * set,
269                                              int way ,
270                                              double separator)
271  :CbcBranchingObject(model,set->id(),way,0.5)
272{
273  set_ = set;
274  separator_ = separator;
275}
276
277// Copy constructor
278CbcLinkBranchingObject::CbcLinkBranchingObject ( const CbcLinkBranchingObject & rhs) :CbcBranchingObject(rhs)
279{
280  set_=rhs.set_;
281  separator_ = rhs.separator_;
282}
283
284// Assignment operator
285CbcLinkBranchingObject & 
286CbcLinkBranchingObject::operator=( const CbcLinkBranchingObject& rhs)
287{
288  if (this != &rhs) {
289    CbcBranchingObject::operator=(rhs);
290    set_=rhs.set_;
291    separator_ = rhs.separator_;
292  }
293  return *this;
294}
295CbcBranchingObject * 
296CbcLinkBranchingObject::clone() const
297{ 
298  return (new CbcLinkBranchingObject(*this));
299}
300
301
302// Destructor
303CbcLinkBranchingObject::~CbcLinkBranchingObject ()
304{
305}
306double
307CbcLinkBranchingObject::branch(bool normalBranch)
308{
309  if (model_->messageHandler()->logLevel()>2&&normalBranch)
310    print(normalBranch);
311  numberBranchesLeft_--;
312  int numberMembers = set_->numberMembers();
313  int numberLinks = set_->numberLinks();
314  const double * weights = set_->weights();
315  OsiSolverInterface * solver = model_->solver();
316  //const double * lower = solver->getColLower();
317  //const double * upper = solver->getColUpper();
318  // *** for way - up means fix all those in down section
319  if (way_<0) {
320    int i;
321    for ( i=0;i<numberMembers;i++) {
322      if (weights[i] > separator_)
323        break;
324    }
325    assert (i<numberMembers);
326    int base=set_->first()+i*numberLinks;;
327    for (;i<numberMembers;i++) {
328      for (int k=0;k<numberLinks;k++) {
329        int iColumn = base+k;
330        solver->setColUpper(iColumn,0.0);
331      }
332      base += numberLinks;
333    }
334    way_=1;       // Swap direction
335  } else {
336    int i;
337    int base=set_->first();
338    for ( i=0;i<numberMembers;i++) { 
339      if (weights[i] >= separator_) {
340        break;
341      } else {
342        for (int k=0;k<numberLinks;k++) {
343          int iColumn = base+k;
344          solver->setColUpper(iColumn,0.0);
345        }
346        base += numberLinks;
347      }
348    }
349    assert (i<numberMembers);
350    way_=-1;      // Swap direction
351  }
352  return 0.0;
353}
354// Print what would happen 
355void
356CbcLinkBranchingObject::print(bool normalBranch)
357{
358  int numberMembers = set_->numberMembers();
359  int numberLinks = set_->numberLinks();
360  const double * weights = set_->weights();
361  OsiSolverInterface * solver = model_->solver();
362  const double * upper = solver->getColUpper();
363  int first=numberMembers;
364  int last=-1;
365  int numberFixed=0;
366  int numberOther=0;
367  int i;
368  int base=set_->first();
369  for ( i=0;i<numberMembers;i++) {
370    for (int k=0;k<numberLinks;k++) {
371      int iColumn = base+k;
372      double bound = upper[iColumn];
373      if (bound) {
374        first = CoinMin(first,i);
375        last = CoinMax(last,i);
376      }
377    }
378    base += numberLinks;
379  }
380  // *** for way - up means fix all those in down section
381  base=set_->first();
382  if (way_<0) {
383    printf("SOS Down");
384    for ( i=0;i<numberMembers;i++) {
385      for (int k=0;k<numberLinks;k++) {
386        int iColumn = base+k;
387        double bound = upper[iColumn];
388        if (weights[i] > separator_)
389          break;
390        else if (bound)
391          numberOther++;
392      }
393      base += numberLinks;
394    }
395    assert (i<numberMembers);
396    for (;i<numberMembers;i++) {
397      for (int k=0;k<numberLinks;k++) {
398        int iColumn = base+k;
399        double bound = upper[iColumn];
400        if (bound)
401          numberFixed++;
402      }
403      base += numberLinks;
404    }
405  } else {
406    printf("SOS Up");
407    for ( i=0;i<numberMembers;i++) {
408      for (int k=0;k<numberLinks;k++) {
409        int iColumn = base+k;
410        double bound = upper[iColumn];
411        if (weights[i] >= separator_)
412          break;
413        else if (bound)
414          numberFixed++;
415      }
416      base += numberLinks;
417    }
418    assert (i<numberMembers);
419    for (;i<numberMembers;i++) {
420      for (int k=0;k<numberLinks;k++) {
421        int iColumn = base+k;
422        double bound = upper[iColumn];
423        if (bound)
424          numberOther++;
425      }
426      base += numberLinks;
427    }
428  }
429  assert ((numberFixed%numberLinks)==0);
430  assert ((numberFixed%numberOther)==0);
431  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
432         separator_,first,weights[first],last,weights[last],numberFixed/numberLinks,
433         numberOther/numberLinks);
434}
Note: See TracBrowser for help on using the repository browser.