source: trunk/Cbc/src/CbcSolverAnalyze.cpp

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

script to format sources

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.4 KB
Line 
1/* $Id: CbcSolverAnalyze.cpp 2465 2019-01-03 19:26:52Z forrest $ */
2// Copyright (C) 2007, 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/*! \file CbcSolverAnalyze.cpp
7
8  Look to see if a constraint is all-integer (variables & coeffs), or could be
9  all integer. Consider whether one or two continuous variables can be declared
10  integer. John's comment is that integer preprocessing might do a better job,
11  so we should consider whether this routine should stay.
12
13  No hurry to get rid of it, in my opinion  -- lh, 091210 --
14
15*/
16
17#include "CbcConfig.h"
18#include "CoinPragma.hpp"
19
20#include "OsiClpSolverInterface.hpp"
21
22#include "ClpMessage.hpp"
23
24#include "CbcModel.hpp"
25
26#ifndef CBC_OTHER_SOLVER
27
28int *analyze(OsiClpSolverInterface *solverMod, int &numberChanged,
29  double &increment, bool changeInt,
30  CoinMessageHandler *generalMessageHandler, bool noPrinting)
31{
32  bool noPrinting_ = noPrinting;
33  OsiSolverInterface *solver = solverMod->clone();
34  char generalPrint[200];
35  if (0) {
36    // just get increment
37    CbcModel model(*solver);
38    model.analyzeObjective();
39    double increment2 = model.getCutoffIncrement();
40    printf("initial cutoff increment %g\n", increment2);
41  }
42  const double *objective = solver->getObjCoefficients();
43  const double *lower = solver->getColLower();
44  const double *upper = solver->getColUpper();
45  int numberColumns = solver->getNumCols();
46  int numberRows = solver->getNumRows();
47  double direction = solver->getObjSense();
48  int iRow, iColumn;
49
50  // Row copy
51  CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
52  const double *elementByRow = matrixByRow.getElements();
53  const int *column = matrixByRow.getIndices();
54  const CoinBigIndex *rowStart = matrixByRow.getVectorStarts();
55  const int *rowLength = matrixByRow.getVectorLengths();
56
57  // Column copy
58  CoinPackedMatrix matrixByCol(*solver->getMatrixByCol());
59  const double *element = matrixByCol.getElements();
60  const int *row = matrixByCol.getIndices();
61  const CoinBigIndex *columnStart = matrixByCol.getVectorStarts();
62  const int *columnLength = matrixByCol.getVectorLengths();
63
64  const double *rowLower = solver->getRowLower();
65  const double *rowUpper = solver->getRowUpper();
66
67  char *ignore = new char[numberRows];
68  int *changed = new int[numberColumns];
69  int *which = new int[numberRows];
70  double *changeRhs = new double[numberRows];
71  memset(changeRhs, 0, numberRows * sizeof(double));
72  memset(ignore, 0, numberRows);
73  numberChanged = 0;
74  int numberInteger = 0;
75  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
76    if (upper[iColumn] > lower[iColumn] + 1.0e-8 && solver->isInteger(iColumn))
77      numberInteger++;
78  }
79  bool finished = false;
80  while (!finished) {
81    int saveNumberChanged = numberChanged;
82    for (iRow = 0; iRow < numberRows; iRow++) {
83      int numberContinuous = 0;
84      double value1 = 0.0, value2 = 0.0;
85      bool allIntegerCoeff = true;
86      double sumFixed = 0.0;
87      int jColumn1 = -1, jColumn2 = -1;
88      for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
89        int jColumn = column[j];
90        double value = elementByRow[j];
91        if (upper[jColumn] > lower[jColumn] + 1.0e-8) {
92          if (!solver->isInteger(jColumn)) {
93            if (numberContinuous == 0) {
94              jColumn1 = jColumn;
95              value1 = value;
96            } else {
97              jColumn2 = jColumn;
98              value2 = value;
99            }
100            numberContinuous++;
101          } else {
102            if (fabs(value - floor(value + 0.5)) > 1.0e-12)
103              allIntegerCoeff = false;
104          }
105        } else {
106          sumFixed += lower[jColumn] * value;
107        }
108      }
109      double low = rowLower[iRow];
110      if (low > -1.0e20) {
111        low -= sumFixed;
112        if (fabs(low - floor(low + 0.5)) > 1.0e-12)
113          allIntegerCoeff = false;
114      }
115      double up = rowUpper[iRow];
116      if (up < 1.0e20) {
117        up -= sumFixed;
118        if (fabs(up - floor(up + 0.5)) > 1.0e-12)
119          allIntegerCoeff = false;
120      }
121      if (!allIntegerCoeff)
122        continue; // can't do
123      if (numberContinuous == 1) {
124        // see if really integer
125        // This does not allow for complicated cases
126        if (low == up) {
127          if (fabs(value1) > 1.0e-3) {
128            value1 = 1.0 / value1;
129            if (fabs(value1 - floor(value1 + 0.5)) < 1.0e-12) {
130              // integer
131              changed[numberChanged++] = jColumn1;
132              solver->setInteger(jColumn1);
133              if (upper[jColumn1] > 1.0e20)
134                solver->setColUpper(jColumn1, 1.0e20);
135              if (lower[jColumn1] < -1.0e20)
136                solver->setColLower(jColumn1, -1.0e20);
137            }
138          }
139        } else {
140          if (fabs(value1) > 1.0e-3) {
141            value1 = 1.0 / value1;
142            if (fabs(value1 - floor(value1 + 0.5)) < 1.0e-12) {
143              // This constraint will not stop it being integer
144              ignore[iRow] = 1;
145            }
146          }
147        }
148      } else if (numberContinuous == 2) {
149        if (low == up) {
150          /* need general theory - for now just look at 2 cases -
151                       1 - +- 1 one in column and just costs i.e. matching objective
152                       2 - +- 1 two in column but feeds into G/L row which will try and minimize
153                    */
154          if (fabs(value1) == 1.0 && value1 * value2 == -1.0 && !lower[jColumn1]
155            && !lower[jColumn2]) {
156            int n = 0;
157            CoinBigIndex i;
158            double objChange = direction * (objective[jColumn1] + objective[jColumn2]);
159            double bound = CoinMin(upper[jColumn1], upper[jColumn2]);
160            bound = CoinMin(bound, 1.0e20);
161            for (i = columnStart[jColumn1]; i < columnStart[jColumn1] + columnLength[jColumn1]; i++) {
162              int jRow = row[i];
163              double value = element[i];
164              if (jRow != iRow) {
165                which[n++] = jRow;
166                changeRhs[jRow] = value;
167              }
168            }
169            for (i = columnStart[jColumn1]; i < columnStart[jColumn1] + columnLength[jColumn1]; i++) {
170              int jRow = row[i];
171              double value = element[i];
172              if (jRow != iRow) {
173                if (!changeRhs[jRow]) {
174                  which[n++] = jRow;
175                  changeRhs[jRow] = value;
176                } else {
177                  changeRhs[jRow] += value;
178                }
179              }
180            }
181            if (objChange >= 0.0) {
182              // see if all rows OK
183              bool good = true;
184              for (i = 0; i < n; i++) {
185                int jRow = which[i];
186                double value = changeRhs[jRow];
187                if (value) {
188                  value *= bound;
189                  if (rowLength[jRow] == 1) {
190                    if (value > 0.0) {
191                      double rhs = rowLower[jRow];
192                      if (rhs > 0.0) {
193                        double ratio = rhs / value;
194                        if (fabs(ratio - floor(ratio + 0.5)) > 1.0e-12)
195                          good = false;
196                      }
197                    } else {
198                      double rhs = rowUpper[jRow];
199                      if (rhs < 0.0) {
200                        double ratio = rhs / value;
201                        if (fabs(ratio - floor(ratio + 0.5)) > 1.0e-12)
202                          good = false;
203                      }
204                    }
205                  } else if (rowLength[jRow] == 2) {
206                    if (value > 0.0) {
207                      if (rowLower[jRow] > -1.0e20)
208                        good = false;
209                    } else {
210                      if (rowUpper[jRow] < 1.0e20)
211                        good = false;
212                    }
213                  } else {
214                    good = false;
215                  }
216                }
217              }
218              if (good) {
219                // both can be integer
220                changed[numberChanged++] = jColumn1;
221                solver->setInteger(jColumn1);
222                if (upper[jColumn1] > 1.0e20)
223                  solver->setColUpper(jColumn1, 1.0e20);
224                if (lower[jColumn1] < -1.0e20)
225                  solver->setColLower(jColumn1, -1.0e20);
226                changed[numberChanged++] = jColumn2;
227                solver->setInteger(jColumn2);
228                if (upper[jColumn2] > 1.0e20)
229                  solver->setColUpper(jColumn2, 1.0e20);
230                if (lower[jColumn2] < -1.0e20)
231                  solver->setColLower(jColumn2, -1.0e20);
232              }
233            }
234            // clear
235            for (i = 0; i < n; i++) {
236              changeRhs[which[i]] = 0.0;
237            }
238          }
239        }
240      }
241    }
242    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
243      if (upper[iColumn] > lower[iColumn] + 1.0e-8 && !solver->isInteger(iColumn)) {
244        double value;
245        value = upper[iColumn];
246        if (value < 1.0e20 && fabs(value - floor(value + 0.5)) > 1.0e-12)
247          continue;
248        value = lower[iColumn];
249        if (value > -1.0e20 && fabs(value - floor(value + 0.5)) > 1.0e-12)
250          continue;
251        bool integer = true;
252        for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
253          int iRow = row[j];
254          if (!ignore[iRow]) {
255            integer = false;
256            break;
257          }
258        }
259        if (integer) {
260          // integer
261          changed[numberChanged++] = iColumn;
262          solver->setInteger(iColumn);
263          if (upper[iColumn] > 1.0e20)
264            solver->setColUpper(iColumn, 1.0e20);
265          if (lower[iColumn] < -1.0e20)
266            solver->setColLower(iColumn, -1.0e20);
267        }
268      }
269    }
270    finished = numberChanged == saveNumberChanged;
271  }
272  delete[] which;
273  delete[] changeRhs;
274  delete[] ignore;
275  //if (numberInteger&&!noPrinting_)
276  //printf("%d integer variables",numberInteger);
277  if (changeInt) {
278    //if (!noPrinting_) {
279    //if (numberChanged)
280    //  printf(" and %d variables made integer\n",numberChanged);
281    //else
282    //  printf("\n");
283    //}
284    //increment=0.0;
285    if (!numberChanged) {
286      delete[] changed;
287      delete solver;
288      return NULL;
289    } else {
290      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
291        if (solver->isInteger(iColumn))
292          solverMod->setInteger(iColumn);
293      }
294      delete solver;
295      return changed;
296    }
297  } else {
298    //if (!noPrinting_) {
299    //if (numberChanged)
300    //  printf(" and %d variables could be made integer\n",numberChanged);
301    //else
302    //  printf("\n");
303    //}
304    // just get increment
305    int logLevel = generalMessageHandler->logLevel();
306    CbcModel model(*solver);
307    if (!model.defaultHandler())
308      model.passInMessageHandler(generalMessageHandler);
309    if (noPrinting_)
310      model.setLogLevel(0);
311    model.analyzeObjective();
312    generalMessageHandler->setLogLevel(logLevel);
313    double increment2 = model.getCutoffIncrement();
314    if (increment2 > increment && increment2 > 0.0) {
315      if (!noPrinting_) {
316        sprintf(generalPrint, "Cutoff increment increased from %g to %g", increment, increment2);
317        CoinMessages generalMessages = solverMod->getModelPtr()->messages();
318        generalMessageHandler->message(CLP_GENERAL, generalMessages)
319          << generalPrint
320          << CoinMessageEol;
321      }
322      increment = increment2;
323    }
324    delete solver;
325    numberChanged = 0;
326    delete[] changed;
327    return NULL;
328  }
329}
330#endif // ifndef CBC_OTHER_SOLVER
331
332/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
333*/
Note: See TracBrowser for help on using the repository browser.