source: branches/devel/Cbc/examples/CbcBranchLink.cpp @ 424

Last change on this file since 424 was 424, checked in by forrest, 13 years ago

many changes

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