source: trunk/Cbc/src/CbcSolverExpandKnapsack.cpp

Last change on this file was 2467, checked in by unxusr, 6 months ago

spaces after angles

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.3 KB
Line 
1/* $Id: CbcSolverExpandKnapsack.cpp 2467 2019-01-03 21:26:29Z tkr $ */
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 = 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 (CoinBigIndex 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 (CoinBigIndex 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 = dynamic_cast< OsiSimpleInteger * >(objects[iObj]);
380#endif
381            assert(obj);
382            int iPriority = priorities[whichColumn[iColumn]];
383            if (iPriority > 0)
384              objects[iObj]->setPriority(iPriority);
385          }
386        }
387      }
388      for (iRow = 0; iRow < numberRows; iRow++) {
389        whichRow[iRow] = iRow;
390      }
391      int numberOther = finalModel->getNumCols();
392      int nLargest = 0;
393      int nelLargest = 0;
394      int nTotal = 0;
395      for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
396        iRow = knapsackRow[iKnapsack];
397        int nCreate = maxTotal;
398        int nelCreate = coinModel.expandKnapsack(iRow, nCreate, NULL, NULL, NULL, NULL);
399        if (nelCreate < 0)
400          badModel = true;
401        nTotal += nCreate;
402        nLargest = CoinMax(nLargest, nCreate);
403        nelLargest = CoinMax(nelLargest, nelCreate);
404      }
405      if (nTotal > maxTotal)
406        badModel = true;
407      if (!badModel) {
408        // Now arrays for building
409        nelLargest = CoinMax(nelLargest, nLargest) + 1;
410        double *buildObj = new double[nLargest];
411        double *buildElement = new double[nelLargest];
412        CoinBigIndex *buildStart = new CoinBigIndex[nLargest + 1];
413        int *buildRow = new int[nelLargest];
414        // alow for integers in knapsacks
415        OsiObject **object = new OsiObject *[numberKnapsack + nTotal];
416        int nSOS = 0;
417        int nObj = numberKnapsack;
418        for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
419          knapsackStart[iKnapsack] = finalModel->getNumCols();
420          iRow = knapsackRow[iKnapsack];
421          int nCreate = 10000;
422          coinModel.expandKnapsack(iRow, nCreate, buildObj, buildStart, buildRow, buildElement);
423          // Redo row numbers
424          for (iColumn = 0; iColumn < nCreate; iColumn++) {
425            for (CoinBigIndex j = buildStart[iColumn]; j < buildStart[iColumn + 1]; j++) {
426              int jRow = buildRow[j];
427              jRow = lookupRow[jRow];
428              assert(jRow >= 0 && jRow < nRow);
429              buildRow[j] = jRow;
430            }
431          }
432          finalModel->addCols(nCreate, buildStart, buildRow, buildElement, NULL, NULL, buildObj);
433          int numberFinal = finalModel->getNumCols();
434          for (iColumn = numberOther; iColumn < numberFinal; iColumn++) {
435            if (markKnapsack[iKnapsack] < 0) {
436              finalModel->setColUpper(iColumn, maxCoefficient);
437              finalModel->setInteger(iColumn);
438            } else {
439              finalModel->setColUpper(iColumn, maxCoefficient + 1.0);
440              finalModel->setInteger(iColumn);
441            }
442            OsiSimpleInteger *sosObject = new OsiSimpleInteger(finalModel, iColumn);
443            sosObject->setPriority(1000000);
444            object[nObj++] = sosObject;
445            buildRow[iColumn - numberOther] = iColumn;
446            buildElement[iColumn - numberOther] = 1.0;
447          }
448          if (markKnapsack[iKnapsack] < 0) {
449            // convexity row
450            finalModel->addRow(numberFinal - numberOther, buildRow, buildElement, 1.0, 1.0);
451          } else {
452            int iColumn = markKnapsack[iKnapsack];
453            int n = numberFinal - numberOther;
454            buildRow[n] = iColumn;
455            buildElement[n++] = -fabs(coefficient[iKnapsack]);
456            // convexity row (sort of)
457            finalModel->addRow(n, buildRow, buildElement, 0.0, 0.0);
458            OsiSOS *sosObject = new OsiSOS(finalModel, n - 1, buildRow, NULL, 1);
459            sosObject->setPriority(iKnapsack + SOSPriority);
460            // Say not integral even if is (switch off heuristics)
461            sosObject->setIntegerValued(false);
462            object[nSOS++] = sosObject;
463          }
464          numberOther = numberFinal;
465        }
466        finalModel->addObjects(nObj, object);
467        for (iKnapsack = 0; iKnapsack < nObj; iKnapsack++)
468          delete object[iKnapsack];
469        delete[] object;
470        // Can we move any rows to cuts
471        const int *cutMarker = coinModel.cutMarker();
472        if (cutMarker && 0) {
473          printf("AMPL CUTS OFF until global cuts fixed\n");
474          cutMarker = NULL;
475        }
476        if (cutMarker) {
477          // Row copy
478          const CoinPackedMatrix *matrixByRow = finalModel->getMatrixByRow();
479          const double *elementByRow = matrixByRow->getElements();
480          const int *column = matrixByRow->getIndices();
481          const CoinBigIndex *rowStart = matrixByRow->getVectorStarts();
482          const int *rowLength = matrixByRow->getVectorLengths();
483
484          const double *rowLower = finalModel->getRowLower();
485          const double *rowUpper = finalModel->getRowUpper();
486          int nDelete = 0;
487          for (iRow = 0; iRow < numberRows; iRow++) {
488            if (cutMarker[iRow] && lookupRow[iRow] >= 0) {
489              int jRow = lookupRow[iRow];
490              whichRow[nDelete++] = jRow;
491              CoinBigIndex start = rowStart[jRow];
492              stored.addCut(rowLower[jRow], rowUpper[jRow],
493                rowLength[jRow], column + start, elementByRow + start);
494            }
495          }
496          finalModel->deleteRows(nDelete, whichRow);
497        }
498        knapsackStart[numberKnapsack] = finalModel->getNumCols();
499        delete[] buildObj;
500        delete[] buildElement;
501        delete[] buildStart;
502        delete[] buildRow;
503        finalModel->writeMps("full");
504      }
505    }
506  }
507  delete[] whichKnapsack;
508  delete[] markRow;
509  delete[] markKnapsack;
510  delete[] coefficient;
511  delete[] linear;
512  delete[] whichRow;
513  delete[] lookupRow;
514  delete si;
515  si = NULL;
516  if (!badModel && finalModel) {
517    finalModel->setDblParam(OsiObjOffset, coinModel.objectiveOffset());
518    return finalModel;
519  } else {
520    delete finalModel;
521    printf("can't make knapsacks - did you set fixedPriority (extra1)\n");
522    return NULL;
523  }
524}
525#endif //COIN_HAS_LINK
526
527// Fills in original solution (coinModel length)
528void afterKnapsack(const CoinModel &coinModel2, const int *whichColumn, const int *knapsackStart,
529  const int *knapsackRow, int numberKnapsack,
530  const double *knapsackSolution, double *solution, int logLevel)
531{
532  CoinModel coinModel = coinModel2;
533  int numberColumns = coinModel.numberColumns();
534  int iColumn;
535  // associate all columns to stop possible error messages
536  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
537    coinModel.associateElement(coinModel.columnName(iColumn), 1.0);
538  }
539  CoinZeroN(solution, numberColumns);
540  int nCol = knapsackStart[0];
541  for (iColumn = 0; iColumn < nCol; iColumn++) {
542    int jColumn = whichColumn[iColumn];
543    solution[jColumn] = knapsackSolution[iColumn];
544  }
545  int *buildRow = new int[numberColumns]; // wild overkill
546  double *buildElement = new double[numberColumns];
547  int iKnapsack;
548  for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
549    int k = -1;
550    for (iColumn = knapsackStart[iKnapsack]; iColumn < knapsackStart[iKnapsack + 1]; iColumn++) {
551      if (knapsackSolution[iColumn] > 1.0e-5) {
552        if (k >= 0) {
553          printf("Two nonzero values for knapsack %d at (%d,%g) and (%d,%g)\n", iKnapsack,
554            k, knapsackSolution[k], iColumn, knapsackSolution[iColumn]);
555          abort();
556        }
557        k = iColumn;
558        assert(fabs(floor(knapsackSolution[iColumn] + 0.5) - knapsackSolution[iColumn]) < 1.0e-5);
559      }
560    }
561    if (k >= 0) {
562      int iRow = knapsackRow[iKnapsack];
563      int nCreate = 10000;
564      int nel = coinModel.expandKnapsack(iRow, nCreate, NULL, NULL, buildRow, buildElement, k - knapsackStart[iKnapsack]);
565      assert(nel);
566      if (logLevel > 0)
567        printf("expanded column %d in knapsack %d has %d nonzero entries:\n",
568          k - knapsackStart[iKnapsack], iKnapsack, nel);
569      for (int i = 0; i < nel; i++) {
570        int jColumn = buildRow[i];
571        double value = buildElement[i];
572        if (logLevel > 0)
573          printf("%d - original %d has value %g\n", i, jColumn, value);
574        solution[jColumn] = value;
575      }
576    }
577  }
578  delete[] buildRow;
579  delete[] buildElement;
580#if 0
581   for (iColumn=0;iColumn<numberColumns;iColumn++) {
582      if (solution[iColumn]>1.0e-5&&coinModel.isInteger(iColumn))
583         printf("%d %g\n",iColumn,solution[iColumn]);
584   }
585#endif
586}
587
588/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
589*/
Note: See TracBrowser for help on using the repository browser.