source: stable/2.8/Cbc/examples/CbcBranchFollow2.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: 13.0 KB
Line 
1// $Id: CbcBranchFollow2.cpp 1574 2011-01-05 01:13:55Z stefan $
2// Copyright (C) 2004, 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 "CbcBranchFollow2.hpp"
19#include "CoinSort.hpp"
20#include "CoinError.hpp"
21// Default Constructor
22CbcFollowOn2::CbcFollowOn2 ()
23  : CbcObject(),
24    rhs_(NULL),
25    maximumRhs_(1)
26{
27}
28
29// Useful constructor
30CbcFollowOn2::CbcFollowOn2 (CbcModel * model)
31  : CbcObject(model)
32{
33  assert (model);
34  OsiSolverInterface * solver = model_->solver();
35  matrix_ = *solver->getMatrixByCol();
36  matrix_.removeGaps();
37  matrixByRow_ = *solver->getMatrixByRow();
38  int numberRows = matrix_.getNumRows();
39  maximumRhs_ =1;
40 
41  rhs_ = new int[numberRows];
42  int i;
43  const double * rowLower = solver->getRowLower();
44  const double * rowUpper = solver->getRowUpper();
45  // Row copy
46  const double * elementByRow = matrixByRow_.getElements();
47  const int * column = matrixByRow_.getIndices();
48  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
49  const int * rowLength = matrixByRow_.getVectorLengths();
50  for (i=0;i<numberRows;i++) {
51    rhs_[i]=0;
52    double value = rowLower[i];
53    if (value==rowUpper[i]) {
54      if (floor(value)==value&&value>=1.0&&value<100.0) {
55        // check elements
56        bool good=true;
57        for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
58          int iColumn = column[j];
59          if (!solver->isInteger(iColumn))
60            good=false;
61          double elValue = elementByRow[j];
62          if (floor(elValue)!=elValue||elValue<1.0)
63            good=false;
64        }
65        if (good)
66          rhs_[i]=(int) value;
67      }
68    }
69  }
70}
71
72// Copy constructor
73CbcFollowOn2::CbcFollowOn2 ( const CbcFollowOn2 & rhs)
74  :CbcObject(rhs),
75   matrix_(rhs.matrix_),
76   matrixByRow_(rhs.matrixByRow_),
77   maximumRhs_(rhs.maximumRhs_)
78{
79  int numberRows = matrix_.getNumRows();
80  rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
81}
82
83// Clone
84CbcObject *
85CbcFollowOn2::clone() const
86{
87  return new CbcFollowOn2(*this);
88}
89
90// Assignment operator
91CbcFollowOn2 & 
92CbcFollowOn2::operator=( const CbcFollowOn2& rhs)
93{
94  if (this!=&rhs) {
95    CbcObject::operator=(rhs);
96    delete [] rhs_;
97    matrix_ = rhs.matrix_;
98    matrixByRow_ = rhs.matrixByRow_;
99    int numberRows = matrix_.getNumRows();
100    rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
101    maximumRhs_ = rhs.maximumRhs_;
102  }
103  return *this;
104}
105
106// Destructor
107CbcFollowOn2::~CbcFollowOn2 ()
108{
109  delete [] rhs_;
110}
111/* As some computation is needed in more than one place - returns row.
112   Also returns other row and effective rhs (so we can know if cut)
113*/
114int 
115CbcFollowOn2::gutsOfFollowOn2(int & otherRow, int & preferredWay,
116                              int & effectiveRhs) const
117{
118  int whichRow=-1;
119  otherRow=-1;
120  int numberRows = matrix_.getNumRows();
121 
122  int i;
123  // For sorting
124  int * sort = new int [numberRows];
125  int * isort = new int [numberRows];
126  // Column copy
127  //const double * element = matrix_.getElements();
128  const int * row = matrix_.getIndices();
129  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
130  const int * columnLength = matrix_.getVectorLengths();
131  // Row copy
132  const double * elementByRow = matrixByRow_.getElements();
133  const int * column = matrixByRow_.getIndices();
134  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
135  const int * rowLength = matrixByRow_.getVectorLengths();
136  OsiSolverInterface * solver = model_->solver();
137  const double * columnLower = solver->getColLower();
138  const double * columnUpper = solver->getColUpper();
139  const double * solution = solver->getColSolution();
140  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
141  int nSort=0;
142  for (i=0;i<numberRows;i++) {
143    if (rhs_[i]) {
144      // check elements
145      double smallest=1.0e10;
146      double largest=0.0;
147      int rhsValue=rhs_[i];
148      int number1=0;
149      int numberOther=0;
150      int numberUnsatisfied=0;
151      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
152        int iColumn = column[j];
153        double value = elementByRow[j];
154        double solValue = solution[iColumn];
155        if (columnLower[iColumn]!=columnUpper[iColumn]) {
156          smallest = CoinMin(smallest,value);
157          largest = CoinMax(largest,value);
158          if (value==1.0)
159            number1++;
160          else
161            numberOther++;
162          if (fabs(floor(solValue+0.5)-solValue)>integerTolerance)
163            numberUnsatisfied++;
164        } else {
165          rhsValue -= (int)(value*floor(solValue+0.5));
166        }
167      }
168      if (numberUnsatisfied>1) {
169        if (smallest<largest) {
170#if 0
171          if  (largest>rhsValue)
172            printf("could fix\n");
173          if (number1==1&&largest==rhsValue)
174            printf("could fix\n");
175#endif
176          if (rhsValue<=maximumRhs_&&0) {
177            // will mean a cut but worth trying
178            sort[nSort]=i;
179            isort[nSort++]=100000-numberUnsatisfied;
180          }
181        } else if (largest==rhsValue) {
182          sort[nSort]=i;
183          isort[nSort++]=-numberUnsatisfied;
184        }
185      }
186    }
187  }
188  if (nSort>1) {
189    CoinSort_2(isort,isort+nSort,sort);
190    assert (isort[1]<0);
191    CoinZeroN(isort,numberRows);
192    double * other = new double[numberRows];
193    CoinZeroN(other,numberRows);
194    int * which = new int[numberRows];
195    //#define COUNT
196#ifndef COUNT
197    bool beforeSolution = model_->getSolutionCount()==0;
198#endif
199    for (int k=0;k<nSort-1;k++) {
200      i=sort[k];
201      int numberUnsatisfied = 0;
202      int n=0;
203      int j;
204      for (j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
205        int iColumn = column[j];
206        if (columnLower[iColumn]!=columnUpper[iColumn]) {
207          double solValue = solution[iColumn]-columnLower[iColumn];
208          if (fabs(floor(solValue+0.5)-solValue)>integerTolerance) {
209            numberUnsatisfied++;
210            for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
211              int iRow = row[jj];
212              if (rhs_[iRow]) {
213                other[iRow]+=solValue;
214                if (isort[iRow]) {
215                  isort[iRow]++;
216                } else {
217                  isort[iRow]=1;
218                  which[n++]=iRow;
219                }
220              }
221            }
222          }
223        }
224      }
225      double total=0.0;
226      // Take out row
227      double sumThis=other[i];
228      other[i]=0.0;
229      assert (numberUnsatisfied==isort[i]);
230      // find one nearest half if solution, one if before solution
231      int iBest=-1;
232      double dtarget=0.5*total;
233#ifdef COUNT
234      int target = (numberUnsatisfied+1)>>1;
235      int best=numberUnsatisfied;
236#else
237      double best;
238      if (beforeSolution)
239        best=dtarget;
240      else
241        best=1.0e30;
242#endif
243      for (j=0;j<n;j++) {
244        int iRow = which[j];
245        double dvalue=other[iRow];
246        other[iRow]=0.0;
247#ifdef COUNT
248        int value = isort[iRow];
249#endif
250        isort[iRow]=0;
251        if (fabs(dvalue)<1.0e-8||fabs(sumThis-dvalue)<1.0e-8)
252          continue;
253        if (fabs(floor(dvalue+0.5)-dvalue)<integerTolerance)
254          continue;
255        dvalue -= floor(dvalue);
256#ifdef COUNT
257        if (abs(value-target)<best&&value!=numberUnsatisfied) {
258          best=abs(value-target);
259          iBest=iRow;
260          if (dvalue<dtarget)
261            preferredWay=1;
262          else
263            preferredWay=-1;
264        }
265#else
266        if (beforeSolution) {
267          if (fabs(dvalue-dtarget)>best) {
268            best = fabs(dvalue-dtarget);
269            iBest=iRow;
270            if (dvalue<dtarget)
271              preferredWay=1;
272            else
273              preferredWay=-1;
274          }
275        } else {
276          if (fabs(dvalue-dtarget)<best) {
277            best = fabs(dvalue-dtarget);
278            iBest=iRow;
279            if (dvalue<dtarget)
280              preferredWay=1;
281            else
282              preferredWay=-1;
283          }
284        }
285#endif
286      }
287      if (iBest>=0) {
288        whichRow=i;
289        otherRow=iBest;
290        //printf("Rows %d (%d) and %d (%d)\n",whichRow,rhs_[whichRow],
291        //     otherRow,rhs_[otherRow]);
292        break;
293      }
294    }
295    delete [] which;
296    delete [] other;
297  }
298  delete [] sort;
299  delete [] isort;
300  return whichRow;
301}
302
303// Infeasibility - large is 0.5
304double 
305CbcFollowOn2::infeasibility(int & preferredWay) const
306{
307  int otherRow=0;
308  int effectiveRhs;
309  int whichRow = gutsOfFollowOn2(otherRow,preferredWay,effectiveRhs);
310  if (whichRow<0) {
311    return 0.0;
312  } else {
313    assert (whichRow!=otherRow);
314    return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
315  }
316}
317
318// This looks at solution and sets bounds to contain solution
319void 
320CbcFollowOn2::feasibleRegion()
321{
322}
323
324
325// Creates a branching object
326CbcBranchingObject * 
327CbcFollowOn2::createBranch(int way) 
328{
329  int otherRow=0;
330  int preferredWay;
331  int effectiveRhs;
332  int whichRow = gutsOfFollowOn2(otherRow,preferredWay,effectiveRhs);
333  assert(way==preferredWay);
334  assert (whichRow>=0);
335  int numberColumns = matrix_.getNumCols();
336 
337  // Column copy
338  //const double * element = matrix_.getElements();
339  const int * row = matrix_.getIndices();
340  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
341  const int * columnLength = matrix_.getVectorLengths();
342  // Row copy
343  //const double * elementByRow = matrixByRow_.getElements();
344  const int * column = matrixByRow_.getIndices();
345  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
346  const int * rowLength = matrixByRow_.getVectorLengths();
347  OsiSolverInterface * solver = model_->solver();
348  const double * columnLower = solver->getColLower();
349  const double * columnUpper = solver->getColUpper();
350  //const double * solution = solver->getColSolution();
351#if 0
352  //printf("Rows %d (%d) and %d (%d)\n",whichRow,rhs_[whichRow],
353  //     otherRow,rhs_[otherRow]);
354  int nFree=0;
355  int nOut=0;
356  int nImplicit=0;
357  int i;
358  int rhsx[100];
359  double * colUpper2 = new double [numberColumns];
360  memcpy(rhsx,rhs_,matrix_.getNumRows()*sizeof(int));
361  for ( i=0;i<numberColumns;i++) {
362    colUpper2[i]=columnUpper[i];
363    if (columnLower[i]==columnUpper[i]) {
364      for (int jj=columnStart[i];jj<columnStart[i]+columnLength[i];jj++) {
365        int iRow = row[jj];
366        nOut += (int) floor(element[jj]*solution[i]+0.5);
367        rhsx[iRow] -= (int) floor(element[jj]*solution[i]+0.5);
368      }
369    }
370  }
371  int nFixedBut=0;
372  for ( i=0;i<numberColumns;i++) {
373    if (columnLower[i]!=columnUpper[i]) {
374      nFree++;
375      bool nonzero=false;
376      if (fabs(solution[i])>1.0e-5) {
377        nonzero=true;
378        //printf("column %d value %g ",i,solution[i]);
379        for (int jj=columnStart[i];jj<columnStart[i]+columnLength[i];jj++) {
380          //int iRow = row[jj];
381          //printf("(%d,%g) ",iRow,element[jj]);
382        }
383        //printf("\n");
384      }
385      bool fixed=false;
386      for (int jj=columnStart[i];jj<columnStart[i]+columnLength[i];jj++) {
387        int iRow = row[jj];
388        if (element[jj]>rhsx[iRow])
389          fixed=true;
390      }
391      if (fixed) {
392        nImplicit++;
393        colUpper2[i]=0.0;
394        if (nonzero)
395          nFixedBut++;
396        assert (!columnLower[i]);
397      }
398    }
399  }
400  // See if anything odd
401  char * check = new char[numberColumns];
402  memset(check,0,numberColumns);
403  int * which2 = new int[numberColumns];
404  int numberRows=matrix_.getNumRows();
405  for (i=0;i<numberRows;i++) {
406    if (rhsx[i]==1) {
407      int nn=0;
408      int j,k;
409      for (j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
410        int iColumn = column[j];
411        double value = elementByRow[j];
412        if (columnLower[iColumn]!=colUpper2[iColumn]) {
413          assert (value==1.0);
414          check[iColumn]=1;
415          which2[nn++]=iColumn;
416        }
417      }
418      for ( k=i+1;k<numberRows;k++) {
419        if (rhsx[k]==1) {
420          int nn2=0;
421          int nnsame=0;
422          for (int j=rowStart[k];j<rowStart[k]+rowLength[k];j++) {
423            int iColumn = column[j];
424            double value = elementByRow[j];
425            if (columnLower[iColumn]!=colUpper2[iColumn]) {
426              assert (value==1.0);
427              nn2++;
428              if (check[iColumn])
429                nnsame++;
430            }
431          }
432          if (nnsame==nn2) {
433            if (nn2<nn)
434              printf("row %d strict subset of row %d, fix some in row %d\n",
435                     k,i,i);
436            else if (nn2==nn)
437              printf("row %d identical to row %d\n",
438                     k,i);
439            else if (nn2>=nn)
440              abort();
441          } else if (nnsame==nn&&nn2>nn) {
442              printf("row %d strict superset of row %d, fix some in row %d\n",
443                     k,i,k);
444          }
445        }
446      }
447      for (k=0;k<nn;k++)
448        check[which2[k]]=0;
449     
450    }
451  }
452  delete [] check;
453  delete [] which2;
454  delete [] colUpper2;
455  printf("%d free (but %d implicitly fixed of which %d nonzero), %d out of rhs\n",nFree,nImplicit,nFixedBut,nOut);
456#endif
457  int nUp=0;
458  int nDown=0;
459  int * upList = new int[numberColumns];
460  int * downList = new int[numberColumns];
461  int j;
462  for (j=rowStart[whichRow];j<rowStart[whichRow]+rowLength[whichRow];j++) {
463    int iColumn = column[j];
464    if (columnLower[iColumn]!=columnUpper[iColumn]) {
465      bool up=true;
466      for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
467        int iRow = row[jj];
468        if (iRow==otherRow) {
469          up=false;
470          break;
471        }
472      }
473      if (up)
474        upList[nUp++]=iColumn;
475      else
476        downList[nDown++]=iColumn;
477    }
478  }
479  //printf("way %d\n",way);
480  // create object
481  CbcBranchingObject * branch
482     = new CbcFixingBranchingObject(model_,way,
483                                         nDown,downList,nUp,upList);
484  delete [] upList;
485  delete [] downList;
486  return branch;
487}
Note: See TracBrowser for help on using the repository browser.