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

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

for sos

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.3 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(bool normalBranch)
567{
568  if (model_->messageHandler()->logLevel()>2&&normalBranch)
569    print(normalBranch);
570  numberBranchesLeft_--;
571  int numberMembers = set_->numberMembers();
572  int numberLinks = set_->numberLinks();
573  const double * weights = set_->weights();
574  const int * which = set_->which();
575  OsiSolverInterface * solver = model_->solver();
576  //const double * lower = solver->getColLower();
577  //const double * upper = solver->getColUpper();
578  // *** for way - up means fix all those in down section
579  if (way_<0) {
580    int i;
581    for ( i=0;i<numberMembers;i++) {
582      if (weights[i] > separator_)
583        break;
584    }
585    assert (i<numberMembers);
586    int base=i*numberLinks;;
587    for (;i<numberMembers;i++) {
588      for (int k=0;k<numberLinks;k++) {
589        int iColumn = which[base+k];
590        solver->setColUpper(iColumn,0.0);
591      }
592      base += numberLinks;
593    }
594    way_=1;       // Swap direction
595  } else {
596    int i;
597    int base=0;
598    for ( i=0;i<numberMembers;i++) { 
599      if (weights[i] >= separator_) {
600        break;
601      } else {
602        for (int k=0;k<numberLinks;k++) {
603          int iColumn = which[base+k];
604          solver->setColUpper(iColumn,0.0);
605        }
606        base += numberLinks;
607      }
608    }
609    assert (i<numberMembers);
610    way_=-1;      // Swap direction
611  }
612  return 0.0;
613}
614// Print what would happen 
615void
616CbcLinkBranchingObject::print(bool normalBranch)
617{
618  int numberMembers = set_->numberMembers();
619  int numberLinks = set_->numberLinks();
620  const double * weights = set_->weights();
621  const int * which = set_->which();
622  OsiSolverInterface * solver = model_->solver();
623  const double * upper = solver->getColUpper();
624  int first=numberMembers;
625  int last=-1;
626  int numberFixed=0;
627  int numberOther=0;
628  int i;
629  int base=0;
630  for ( i=0;i<numberMembers;i++) {
631    for (int k=0;k<numberLinks;k++) {
632      int iColumn = which[base+k];
633      double bound = upper[iColumn];
634      if (bound) {
635        first = CoinMin(first,i);
636        last = CoinMax(last,i);
637      }
638    }
639    base += numberLinks;
640  }
641  // *** for way - up means fix all those in down section
642  base=0;
643  if (way_<0) {
644    printf("SOS Down");
645    for ( i=0;i<numberMembers;i++) {
646      if (weights[i] > separator_) 
647        break;
648      for (int k=0;k<numberLinks;k++) {
649        int iColumn = which[base+k];
650        double bound = upper[iColumn];
651        if (bound)
652          numberOther++;
653      }
654      base += numberLinks;
655    }
656    assert (i<numberMembers);
657    for (;i<numberMembers;i++) {
658      for (int k=0;k<numberLinks;k++) {
659        int iColumn = which[base+k];
660        double bound = upper[iColumn];
661        if (bound)
662          numberFixed++;
663      }
664      base += numberLinks;
665    }
666  } else {
667    printf("SOS Up");
668    for ( i=0;i<numberMembers;i++) {
669      if (weights[i] >= separator_)
670        break;
671      for (int k=0;k<numberLinks;k++) {
672        int iColumn = which[base+k];
673        double bound = upper[iColumn];
674        if (bound)
675          numberFixed++;
676      }
677      base += numberLinks;
678    }
679    assert (i<numberMembers);
680    for (;i<numberMembers;i++) {
681      for (int k=0;k<numberLinks;k++) {
682        int iColumn = which[base+k];
683        double bound = upper[iColumn];
684        if (bound)
685          numberOther++;
686      }
687      base += numberLinks;
688    }
689  }
690  assert ((numberFixed%numberLinks)==0);
691  assert ((numberOther%numberLinks)==0);
692  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
693         separator_,first,weights[first],last,weights[last],numberFixed/numberLinks,
694         numberOther/numberLinks);
695}
Note: See TracBrowser for help on using the repository browser.