source: trunk/Cbc/examples/CbcBranchFollow2.cpp

Last change on this file was 1898, checked in by stefan, 5 years ago

fixup examples

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