source: stable/2.8/Cbc/examples/CbcBranchLink.cpp @ 1856

Last change on this file since 1856 was 1574, checked in by lou, 8 years ago

Change to EPL license notice.

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