source: trunk/Cbc/src/CbcSolverExpandKnapsack.cpp @ 1566

Last change on this file since 1566 was 1398, checked in by bjarni, 10 years ago

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

File size: 22.7 KB
Line 
1/* $Id: CbcSolverExpandKnapsack.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 CbcSolverExpandKnapsack.cpp
6
7  Returns OsiSolverInterface (User should delete)
8  On entry numberKnapsack is maximum number of Total entries
9
10  Expanding possibilities of x*y, where x*y are both integers, constructing
11  a knapsack constraint. Results in a tighter model.
12
13*/
14
15#include "CbcConfig.h"
16#include "CoinPragma.hpp"
17
18#include "OsiSolverInterface.hpp"
19
20#include "CglStored.hpp"
21
22#ifndef COIN_HAS_LINK
23#define COIN_HAS_LINK
24#endif
25#ifdef COIN_HAS_LINK
26#include "CbcLinked.hpp"
27#endif
28
29#ifdef COIN_HAS_LINK
30
31OsiSolverInterface *
32expandKnapsack(CoinModel & model, int * whichColumn, int * knapsackStart,
33               int * knapsackRow, int &numberKnapsack,
34               CglStored & stored, int logLevel,
35               int fixedPriority, int SOSPriority, CoinModel & tightenedModel)
36{
37    int maxTotal = numberKnapsack;
38    // load from coin model
39    OsiSolverLink *si = new OsiSolverLink();
40    OsiSolverInterface * finalModel = NULL;
41    si->setDefaultMeshSize(0.001);
42    // need some relative granularity
43    si->setDefaultBound(100.0);
44    si->setDefaultMeshSize(0.01);
45    si->setDefaultBound(100000.0);
46    si->setIntegerPriority(1000);
47    si->setBiLinearPriority(10000);
48    si->load(model, true, logLevel);
49    // get priorities
50    const int * priorities = model.priorities();
51    int numberColumns = model.numberColumns();
52    if (priorities) {
53        OsiObject ** objects = si->objects();
54        int numberObjects = si->numberObjects();
55        for (int iObj = 0; iObj < numberObjects; iObj++) {
56            int iColumn = objects[iObj]->columnNumber();
57            if (iColumn >= 0 && iColumn < numberColumns) {
58#ifndef NDEBUG
59                OsiSimpleInteger * obj =
60                    dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
61#endif
62                assert (obj);
63                int iPriority = priorities[iColumn];
64                if (iPriority > 0)
65                    objects[iObj]->setPriority(iPriority);
66            }
67        }
68        if (fixedPriority > 0) {
69            si->setFixedPriority(fixedPriority);
70        }
71        if (SOSPriority < 0)
72            SOSPriority = 100000;
73    }
74    CoinModel coinModel = *si->coinModel();
75    assert(coinModel.numberRows() > 0);
76    tightenedModel = coinModel;
77    int numberRows = coinModel.numberRows();
78    // Mark variables
79    int * whichKnapsack = new int [numberColumns];
80    int iRow, iColumn;
81    for (iColumn = 0; iColumn < numberColumns; iColumn++)
82        whichKnapsack[iColumn] = -1;
83    int kRow;
84    bool badModel = false;
85    // analyze
86    if (logLevel > 1) {
87        for (iRow = 0; iRow < numberRows; iRow++) {
88            /* Just obvious one at first
89            positive non unit coefficients
90            all integer
91            positive rowUpper
92            for now - linear (but further down in code may use nonlinear)
93            column bounds should be tight
94            */
95            //double lower = coinModel.getRowLower(iRow);
96            double upper = coinModel.getRowUpper(iRow);
97            if (upper < 1.0e10) {
98                CoinModelLink triple = coinModel.firstInRow(iRow);
99                bool possible = true;
100                int n = 0;
101                int n1 = 0;
102                while (triple.column() >= 0) {
103                    int iColumn = triple.column();
104                    const char *  el = coinModel.getElementAsString(iRow, iColumn);
105                    if (!strcmp("Numeric", el)) {
106                        if (coinModel.columnLower(iColumn) == coinModel.columnUpper(iColumn)) {
107                            triple = coinModel.next(triple);
108                            continue; // fixed
109                        }
110                        double value = coinModel.getElement(iRow, iColumn);
111                        if (value < 0.0) {
112                            possible = false;
113                        } else {
114                            n++;
115                            if (value == 1.0)
116                                n1++;
117                            if (coinModel.columnLower(iColumn) < 0.0)
118                                possible = false;
119                            if (!coinModel.isInteger(iColumn))
120                                possible = false;
121                            if (whichKnapsack[iColumn] >= 0)
122                                possible = false;
123                        }
124                    } else {
125                        possible = false; // non linear
126                    }
127                    triple = coinModel.next(triple);
128                }
129                if (n - n1 > 1 && possible) {
130                    double lower = coinModel.getRowLower(iRow);
131                    double upper = coinModel.getRowUpper(iRow);
132                    CoinModelLink triple = coinModel.firstInRow(iRow);
133                    while (triple.column() >= 0) {
134                        int iColumn = triple.column();
135                        lower -= coinModel.columnLower(iColumn) * triple.value();
136                        upper -= coinModel.columnLower(iColumn) * triple.value();
137                        triple = coinModel.next(triple);
138                    }
139                    printf("%d is possible %g <=", iRow, lower);
140                    // print
141                    triple = coinModel.firstInRow(iRow);
142                    while (triple.column() >= 0) {
143                        int iColumn = triple.column();
144                        if (coinModel.columnLower(iColumn) != coinModel.columnUpper(iColumn))
145                            printf(" (%d,el %g up %g)", iColumn, triple.value(),
146                                   coinModel.columnUpper(iColumn) - coinModel.columnLower(iColumn));
147                        triple = coinModel.next(triple);
148                    }
149                    printf(" <= %g\n", upper);
150                }
151            }
152        }
153    }
154    numberKnapsack = 0;
155    for (kRow = 0; kRow < numberRows; kRow++) {
156        iRow = kRow;
157        /* Just obvious one at first
158           positive non unit coefficients
159           all integer
160           positive rowUpper
161           for now - linear (but further down in code may use nonlinear)
162           column bounds should be tight
163        */
164        //double lower = coinModel.getRowLower(iRow);
165        double upper = coinModel.getRowUpper(iRow);
166        if (upper < 1.0e10) {
167            CoinModelLink triple = coinModel.firstInRow(iRow);
168            bool possible = true;
169            int n = 0;
170            int n1 = 0;
171            while (triple.column() >= 0) {
172                int iColumn = triple.column();
173                const char *  el = coinModel.getElementAsString(iRow, iColumn);
174                if (!strcmp("Numeric", el)) {
175                    if (coinModel.columnLower(iColumn) == coinModel.columnUpper(iColumn)) {
176                        triple = coinModel.next(triple);
177                        continue; // fixed
178                    }
179                    double value = coinModel.getElement(iRow, iColumn);
180                    if (value < 0.0) {
181                        possible = false;
182                    } else {
183                        n++;
184                        if (value == 1.0)
185                            n1++;
186                        if (coinModel.columnLower(iColumn) < 0.0)
187                            possible = false;
188                        if (!coinModel.isInteger(iColumn))
189                            possible = false;
190                        if (whichKnapsack[iColumn] >= 0)
191                            possible = false;
192                    }
193                } else {
194                    possible = false; // non linear
195                }
196                triple = coinModel.next(triple);
197            }
198            if (n - n1 > 1 && possible) {
199                // try
200                CoinModelLink triple = coinModel.firstInRow(iRow);
201                while (triple.column() >= 0) {
202                    int iColumn = triple.column();
203                    if (coinModel.columnLower(iColumn) != coinModel.columnUpper(iColumn))
204                        whichKnapsack[iColumn] = numberKnapsack;
205                    triple = coinModel.next(triple);
206                }
207                knapsackRow[numberKnapsack++] = iRow;
208            }
209        }
210    }
211    if (logLevel > 0)
212        printf("%d out of %d candidate rows are possible\n", numberKnapsack, numberRows);
213    // Check whether we can get rid of nonlinearities
214    /* mark rows
215       -2 in knapsack and other variables
216       -1 not involved
217       n only in knapsack n
218    */
219    int * markRow = new int [numberRows];
220    for (iRow = 0; iRow < numberRows; iRow++)
221        markRow[iRow] = -1;
222    int canDo = 1; // OK and linear
223    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
224        CoinModelLink triple = coinModel.firstInColumn(iColumn);
225        int iKnapsack = whichKnapsack[iColumn];
226        bool linear = true;
227        // See if quadratic objective
228        const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
229        if (strcmp(expr, "Numeric")) {
230            linear = false;
231        }
232        while (triple.row() >= 0) {
233            int iRow = triple.row();
234            if (iKnapsack >= 0) {
235                if (markRow[iRow] == -1) {
236                    markRow[iRow] = iKnapsack;
237                } else if (markRow[iRow] != iKnapsack) {
238                    markRow[iRow] = -2;
239                }
240            }
241            const char * expr = coinModel.getElementAsString(iRow, iColumn);
242            if (strcmp(expr, "Numeric")) {
243                linear = false;
244            }
245            triple = coinModel.next(triple);
246        }
247        if (!linear) {
248            if (whichKnapsack[iColumn] < 0) {
249                canDo = 0;
250                break;
251            } else {
252                canDo = 2;
253            }
254        }
255    }
256    int * markKnapsack = NULL;
257    double * coefficient = NULL;
258    double * linear = NULL;
259    int * whichRow = NULL;
260    int * lookupRow = NULL;
261    badModel = (canDo == 0);
262    if (numberKnapsack && canDo) {
263        /* double check - OK if
264           no nonlinear
265           nonlinear only on columns in knapsack
266           nonlinear only on columns in knapsack * ONE other - same for all in knapsack
267           AND that is only row connected to knapsack
268           (theoretically could split knapsack if two other and small numbers)
269           also ONE could be ONE expression - not just a variable
270        */
271        int iKnapsack;
272        markKnapsack = new int [numberKnapsack];
273        coefficient = new double [numberKnapsack];
274        linear = new double [numberColumns];
275        for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++)
276            markKnapsack[iKnapsack] = -1;
277        if (canDo == 2) {
278            for (iRow = -1; iRow < numberRows; iRow++) {
279                int numberOdd;
280                CoinPackedMatrix * row = coinModel.quadraticRow(iRow, linear, numberOdd);
281                if (row) {
282                    // see if valid
283                    const double * element = row->getElements();
284                    const int * column = row->getIndices();
285                    const CoinBigIndex * columnStart = row->getVectorStarts();
286                    const int * columnLength = row->getVectorLengths();
287                    int numberLook = row->getNumCols();
288                    for (int i = 0; i < numberLook; i++) {
289                        int iKnapsack = whichKnapsack[i];
290                        if (iKnapsack < 0) {
291                            // might be able to swap - but for now can't have knapsack in
292                            for (int j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
293                                int iColumn = column[j];
294                                if (whichKnapsack[iColumn] >= 0) {
295                                    canDo = 0; // no good
296                                    badModel = true;
297                                    break;
298                                }
299                            }
300                        } else {
301                            // OK if in same knapsack - or maybe just one
302                            int marked = markKnapsack[iKnapsack];
303                            for (int j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
304                                int iColumn = column[j];
305                                if (whichKnapsack[iColumn] != iKnapsack && whichKnapsack[iColumn] >= 0) {
306                                    canDo = 0; // no good
307                                    badModel = true;
308                                    break;
309                                } else if (marked == -1) {
310                                    markKnapsack[iKnapsack] = iColumn;
311                                    marked = iColumn;
312                                    coefficient[iKnapsack] = element[j];
313                                    coinModel.associateElement(coinModel.columnName(iColumn), 1.0);
314                                } else if (marked != iColumn) {
315                                    badModel = true;
316                                    canDo = 0; // no good
317                                    break;
318                                } else {
319                                    // could manage with different coefficients - but for now ...
320                                    assert(coefficient[iKnapsack] == element[j]);
321                                }
322                            }
323                        }
324                    }
325                    delete row;
326                }
327            }
328        }
329        if (canDo) {
330            // for any rows which are cuts
331            whichRow = new int [numberRows];
332            lookupRow = new int [numberRows];
333            bool someNonlinear = false;
334            double maxCoefficient = 1.0;
335            for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
336                if (markKnapsack[iKnapsack] >= 0) {
337                    someNonlinear = true;
338                    int iColumn = markKnapsack[iKnapsack];
339                    maxCoefficient = CoinMax(maxCoefficient, fabs(coefficient[iKnapsack] * coinModel.columnUpper(iColumn)));
340                }
341            }
342            if (someNonlinear) {
343                // associate all columns to stop possible error messages
344                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
345                    coinModel.associateElement(coinModel.columnName(iColumn), 1.0);
346                }
347            }
348            ClpSimplex tempModel;
349            tempModel.loadProblem(coinModel);
350            // Create final model - first without knapsacks
351            int nCol = 0;
352            int nRow = 0;
353            for (iRow = 0; iRow < numberRows; iRow++) {
354                if (markRow[iRow] < 0) {
355                    lookupRow[iRow] = nRow;
356                    whichRow[nRow++] = iRow;
357                } else {
358                    lookupRow[iRow] = -1;
359                }
360            }
361            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
362                if (whichKnapsack[iColumn] < 0)
363                    whichColumn[nCol++] = iColumn;
364            }
365            ClpSimplex finalModelX(&tempModel, nRow, whichRow, nCol, whichColumn, false, false, false);
366            OsiClpSolverInterface finalModelY(&finalModelX, true);
367            finalModel = finalModelY.clone();
368            finalModelY.releaseClp();
369            // Put back priorities
370            const int * priorities = model.priorities();
371            if (priorities) {
372                finalModel->findIntegers(false);
373                OsiObject ** objects = finalModel->objects();
374                int numberObjects = finalModel->numberObjects();
375                for (int iObj = 0; iObj < numberObjects; iObj++) {
376                    int iColumn = objects[iObj]->columnNumber();
377                    if (iColumn >= 0 && iColumn < nCol) {
378#ifndef NDEBUG
379                        OsiSimpleInteger * obj =
380                            dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
381#endif
382                        assert (obj);
383                        int iPriority = priorities[whichColumn[iColumn]];
384                        if (iPriority > 0)
385                            objects[iObj]->setPriority(iPriority);
386                    }
387                }
388            }
389            for (iRow = 0; iRow < numberRows; iRow++) {
390                whichRow[iRow] = iRow;
391            }
392            int numberOther = finalModel->getNumCols();
393            int nLargest = 0;
394            int nelLargest = 0;
395            int nTotal = 0;
396            for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
397                iRow = knapsackRow[iKnapsack];
398                int nCreate = maxTotal;
399                int nelCreate = coinModel.expandKnapsack(iRow, nCreate, NULL, NULL, NULL, NULL);
400                if (nelCreate < 0)
401                    badModel = true;
402                nTotal += nCreate;
403                nLargest = CoinMax(nLargest, nCreate);
404                nelLargest = CoinMax(nelLargest, nelCreate);
405            }
406            if (nTotal > maxTotal)
407                badModel = true;
408            if (!badModel) {
409                // Now arrays for building
410                nelLargest = CoinMax(nelLargest, nLargest) + 1;
411                double * buildObj = new double [nLargest];
412                double * buildElement = new double [nelLargest];
413                int * buildStart = new int[nLargest+1];
414                int * buildRow = new int[nelLargest];
415                // alow for integers in knapsacks
416                OsiObject ** object = new OsiObject * [numberKnapsack+nTotal];
417                int nSOS = 0;
418                int nObj = numberKnapsack;
419                for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
420                    knapsackStart[iKnapsack] = finalModel->getNumCols();
421                    iRow = knapsackRow[iKnapsack];
422                    int nCreate = 10000;
423                    coinModel.expandKnapsack(iRow, nCreate, buildObj, buildStart, buildRow, buildElement);
424                    // Redo row numbers
425                    for (iColumn = 0; iColumn < nCreate; iColumn++) {
426                        for (int j = buildStart[iColumn]; j < buildStart[iColumn+1]; j++) {
427                            int jRow = buildRow[j];
428                            jRow = lookupRow[jRow];
429                            assert (jRow >= 0 && jRow < nRow);
430                            buildRow[j] = jRow;
431                        }
432                    }
433                    finalModel->addCols(nCreate, buildStart, buildRow, buildElement, NULL, NULL, buildObj);
434                    int numberFinal = finalModel->getNumCols();
435                    for (iColumn = numberOther; iColumn < numberFinal; iColumn++) {
436                        if (markKnapsack[iKnapsack] < 0) {
437                            finalModel->setColUpper(iColumn, maxCoefficient);
438                            finalModel->setInteger(iColumn);
439                        } else {
440                            finalModel->setColUpper(iColumn, maxCoefficient + 1.0);
441                            finalModel->setInteger(iColumn);
442                        }
443                        OsiSimpleInteger * sosObject = new OsiSimpleInteger(finalModel, iColumn);
444                        sosObject->setPriority(1000000);
445                        object[nObj++] = sosObject;
446                        buildRow[iColumn-numberOther] = iColumn;
447                        buildElement[iColumn-numberOther] = 1.0;
448                    }
449                    if (markKnapsack[iKnapsack] < 0) {
450                        // convexity row
451                        finalModel->addRow(numberFinal - numberOther, buildRow, buildElement, 1.0, 1.0);
452                    } else {
453                        int iColumn = markKnapsack[iKnapsack];
454                        int n = numberFinal - numberOther;
455                        buildRow[n] = iColumn;
456                        buildElement[n++] = -fabs(coefficient[iKnapsack]);
457                        // convexity row (sort of)
458                        finalModel->addRow(n, buildRow, buildElement, 0.0, 0.0);
459                        OsiSOS * sosObject = new OsiSOS(finalModel, n - 1, buildRow, NULL, 1);
460                        sosObject->setPriority(iKnapsack + SOSPriority);
461                        // Say not integral even if is (switch off heuristics)
462                        sosObject->setIntegerValued(false);
463                        object[nSOS++] = sosObject;
464                    }
465                    numberOther = numberFinal;
466                }
467                finalModel->addObjects(nObj, object);
468                for (iKnapsack = 0; iKnapsack < nObj; iKnapsack++)
469                    delete object[iKnapsack];
470                delete [] object;
471                // Can we move any rows to cuts
472                const int * cutMarker = coinModel.cutMarker();
473                if (cutMarker && 0) {
474                    printf("AMPL CUTS OFF until global cuts fixed\n");
475                    cutMarker = NULL;
476                }
477                if (cutMarker) {
478                    // Row copy
479                    const CoinPackedMatrix * matrixByRow = finalModel->getMatrixByRow();
480                    const double * elementByRow = matrixByRow->getElements();
481                    const int * column = matrixByRow->getIndices();
482                    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
483                    const int * rowLength = matrixByRow->getVectorLengths();
484
485                    const double * rowLower = finalModel->getRowLower();
486                    const double * rowUpper = finalModel->getRowUpper();
487                    int nDelete = 0;
488                    for (iRow = 0; iRow < numberRows; iRow++) {
489                        if (cutMarker[iRow] && lookupRow[iRow] >= 0) {
490                            int jRow = lookupRow[iRow];
491                            whichRow[nDelete++] = jRow;
492                            int start = rowStart[jRow];
493                            stored.addCut(rowLower[jRow], rowUpper[jRow],
494                                          rowLength[jRow], column + start, elementByRow + start);
495                        }
496                    }
497                    finalModel->deleteRows(nDelete, whichRow);
498                }
499                knapsackStart[numberKnapsack] = finalModel->getNumCols();
500                delete [] buildObj;
501                delete [] buildElement;
502                delete [] buildStart;
503                delete [] buildRow;
504                finalModel->writeMps("full");
505            }
506        }
507    }
508    delete [] whichKnapsack;
509    delete [] markRow;
510    delete [] markKnapsack;
511    delete [] coefficient;
512    delete [] linear;
513    delete [] whichRow;
514    delete [] lookupRow;
515    delete si;
516    si = NULL;
517    if (!badModel && finalModel) {
518        finalModel->setDblParam(OsiObjOffset, coinModel.objectiveOffset());
519        return finalModel;
520    } else {
521        delete finalModel;
522        printf("can't make knapsacks - did you set fixedPriority (extra1)\n");
523        return NULL;
524    }
525}
526#endif  //COIN_HAS_LINK
527
528
Note: See TracBrowser for help on using the repository browser.