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

Last change on this file since 1899 was 1899, checked in by stefan, 6 years ago

fixup svn properties

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