source: trunk/Cbc/src/CbcSolverAnalyze.cpp @ 1422

Last change on this file since 1422 was 1395, checked in by bjarni, 9 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.