source: trunk/Cbc/examples/CbcBranchFollow2.cpp

Last change on this file was 2469, checked in by unxusr, 5 months ago

formatting

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