source: trunk/Cbc/src/CbcBranchToFixLots.cpp

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

spaces after angles

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.6 KB
Line 
1// $Id: CbcBranchToFixLots.cpp 2467 2019-01-03 21:26:29Z 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// Edwin 11/13/2009-- carved out of CbcBranchCut
7
8#if defined(_MSC_VER)
9// Turn off compiler warning about long names
10#pragma warning(disable : 4786)
11#endif
12#include <cassert>
13#include <cstdlib>
14#include <cmath>
15#include <cfloat>
16//#define CBC_DEBUG
17
18#include "OsiSolverInterface.hpp"
19#include "CbcModel.hpp"
20#include "CbcMessage.hpp"
21#include "CbcBranchCut.hpp"
22#include "CoinSort.hpp"
23#include "CoinError.hpp"
24#include "CbcBranchToFixLots.hpp"
25
26/** Default Constructor
27
28  Equivalent to an unspecified binary variable.
29*/
30CbcBranchToFixLots::CbcBranchToFixLots()
31  : CbcBranchCut()
32  , djTolerance_(COIN_DBL_MAX)
33  , fractionFixed_(1.0)
34  , mark_(NULL)
35  , depth_(-1)
36  , numberClean_(0)
37  , alwaysCreate_(false)
38{
39}
40
41/* Useful constructor - passed reduced cost tolerance and fraction we would like fixed.
42   Also depth level to do at.
43   Also passed number of 1 rows which when clean triggers fix
44   Always does if all 1 rows cleaned up and number>0 or if fraction columns reached
45   Also whether to create branch if can't reach fraction.
46*/
47CbcBranchToFixLots::CbcBranchToFixLots(CbcModel *model, double djTolerance,
48  double fractionFixed, int depth,
49  int numberClean,
50  const char *mark, bool alwaysCreate)
51  : CbcBranchCut(model)
52{
53  djTolerance_ = djTolerance;
54  fractionFixed_ = fractionFixed;
55  if (mark) {
56    int numberColumns = model->getNumCols();
57    mark_ = new char[numberColumns];
58    memcpy(mark_, mark, numberColumns);
59  } else {
60    mark_ = NULL;
61  }
62  depth_ = depth;
63  assert(model);
64  OsiSolverInterface *solver = model_->solver();
65  matrixByRow_ = *solver->getMatrixByRow();
66  numberClean_ = numberClean;
67  alwaysCreate_ = alwaysCreate;
68}
69// Copy constructor
70CbcBranchToFixLots::CbcBranchToFixLots(const CbcBranchToFixLots &rhs)
71  : CbcBranchCut(rhs)
72{
73  djTolerance_ = rhs.djTolerance_;
74  fractionFixed_ = rhs.fractionFixed_;
75  int numberColumns = model_->getNumCols();
76  mark_ = CoinCopyOfArray(rhs.mark_, numberColumns);
77  matrixByRow_ = rhs.matrixByRow_;
78  depth_ = rhs.depth_;
79  numberClean_ = rhs.numberClean_;
80  alwaysCreate_ = rhs.alwaysCreate_;
81}
82
83// Clone
84CbcObject *
85CbcBranchToFixLots::clone() const
86{
87  return new CbcBranchToFixLots(*this);
88}
89
90// Assignment operator
91CbcBranchToFixLots &
92CbcBranchToFixLots::operator=(const CbcBranchToFixLots &rhs)
93{
94  if (this != &rhs) {
95    CbcBranchCut::operator=(rhs);
96    djTolerance_ = rhs.djTolerance_;
97    fractionFixed_ = rhs.fractionFixed_;
98    int numberColumns = model_->getNumCols();
99    delete[] mark_;
100    mark_ = CoinCopyOfArray(rhs.mark_, numberColumns);
101    matrixByRow_ = rhs.matrixByRow_;
102    depth_ = rhs.depth_;
103    numberClean_ = rhs.numberClean_;
104    alwaysCreate_ = rhs.alwaysCreate_;
105  }
106  return *this;
107}
108
109// Destructor
110CbcBranchToFixLots::~CbcBranchToFixLots()
111{
112  delete[] mark_;
113}
114CbcBranchingObject *
115CbcBranchToFixLots::createCbcBranch(OsiSolverInterface *solver, const OsiBranchingInformation * /*info*/, int /*way*/)
116{
117  // by default way must be -1
118  //assert (way==-1);
119  //OsiSolverInterface * solver = model_->solver();
120  const double *solution = model_->testSolution();
121  const double *lower = solver->getColLower();
122  const double *upper = solver->getColUpper();
123  const double *dj = solver->getReducedCost();
124  int i;
125  int numberIntegers = model_->numberIntegers();
126  const int *integerVariable = model_->integerVariable();
127  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
128  // make smaller ?
129  double tolerance = CoinMin(1.0e-8, integerTolerance);
130  // How many fixed are we aiming at
131  int wantedFixed = static_cast< int >(static_cast< double >(numberIntegers) * fractionFixed_);
132  int nSort = 0;
133  int numberFixed = 0;
134  int numberColumns = solver->getNumCols();
135  int *sort = new int[numberColumns];
136  double *dsort = new double[numberColumns];
137  if (djTolerance_ != -1.234567) {
138    int type = shallWe();
139    assert(type);
140    // Take clean first
141    if (type == 1) {
142      for (i = 0; i < numberIntegers; i++) {
143        int iColumn = integerVariable[i];
144        if (upper[iColumn] > lower[iColumn]) {
145          if (!mark_ || !mark_[iColumn]) {
146            if (solution[iColumn] < lower[iColumn] + tolerance) {
147              if (dj[iColumn] > djTolerance_) {
148                dsort[nSort] = -dj[iColumn];
149                sort[nSort++] = iColumn;
150              }
151            } else if (solution[iColumn] > upper[iColumn] - tolerance) {
152              if (dj[iColumn] < -djTolerance_) {
153                dsort[nSort] = dj[iColumn];
154                sort[nSort++] = iColumn;
155              }
156            }
157          }
158        } else {
159          numberFixed++;
160        }
161      }
162      // sort
163      CoinSort_2(dsort, dsort + nSort, sort);
164      nSort = CoinMin(nSort, wantedFixed - numberFixed);
165    } else if (type < 10) {
166      int i;
167      //const double * rowLower = solver->getRowLower();
168      const double *rowUpper = solver->getRowUpper();
169      // Row copy
170      const double *elementByRow = matrixByRow_.getElements();
171      const int *column = matrixByRow_.getIndices();
172      const CoinBigIndex *rowStart = matrixByRow_.getVectorStarts();
173      const int *rowLength = matrixByRow_.getVectorLengths();
174      const double *columnLower = solver->getColLower();
175      const double *columnUpper = solver->getColUpper();
176      const double *solution = solver->getColSolution();
177      int numberColumns = solver->getNumCols();
178      int numberRows = solver->getNumRows();
179      for (i = 0; i < numberColumns; i++) {
180        sort[i] = i;
181        if (columnLower[i] != columnUpper[i]) {
182          dsort[i] = 1.0e100;
183        } else {
184          dsort[i] = 1.0e50;
185          numberFixed++;
186        }
187      }
188      for (i = 0; i < numberRows; i++) {
189        double rhsValue = rowUpper[i];
190        bool oneRow = true;
191        // check elements
192        int numberUnsatisfied = 0;
193        for (CoinBigIndex j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
194          int iColumn = column[j];
195          double value = elementByRow[j];
196          double solValue = solution[iColumn];
197          if (columnLower[iColumn] != columnUpper[iColumn]) {
198            if (solValue < 1.0 - integerTolerance && solValue > integerTolerance)
199              numberUnsatisfied++;
200            if (value != 1.0) {
201              oneRow = false;
202              break;
203            }
204          } else {
205            rhsValue -= value * floor(solValue + 0.5);
206          }
207        }
208        if (oneRow && rhsValue <= 1.0 + tolerance) {
209          if (!numberUnsatisfied) {
210            for (CoinBigIndex j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
211              int iColumn = column[j];
212              if (dsort[iColumn] > 1.0e50) {
213                dsort[iColumn] = 0;
214                nSort++;
215              }
216            }
217          }
218        }
219      }
220      // sort
221      CoinSort_2(dsort, dsort + numberColumns, sort);
222    } else {
223      // new way
224      for (i = 0; i < numberIntegers; i++) {
225        int iColumn = integerVariable[i];
226        if (upper[iColumn] > lower[iColumn]) {
227          if (!mark_ || !mark_[iColumn]) {
228            double distanceDown = solution[iColumn] - lower[iColumn];
229            double distanceUp = upper[iColumn] - solution[iColumn];
230            double distance = CoinMin(distanceDown, distanceUp);
231            if (distance > 0.001 && distance < 0.5) {
232              dsort[nSort] = distance;
233              sort[nSort++] = iColumn;
234            }
235          }
236        }
237      }
238      // sort
239      CoinSort_2(dsort, dsort + nSort, sort);
240      int n = 0;
241      double sum = 0.0;
242      for (int k = 0; k < nSort; k++) {
243        sum += dsort[k];
244        if (sum <= djTolerance_)
245          n = k;
246        else
247          break;
248      }
249      nSort = CoinMin(n, numberClean_ / 1000000);
250    }
251  } else {
252#define FIX_IF_LESS -0.1
253    // 3 in same row and sum <FIX_IF_LESS?
254    int numberRows = matrixByRow_.getNumRows();
255    const double *solution = model_->testSolution();
256    const int *column = matrixByRow_.getIndices();
257    const CoinBigIndex *rowStart = matrixByRow_.getVectorStarts();
258    const int *rowLength = matrixByRow_.getVectorLengths();
259    double bestSum = 1.0;
260    int nBest = -1;
261    int kRow = -1;
262    OsiSolverInterface *solver = model_->solver();
263    for (int i = 0; i < numberRows; i++) {
264      int numberUnsatisfied = 0;
265      double sum = 0.0;
266      for (CoinBigIndex j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
267        int iColumn = column[j];
268        if (solver->isInteger(iColumn)) {
269          double solValue = solution[iColumn];
270          if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
271            numberUnsatisfied++;
272            sum += solValue;
273          }
274        }
275      }
276      if (numberUnsatisfied >= 3 && sum < FIX_IF_LESS) {
277        // possible
278        if (numberUnsatisfied > nBest || (numberUnsatisfied == nBest && sum < bestSum)) {
279          nBest = numberUnsatisfied;
280          bestSum = sum;
281          kRow = i;
282        }
283      }
284    }
285    assert(nBest > 0);
286    for (CoinBigIndex j = rowStart[kRow]; j < rowStart[kRow] + rowLength[kRow]; j++) {
287      int iColumn = column[j];
288      if (solver->isInteger(iColumn)) {
289        double solValue = solution[iColumn];
290        if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
291          sort[nSort++] = iColumn;
292        }
293      }
294    }
295  }
296  OsiRowCut down;
297  down.setLb(-COIN_DBL_MAX);
298  double rhs = 0.0;
299  for (i = 0; i < nSort; i++) {
300    int iColumn = sort[i];
301    double distanceDown = solution[iColumn] - lower[iColumn];
302    double distanceUp = upper[iColumn] - solution[iColumn];
303    if (distanceDown < distanceUp) {
304      rhs += lower[iColumn];
305      dsort[i] = 1.0;
306    } else {
307      rhs -= upper[iColumn];
308      dsort[i] = -1.0;
309    }
310  }
311  down.setUb(rhs);
312  down.setRow(nSort, sort, dsort);
313  down.setEffectiveness(COIN_DBL_MAX); // so will persist
314  delete[] sort;
315  delete[] dsort;
316  // up is same - just with rhs changed
317  OsiRowCut up = down;
318  up.setLb(rhs + 1.0);
319  up.setUb(COIN_DBL_MAX);
320  // Say can fix one way
321  CbcCutBranchingObject *newObject = new CbcCutBranchingObject(model_, down, up, true);
322  if (model_->messageHandler()->logLevel() > 1)
323    printf("creating cut in CbcBranchCut\n");
324  return newObject;
325}
326/* Does a lot of the work,
327   Returns 0 if no good, 1 if dj, 2 if clean, 3 if both
328   10 if branching on ones away from bound
329*/
330int CbcBranchToFixLots::shallWe() const
331{
332  int returnCode = 0;
333  OsiSolverInterface *solver = model_->solver();
334  int numberRows = matrixByRow_.getNumRows();
335  //if (numberRows!=solver->getNumRows())
336  //return 0;
337  const double *solution = model_->testSolution();
338  const double *lower = solver->getColLower();
339  const double *upper = solver->getColUpper();
340  const double *dj = solver->getReducedCost();
341  int i;
342  int numberIntegers = model_->numberIntegers();
343  const int *integerVariable = model_->integerVariable();
344  if (numberClean_ > 1000000) {
345    int wanted = numberClean_ % 1000000;
346    int *sort = new int[numberIntegers];
347    double *dsort = new double[numberIntegers];
348    int nSort = 0;
349    for (i = 0; i < numberIntegers; i++) {
350      int iColumn = integerVariable[i];
351      if (upper[iColumn] > lower[iColumn]) {
352        if (!mark_ || !mark_[iColumn]) {
353          double distanceDown = solution[iColumn] - lower[iColumn];
354          double distanceUp = upper[iColumn] - solution[iColumn];
355          double distance = CoinMin(distanceDown, distanceUp);
356          if (distance > 0.001 && distance < 0.5) {
357            dsort[nSort] = distance;
358            sort[nSort++] = iColumn;
359          }
360        }
361      }
362    }
363    // sort
364    CoinSort_2(dsort, dsort + nSort, sort);
365    int n = 0;
366    double sum = 0.0;
367    for (int k = 0; k < nSort; k++) {
368      sum += dsort[k];
369      if (sum <= djTolerance_)
370        n = k;
371      else
372        break;
373    }
374    delete[] sort;
375    delete[] dsort;
376    return (n >= wanted) ? 10 : 0;
377  }
378  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
379  // make smaller ?
380  double tolerance = CoinMin(1.0e-8, integerTolerance);
381  // How many fixed are we aiming at
382  int wantedFixed = static_cast< int >(static_cast< double >(numberIntegers) * fractionFixed_);
383  if (djTolerance_ < 1.0e10) {
384    int nSort = 0;
385    int numberFixed = 0;
386    for (i = 0; i < numberIntegers; i++) {
387      int iColumn = integerVariable[i];
388      if (upper[iColumn] > lower[iColumn]) {
389        if (!mark_ || !mark_[iColumn]) {
390          if (solution[iColumn] < lower[iColumn] + tolerance) {
391            if (dj[iColumn] > djTolerance_) {
392              nSort++;
393            }
394          } else if (solution[iColumn] > upper[iColumn] - tolerance) {
395            if (dj[iColumn] < -djTolerance_) {
396              nSort++;
397            }
398          }
399        }
400      } else {
401        numberFixed++;
402      }
403    }
404    if (numberFixed + nSort < wantedFixed && !alwaysCreate_) {
405      returnCode = 0;
406    } else if (numberFixed < wantedFixed) {
407      returnCode = 1;
408    } else {
409      returnCode = 0;
410    }
411  }
412  if (numberClean_) {
413    // see how many rows clean
414    int i;
415    //const double * rowLower = solver->getRowLower();
416    const double *rowUpper = solver->getRowUpper();
417    // Row copy
418    const double *elementByRow = matrixByRow_.getElements();
419    const int *column = matrixByRow_.getIndices();
420    const CoinBigIndex *rowStart = matrixByRow_.getVectorStarts();
421    const int *rowLength = matrixByRow_.getVectorLengths();
422    const double *columnLower = solver->getColLower();
423    const double *columnUpper = solver->getColUpper();
424    const double *solution = solver->getColSolution();
425    int numberClean = 0;
426    bool someToDoYet = false;
427    int numberColumns = solver->getNumCols();
428    char *mark = new char[numberColumns];
429    int numberFixed = 0;
430    for (i = 0; i < numberColumns; i++) {
431      if (columnLower[i] != columnUpper[i]) {
432        mark[i] = 0;
433      } else {
434        mark[i] = 1;
435        numberFixed++;
436      }
437    }
438    int numberNewFixed = 0;
439    for (i = 0; i < numberRows; i++) {
440      double rhsValue = rowUpper[i];
441      bool oneRow = true;
442      // check elements
443      int numberUnsatisfied = 0;
444      for (CoinBigIndex j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
445        int iColumn = column[j];
446        double value = elementByRow[j];
447        double solValue = solution[iColumn];
448        if (columnLower[iColumn] != columnUpper[iColumn]) {
449          if (solValue < 1.0 - integerTolerance && solValue > integerTolerance)
450            numberUnsatisfied++;
451          if (value != 1.0) {
452            oneRow = false;
453            break;
454          }
455        } else {
456          rhsValue -= value * floor(solValue + 0.5);
457        }
458      }
459      if (oneRow && rhsValue <= 1.0 + tolerance) {
460        if (numberUnsatisfied) {
461          someToDoYet = true;
462        } else {
463          numberClean++;
464          for (CoinBigIndex j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
465            int iColumn = column[j];
466            if (columnLower[iColumn] != columnUpper[iColumn] && !mark[iColumn]) {
467              mark[iColumn] = 1;
468              numberNewFixed++;
469            }
470          }
471        }
472      }
473    }
474    delete[] mark;
475    //printf("%d clean, %d old fixed, %d new fixed\n",
476    //   numberClean,numberFixed,numberNewFixed);
477    if (someToDoYet && numberClean < numberClean_
478      && numberNewFixed + numberFixed < wantedFixed) {
479    } else if (numberFixed < wantedFixed) {
480      returnCode |= 2;
481    } else {
482    }
483  }
484  return returnCode;
485}
486double
487CbcBranchToFixLots::infeasibility(const OsiBranchingInformation * /*info*/,
488  int &preferredWay) const
489{
490  preferredWay = -1;
491  CbcNode *node = model_->currentNode();
492  int depth;
493  if (node)
494    depth = CoinMax(node->depth(), 0);
495  else
496    return 0.0;
497  if (depth_ < 0) {
498    return 0.0;
499  } else if (depth_ > 0) {
500    if ((depth % depth_) != 0)
501      return 0.0;
502  }
503  if (djTolerance_ != -1.234567) {
504    if (!shallWe())
505      return 0.0;
506    else
507      return 1.0e20;
508  } else {
509    // See if 3 in same row and sum <FIX_IF_LESS?
510    int numberRows = matrixByRow_.getNumRows();
511    const double *solution = model_->testSolution();
512    const int *column = matrixByRow_.getIndices();
513    const CoinBigIndex *rowStart = matrixByRow_.getVectorStarts();
514    const int *rowLength = matrixByRow_.getVectorLengths();
515    double bestSum = 1.0;
516    int nBest = -1;
517    OsiSolverInterface *solver = model_->solver();
518    for (int i = 0; i < numberRows; i++) {
519      int numberUnsatisfied = 0;
520      double sum = 0.0;
521      for (CoinBigIndex j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
522        int iColumn = column[j];
523        if (solver->isInteger(iColumn)) {
524          double solValue = solution[iColumn];
525          if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
526            numberUnsatisfied++;
527            sum += solValue;
528          }
529        }
530      }
531      if (numberUnsatisfied >= 3 && sum < FIX_IF_LESS) {
532        // possible
533        if (numberUnsatisfied > nBest || (numberUnsatisfied == nBest && sum < bestSum)) {
534          nBest = numberUnsatisfied;
535          bestSum = sum;
536        }
537      }
538    }
539    if (nBest > 0)
540      return 1.0e20;
541    else
542      return 0.0;
543  }
544}
545// Redoes data when sequence numbers change
546void CbcBranchToFixLots::redoSequenceEtc(CbcModel *model, int numberColumns, const int *originalColumns)
547{
548  model_ = model;
549  if (mark_) {
550    OsiSolverInterface *solver = model_->solver();
551    int numberColumnsNow = solver->getNumCols();
552    char *temp = new char[numberColumnsNow];
553    memset(temp, 0, numberColumnsNow);
554    for (int i = 0; i < numberColumns; i++) {
555      int j = originalColumns[i];
556      temp[i] = mark_[j];
557    }
558    delete[] mark_;
559    mark_ = temp;
560  }
561  OsiSolverInterface *solver = model_->solver();
562  matrixByRow_ = *solver->getMatrixByRow();
563}
564
565/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
566*/
Note: See TracBrowser for help on using the repository browser.