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

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

try and keep examples in sync

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.2 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#include "CoinPackedMatrix.hpp"
18
19// Default Constructor
20CbcLink::CbcLink ()
21  : CbcObject(),
22    weights_(NULL),
23    numberMembers_(0),
24    numberLinks_(0),
25    which_(NULL),
26    sosType_(1)
27{
28}
29
30// Useful constructor (which are indices)
31CbcLink::CbcLink (CbcModel * model,  int numberMembers,
32           int numberLinks, int first , const double * weights, int identifier)
33  : CbcObject(model),
34    numberMembers_(numberMembers),
35    numberLinks_(numberLinks),
36    which_(NULL),
37    sosType_(1)
38{
39  id_=identifier;
40  if (numberMembers_) {
41    weights_ = new double[numberMembers_];
42    which_ = new int[numberMembers_*numberLinks_];
43    if (weights) {
44      memcpy(weights_,weights,numberMembers_*sizeof(double));
45    } else {
46      for (int i=0;i<numberMembers_;i++)
47        weights_[i]=i;
48    }
49    // weights must be increasing
50    int i;
51    double last=-COIN_DBL_MAX;
52    for (i=0;i<numberMembers_;i++) {
53      assert (weights_[i]>last+1.0e-12);
54      last=weights_[i];
55    }
56    for (i=0;i<numberMembers_*numberLinks_;i++) {
57      which_[i]=first+i;
58    }
59  } else {
60    weights_ = NULL;
61  }
62}
63
64// Useful constructor (which are indices)
65CbcLink::CbcLink (CbcModel * model,  int numberMembers,
66           int numberLinks, int sosType, const int * which , const double * weights, int identifier)
67  : CbcObject(model),
68    numberMembers_(numberMembers),
69    numberLinks_(numberLinks),
70    which_(NULL),
71    sosType_(sosType)
72{
73  id_=identifier;
74  if (numberMembers_) {
75    weights_ = new double[numberMembers_];
76    which_ = new int[numberMembers_*numberLinks_];
77    if (weights) {
78      memcpy(weights_,weights,numberMembers_*sizeof(double));
79    } else {
80      for (int i=0;i<numberMembers_;i++)
81        weights_[i]=i;
82    }
83    // weights must be increasing
84    int i;
85    double last=-COIN_DBL_MAX;
86    for (i=0;i<numberMembers_;i++) {
87      assert (weights_[i]>last+1.0e-12);
88      last=weights_[i];
89    }
90    for (i=0;i<numberMembers_*numberLinks_;i++) {
91      which_[i]= which[i];
92    }
93  } else {
94    weights_ = NULL;
95  }
96}
97
98// Copy constructor
99CbcLink::CbcLink ( const CbcLink & rhs)
100  :CbcObject(rhs)
101{
102  numberMembers_ = rhs.numberMembers_;
103  numberLinks_ = rhs.numberLinks_;
104  sosType_ = rhs.sosType_;
105  if (numberMembers_) {
106    weights_ = CoinCopyOfArray(rhs.weights_,numberMembers_);
107    which_ = CoinCopyOfArray(rhs.which_,numberMembers_*numberLinks_);
108  } else {
109    weights_ = NULL;
110    which_ = NULL;
111  }
112}
113
114// Clone
115CbcObject *
116CbcLink::clone() const
117{
118  return new CbcLink(*this);
119}
120
121// Assignment operator
122CbcLink & 
123CbcLink::operator=( const CbcLink& rhs)
124{
125  if (this!=&rhs) {
126    CbcObject::operator=(rhs);
127    delete [] weights_;
128    delete [] which_;
129    numberMembers_ = rhs.numberMembers_;
130    numberLinks_ = rhs.numberLinks_;
131    sosType_ = rhs.sosType_;
132    if (numberMembers_) {
133      weights_ = CoinCopyOfArray(rhs.weights_,numberMembers_);
134      which_ = CoinCopyOfArray(rhs.which_,numberMembers_*numberLinks_);
135    } else {
136      weights_ = NULL;
137      which_ = NULL;
138    }
139  }
140  return *this;
141}
142
143// Destructor
144CbcLink::~CbcLink ()
145{
146  delete [] weights_;
147  delete [] which_;
148}
149
150// Infeasibility - large is 0.5
151double 
152CbcLink::infeasibility(int & preferredWay) const
153{
154  int j;
155  int firstNonZero=-1;
156  int lastNonZero = -1;
157  OsiSolverInterface * solver = model_->solver();
158  const double * solution = model_->testSolution();
159  //const double * lower = solver->getColLower();
160  const double * upper = solver->getColUpper();
161  double integerTolerance = 
162    model_->getDblParam(CbcModel::CbcIntegerTolerance);
163  double weight = 0.0;
164  double sum =0.0;
165
166  // check bounds etc
167  double lastWeight=-1.0e100;
168  int base=0;
169  for (j=0;j<numberMembers_;j++) {
170    for (int k=0;k<numberLinks_;k++) {
171      int iColumn = which_[base+k];
172      //if (lower[iColumn])
173      //throw CoinError("Non zero lower bound in CBCLink","infeasibility","CbcLink");
174      if (lastWeight>=weights_[j]-1.0e-7)
175        throw CoinError("Weights too close together in CBCLink","infeasibility","CbcLink");
176      double value = CoinMax(0.0,solution[iColumn]);
177      sum += value;
178      if (value>integerTolerance&&upper[iColumn]) {
179        // Possibly due to scaling a fixed variable might slip through
180        if (value>upper[iColumn]+1.0e-8) {
181          // Could change to #ifdef CBC_DEBUG
182#ifndef NDEBUG
183          if (model_->messageHandler()->logLevel()>1)
184            printf("** Variable %d (%d) has value %g and upper bound of %g\n",
185                   iColumn,j,value,upper[iColumn]);
186#endif
187        } 
188        value = CoinMin(value,upper[iColumn]);
189        weight += weights_[j]*value;
190        if (firstNonZero<0)
191          firstNonZero=j;
192        lastNonZero=j;
193      }
194    }
195    base += numberLinks_;
196  }
197  double valueInfeasibility;
198  preferredWay=1;
199  if (lastNonZero-firstNonZero>=sosType_) {
200    // find where to branch
201    assert (sum>0.0);
202    weight /= sum;
203    valueInfeasibility = lastNonZero-firstNonZero+1;
204    valueInfeasibility *= 0.5/((double) numberMembers_);
205    //#define DISTANCE
206#ifdef DISTANCE
207    assert (sosType_==1); // code up
208    /* may still be satisfied.
209       For LOS type 2 we might wish to move coding around
210       and keep initial info in model_ for speed
211    */
212    int iWhere;
213    bool possible=false;
214    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
215      if (fabs(weight-weights_[iWhere])<1.0e-8) {
216        possible=true;
217        break;
218      }
219    }
220    if (possible) {
221      // One could move some of this (+ arrays) into model_
222      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
223      const double * element = matrix->getMutableElements();
224      const int * row = matrix->getIndices();
225      const CoinBigIndex * columnStart = matrix->getVectorStarts();
226      const int * columnLength = matrix->getVectorLengths();
227      const double * rowSolution = solver->getRowActivity();
228      const double * rowLower = solver->getRowLower();
229      const double * rowUpper = solver->getRowUpper();
230      int numberRows = matrix->getNumRows();
231      double * array = new double [numberRows];
232      CoinZeroN(array,numberRows);
233      int * which = new int [numberRows];
234      int n=0;
235      int base=numberLinks_*firstNonZero;
236      for (j=firstNonZero;j<=lastNonZero;j++) {
237        for (int k=0;k<numberLinks_;k++) {
238          int iColumn = which_[base+k];
239          double value = CoinMax(0.0,solution[iColumn]);
240          if (value>integerTolerance&&upper[iColumn]) {
241            value = CoinMin(value,upper[iColumn]);
242            for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
243              int iRow = row[j];
244              double a = array[iRow];
245              if (a) {
246                a += value*element[j];
247                if (!a)
248                  a = 1.0e-100;
249              } else {
250                which[n++]=iRow;
251                a=value*element[j];
252                assert (a);
253              }
254              array[iRow]=a;
255            }
256          }
257        }
258        base += numberLinks_;
259      }
260      base=numberLinks_*iWhere;
261      for (int k=0;k<numberLinks_;k++) {
262        int iColumn = which_[base+k];
263        const double value = 1.0;
264        for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
265          int iRow = row[j];
266          double a = array[iRow];
267          if (a) {
268            a -= value*element[j];
269            if (!a)
270              a = 1.0e-100;
271          } else {
272            which[n++]=iRow;
273            a=-value*element[j];
274            assert (a);
275          }
276          array[iRow]=a;
277        }
278      }
279      for (j=0;j<n;j++) {
280        int iRow = which[j];
281        // moving to point will increase row solution by this
282        double distance = array[iRow];
283        if (distance>1.0e-8) {
284          if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
285            possible=false;
286            break;
287          }
288        } else if (distance<-1.0e-8) {
289          if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
290            possible=false;
291            break;
292          } 
293        }
294      }
295      for (j=0;j<n;j++)
296        array[which[j]]=0.0;
297      delete [] array;
298      delete [] which;
299      if (possible) {
300        valueInfeasibility=0.0;
301        printf("possible %d %d %d\n",firstNonZero,lastNonZero,iWhere);
302      }
303    }
304#endif
305  } else {
306    valueInfeasibility = 0.0; // satisfied
307  }
308  return valueInfeasibility;
309}
310
311// This looks at solution and sets bounds to contain solution
312void 
313CbcLink::feasibleRegion()
314{
315  int j;
316  int firstNonZero=-1;
317  int lastNonZero = -1;
318  OsiSolverInterface * solver = model_->solver();
319  const double * solution = model_->testSolution();
320  const double * upper = solver->getColUpper();
321  double integerTolerance = 
322    model_->getDblParam(CbcModel::CbcIntegerTolerance);
323  double weight = 0.0;
324  double sum =0.0;
325
326  int base=0;
327  for (j=0;j<numberMembers_;j++) {
328    for (int k=0;k<numberLinks_;k++) {
329      int iColumn = which_[base+k];
330      double value = CoinMax(0.0,solution[iColumn]);
331      sum += value;
332      if (value>integerTolerance&&upper[iColumn]) {
333        weight += weights_[j]*value;
334        if (firstNonZero<0)
335          firstNonZero=j;
336        lastNonZero=j;
337      }
338    }
339    base += numberLinks_;
340  }
341#ifdef DISTANCE
342  if (lastNonZero-firstNonZero>sosType_-1) {
343    /* may still be satisfied.
344       For LOS type 2 we might wish to move coding around
345       and keep initial info in model_ for speed
346    */
347    int iWhere;
348    bool possible=false;
349    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
350      if (fabs(weight-weights_[iWhere])<1.0e-8) {
351        possible=true;
352        break;
353      }
354    }
355    if (possible) {
356      // One could move some of this (+ arrays) into model_
357      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
358      const double * element = matrix->getMutableElements();
359      const int * row = matrix->getIndices();
360      const CoinBigIndex * columnStart = matrix->getVectorStarts();
361      const int * columnLength = matrix->getVectorLengths();
362      const double * rowSolution = solver->getRowActivity();
363      const double * rowLower = solver->getRowLower();
364      const double * rowUpper = solver->getRowUpper();
365      int numberRows = matrix->getNumRows();
366      double * array = new double [numberRows];
367      CoinZeroN(array,numberRows);
368      int * which = new int [numberRows];
369      int n=0;
370      int base=numberLinks_*firstNonZero;
371      for (j=firstNonZero;j<=lastNonZero;j++) {
372        for (int k=0;k<numberLinks_;k++) {
373          int iColumn = which_[base+k];
374          double value = CoinMax(0.0,solution[iColumn]);
375          if (value>integerTolerance&&upper[iColumn]) {
376            value = CoinMin(value,upper[iColumn]);
377            for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
378              int iRow = row[j];
379              double a = array[iRow];
380              if (a) {
381                a += value*element[j];
382                if (!a)
383                  a = 1.0e-100;
384              } else {
385                which[n++]=iRow;
386                a=value*element[j];
387                assert (a);
388              }
389              array[iRow]=a;
390            }
391          }
392        }
393        base += numberLinks_;
394      }
395      base=numberLinks_*iWhere;
396      for (int k=0;k<numberLinks_;k++) {
397        int iColumn = which_[base+k];
398        const double value = 1.0;
399        for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
400          int iRow = row[j];
401          double a = array[iRow];
402          if (a) {
403            a -= value*element[j];
404            if (!a)
405              a = 1.0e-100;
406          } else {
407            which[n++]=iRow;
408            a=-value*element[j];
409            assert (a);
410          }
411          array[iRow]=a;
412        }
413      }
414      for (j=0;j<n;j++) {
415        int iRow = which[j];
416        // moving to point will increase row solution by this
417        double distance = array[iRow];
418        if (distance>1.0e-8) {
419          if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
420            possible=false;
421            break;
422          }
423        } else if (distance<-1.0e-8) {
424          if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
425            possible=false;
426            break;
427          } 
428        }
429      }
430      for (j=0;j<n;j++)
431        array[which[j]]=0.0;
432      delete [] array;
433      delete [] which;
434      if (possible) {
435        printf("possible feas region %d %d %d\n",firstNonZero,lastNonZero,iWhere);
436        firstNonZero=iWhere;
437        lastNonZero=iWhere;
438      }
439    }
440  }
441#else
442  assert (lastNonZero-firstNonZero<sosType_) ;
443#endif
444  base=0;
445  for (j=0;j<firstNonZero;j++) {
446    for (int k=0;k<numberLinks_;k++) {
447      int iColumn = which_[base+k];
448      solver->setColUpper(iColumn,0.0);
449    }
450    base += numberLinks_;
451  }
452  // skip
453  base += numberLinks_;
454  for (j=lastNonZero+1;j<numberMembers_;j++) {
455    for (int k=0;k<numberLinks_;k++) {
456      int iColumn = which_[base+k];
457      solver->setColUpper(iColumn,0.0);
458    }
459    base += numberLinks_;
460  }
461}
462
463
464// Creates a branching object
465CbcBranchingObject * 
466CbcLink::createBranch(int way) 
467{
468  int j;
469  const double * solution = model_->testSolution();
470  double integerTolerance = 
471      model_->getDblParam(CbcModel::CbcIntegerTolerance);
472  OsiSolverInterface * solver = model_->solver();
473  const double * upper = solver->getColUpper();
474  int firstNonFixed=-1;
475  int lastNonFixed=-1;
476  int firstNonZero=-1;
477  int lastNonZero = -1;
478  double weight = 0.0;
479  double sum =0.0;
480  int base=0;
481  for (j=0;j<numberMembers_;j++) {
482    for (int k=0;k<numberLinks_;k++) {
483      int iColumn = which_[base+k];
484      if (upper[iColumn]) {
485        double value = CoinMax(0.0,solution[iColumn]);
486        sum += value;
487        if (firstNonFixed<0)
488          firstNonFixed=j;
489        lastNonFixed=j;
490        if (value>integerTolerance) {
491          weight += weights_[j]*value;
492          if (firstNonZero<0)
493            firstNonZero=j;
494          lastNonZero=j;
495        }
496      }
497    }
498    base += numberLinks_;
499  }
500  assert (lastNonZero-firstNonZero>=sosType_) ;
501  // find where to branch
502  assert (sum>0.0);
503  weight /= sum;
504  int iWhere;
505  double separator=0.0;
506  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
507    if (weight<weights_[iWhere+1])
508      break;
509  if (sosType_==1) {
510    // SOS 1
511    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
512  } else {
513    // SOS 2
514    if (iWhere==firstNonFixed)
515      iWhere++;;
516    if (iWhere==lastNonFixed-1)
517      iWhere = lastNonFixed-2;
518    separator = weights_[iWhere+1];
519  }
520  // create object
521  CbcBranchingObject * branch;
522  branch = new CbcLinkBranchingObject(model_,this,way,separator);
523  return branch;
524}
525// Useful constructor
526CbcLinkBranchingObject::CbcLinkBranchingObject (CbcModel * model,
527                                              const CbcLink * set,
528                                              int way ,
529                                              double separator)
530  :CbcBranchingObject(model,set->id(),way,0.5)
531{
532  set_ = set;
533  separator_ = separator;
534}
535
536// Copy constructor
537CbcLinkBranchingObject::CbcLinkBranchingObject ( const CbcLinkBranchingObject & rhs) :CbcBranchingObject(rhs)
538{
539  set_=rhs.set_;
540  separator_ = rhs.separator_;
541}
542
543// Assignment operator
544CbcLinkBranchingObject & 
545CbcLinkBranchingObject::operator=( const CbcLinkBranchingObject& rhs)
546{
547  if (this != &rhs) {
548    CbcBranchingObject::operator=(rhs);
549    set_=rhs.set_;
550    separator_ = rhs.separator_;
551  }
552  return *this;
553}
554CbcBranchingObject * 
555CbcLinkBranchingObject::clone() const
556{ 
557  return (new CbcLinkBranchingObject(*this));
558}
559
560
561// Destructor
562CbcLinkBranchingObject::~CbcLinkBranchingObject ()
563{
564}
565double
566CbcLinkBranchingObject::branch()
567{
568  decrementNumberBranchesLeft();
569  int numberMembers = set_->numberMembers();
570  int numberLinks = set_->numberLinks();
571  const double * weights = set_->weights();
572  const int * which = set_->which();
573  OsiSolverInterface * solver = model_->solver();
574  //const double * lower = solver->getColLower();
575  //const double * upper = solver->getColUpper();
576  // *** for way - up means fix all those in down section
577  if (way_<0) {
578    int i;
579    for ( i=0;i<numberMembers;i++) {
580      if (weights[i] > separator_)
581        break;
582    }
583    assert (i<numberMembers);
584    int base=i*numberLinks;;
585    for (;i<numberMembers;i++) {
586      for (int k=0;k<numberLinks;k++) {
587        int iColumn = which[base+k];
588        solver->setColUpper(iColumn,0.0);
589      }
590      base += numberLinks;
591    }
592    way_=1;       // Swap direction
593  } else {
594    int i;
595    int base=0;
596    for ( i=0;i<numberMembers;i++) { 
597      if (weights[i] >= separator_) {
598        break;
599      } else {
600        for (int k=0;k<numberLinks;k++) {
601          int iColumn = which[base+k];
602          solver->setColUpper(iColumn,0.0);
603        }
604        base += numberLinks;
605      }
606    }
607    assert (i<numberMembers);
608    way_=-1;      // Swap direction
609  }
610  return 0.0;
611}
612// Print what would happen 
613void
614CbcLinkBranchingObject::print()
615{
616  int numberMembers = set_->numberMembers();
617  int numberLinks = set_->numberLinks();
618  const double * weights = set_->weights();
619  const int * which = set_->which();
620  OsiSolverInterface * solver = model_->solver();
621  const double * upper = solver->getColUpper();
622  int first=numberMembers;
623  int last=-1;
624  int numberFixed=0;
625  int numberOther=0;
626  int i;
627  int base=0;
628  for ( i=0;i<numberMembers;i++) {
629    for (int k=0;k<numberLinks;k++) {
630      int iColumn = which[base+k];
631      double bound = upper[iColumn];
632      if (bound) {
633        first = CoinMin(first,i);
634        last = CoinMax(last,i);
635      }
636    }
637    base += numberLinks;
638  }
639  // *** for way - up means fix all those in down section
640  base=0;
641  if (way_<0) {
642    printf("SOS Down");
643    for ( i=0;i<numberMembers;i++) {
644      if (weights[i] > separator_) 
645        break;
646      for (int k=0;k<numberLinks;k++) {
647        int iColumn = which[base+k];
648        double bound = upper[iColumn];
649        if (bound)
650          numberOther++;
651      }
652      base += numberLinks;
653    }
654    assert (i<numberMembers);
655    for (;i<numberMembers;i++) {
656      for (int k=0;k<numberLinks;k++) {
657        int iColumn = which[base+k];
658        double bound = upper[iColumn];
659        if (bound)
660          numberFixed++;
661      }
662      base += numberLinks;
663    }
664  } else {
665    printf("SOS Up");
666    for ( i=0;i<numberMembers;i++) {
667      if (weights[i] >= separator_)
668        break;
669      for (int k=0;k<numberLinks;k++) {
670        int iColumn = which[base+k];
671        double bound = upper[iColumn];
672        if (bound)
673          numberFixed++;
674      }
675      base += numberLinks;
676    }
677    assert (i<numberMembers);
678    for (;i<numberMembers;i++) {
679      for (int k=0;k<numberLinks;k++) {
680        int iColumn = which[base+k];
681        double bound = upper[iColumn];
682        if (bound)
683          numberOther++;
684      }
685      base += numberLinks;
686    }
687  }
688  assert ((numberFixed%numberLinks)==0);
689  assert ((numberOther%numberLinks)==0);
690  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
691         separator_,first,weights[first],last,weights[last],numberFixed/numberLinks,
692         numberOther/numberLinks);
693}
Note: See TracBrowser for help on using the repository browser.