source: stable/2.5/Cbc/src/CbcSolverAnalyze.cpp @ 1510

Last change on this file since 1510 was 1395, checked in by bjarni, 10 years ago

Extracted analyze() from CbcSolver?.cpp and placed it in CbcSolverAnalyze?.cpp, updated libCbc.vcproj (v9) to include CbcSolverAnalyze?.cpp

File size: 14.0 KB
Line 
1/* $Id: CbcSolverAnalyze.cpp 1240 2009-10-02 18:41:44Z forrest $ */
2// Copyright (C) 2007, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5/*! \file CbcSolverAnalyze.cpp
6
7  Look to see if a constraint is all-integer (variables & coeffs), or could be
8  all integer. Consider whether one or two continuous variables can be declared
9  integer. John's comment is that integer preprocessing might do a better job,
10  so we should consider whether this routine should stay.
11
12  No hurry to get rid of it, in my opinion  -- lh, 091210 --
13
14*/
15
16#include "CbcConfig.h"
17#include "CoinPragma.hpp"
18
19#include "OsiClpSolverInterface.hpp"
20
21#include "ClpMessage.hpp"
22
23#include "CbcModel.hpp"
24
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                        int 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        delete [] ignore;
285        //increment=0.0;
286        if (!numberChanged) {
287            delete [] changed;
288            delete solver;
289            return NULL;
290        } else {
291            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
292                if (solver->isInteger(iColumn))
293                    solverMod->setInteger(iColumn);
294            }
295            delete solver;
296            return changed;
297        }
298    } else {
299        //if (!noPrinting_) {
300        //if (numberChanged)
301        //  printf(" and %d variables could be made integer\n",numberChanged);
302        //else
303        //  printf("\n");
304        //}
305        // just get increment
306        int logLevel = generalMessageHandler->logLevel();
307        CbcModel model(*solver);
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
Note: See TracBrowser for help on using the repository browser.