source: trunk/Cbc/examples/CbcBranchLink.cpp @ 1565

Last change on this file since 1565 was 932, checked in by forrest, 11 years ago

allow Compare

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.7 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}
694/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
695    same type and must have the same original object, but they may have
696    different feasible regions.
697    Return the appropriate CbcRangeCompare value (first argument being the
698    sub/superset if that's the case). In case of overlap (and if \c
699    replaceIfOverlap is true) replace the current branching object with one
700    whose feasible region is the overlap.
701*/
702CbcRangeCompare
703CbcLinkBranchingObject::compareBranchingObject
704(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
705{
706  throw("must implement");
707}
Note: See TracBrowser for help on using the repository browser.