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

Last change on this file since 1173 was 1173, checked in by forrest, 10 years ago

add $id

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