source: trunk/Cbc/examples/CbcBranchLink.cpp

Last change on this file was 1900, checked in by forrest, 5 years ago

fix some examples and a bug on repeated use of same CbcModel?

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