Changeset 1398 for branches/sandbox


Ignore:
Timestamp:
Dec 10, 2009 11:46:08 PM (10 years ago)
Author:
bjarni
Message:

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

Location:
branches/sandbox/Cbc
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/Cbc/MSVisualStudio/v9/libCbc/libCbc.vcproj

    r1395 r1398  
    948948                        </File>
    949949                        <File
     950                                RelativePath="..\..\..\src\CbcSolverExpandKnapsack.cpp"
     951                                >
     952                        </File>
     953                        <File
    950954                                RelativePath="..\..\..\src\CbcSolverHeuristics.cpp"
    951955                                >
     
    13411345                        </File>
    13421346                        <File
     1347                                RelativePath="..\..\..\src\CbcSolverExpandKnapsack.hpp"
     1348                                >
     1349                        </File>
     1350                        <File
    13431351                                RelativePath="..\..\..\src\CbcSolverHeuristics.hpp"
    13441352                                >
  • branches/sandbox/Cbc/src/CbcSolver.cpp

    r1395 r1398  
    202202
    203203#include "CbcSolverAnalyze.hpp"
     204#include "CbcSolverExpandKnapsack.hpp"
    204205
    205206#include "CbcSolver.hpp"
     
    797798}
    798799#endif
    799 
    800 
    801 
    802 #ifdef COIN_HAS_LINK
    803 /*  Returns OsiSolverInterface (User should delete)
    804     On entry numberKnapsack is maximum number of Total entries
    805 
    806     Expanding possibilities of x*y, where x*y are both integers, constructing
    807     a knapsack constraint. Results in a tighter model.
    808 */
    809 static OsiSolverInterface *
    810 expandKnapsack(CoinModel & model, int * whichColumn, int * knapsackStart,
    811                int * knapsackRow, int &numberKnapsack,
    812                CglStored & stored, int logLevel,
    813                int fixedPriority, int SOSPriority, CoinModel & tightenedModel)
    814 {
    815     int maxTotal = numberKnapsack;
    816     // load from coin model
    817     OsiSolverLink *si = new OsiSolverLink();
    818     OsiSolverInterface * finalModel = NULL;
    819     si->setDefaultMeshSize(0.001);
    820     // need some relative granularity
    821     si->setDefaultBound(100.0);
    822     si->setDefaultMeshSize(0.01);
    823     si->setDefaultBound(100000.0);
    824     si->setIntegerPriority(1000);
    825     si->setBiLinearPriority(10000);
    826     si->load(model, true, logLevel);
    827     // get priorities
    828     const int * priorities = model.priorities();
    829     int numberColumns = model.numberColumns();
    830     if (priorities) {
    831         OsiObject ** objects = si->objects();
    832         int numberObjects = si->numberObjects();
    833         for (int iObj = 0; iObj < numberObjects; iObj++) {
    834             int iColumn = objects[iObj]->columnNumber();
    835             if (iColumn >= 0 && iColumn < numberColumns) {
    836 #ifndef NDEBUG
    837                 OsiSimpleInteger * obj =
    838                     dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
    839 #endif
    840                 assert (obj);
    841                 int iPriority = priorities[iColumn];
    842                 if (iPriority > 0)
    843                     objects[iObj]->setPriority(iPriority);
    844             }
    845         }
    846         if (fixedPriority > 0) {
    847             si->setFixedPriority(fixedPriority);
    848         }
    849         if (SOSPriority < 0)
    850             SOSPriority = 100000;
    851     }
    852     CoinModel coinModel = *si->coinModel();
    853     assert(coinModel.numberRows() > 0);
    854     tightenedModel = coinModel;
    855     int numberRows = coinModel.numberRows();
    856     // Mark variables
    857     int * whichKnapsack = new int [numberColumns];
    858     int iRow, iColumn;
    859     for (iColumn = 0; iColumn < numberColumns; iColumn++)
    860         whichKnapsack[iColumn] = -1;
    861     int kRow;
    862     bool badModel = false;
    863     // analyze
    864     if (logLevel > 1) {
    865         for (iRow = 0; iRow < numberRows; iRow++) {
    866             /* Just obvious one at first
    867             positive non unit coefficients
    868             all integer
    869             positive rowUpper
    870             for now - linear (but further down in code may use nonlinear)
    871             column bounds should be tight
    872             */
    873             //double lower = coinModel.getRowLower(iRow);
    874             double upper = coinModel.getRowUpper(iRow);
    875             if (upper < 1.0e10) {
    876                 CoinModelLink triple = coinModel.firstInRow(iRow);
    877                 bool possible = true;
    878                 int n = 0;
    879                 int n1 = 0;
    880                 while (triple.column() >= 0) {
    881                     int iColumn = triple.column();
    882                     const char *  el = coinModel.getElementAsString(iRow, iColumn);
    883                     if (!strcmp("Numeric", el)) {
    884                         if (coinModel.columnLower(iColumn) == coinModel.columnUpper(iColumn)) {
    885                             triple = coinModel.next(triple);
    886                             continue; // fixed
    887                         }
    888                         double value = coinModel.getElement(iRow, iColumn);
    889                         if (value < 0.0) {
    890                             possible = false;
    891                         } else {
    892                             n++;
    893                             if (value == 1.0)
    894                                 n1++;
    895                             if (coinModel.columnLower(iColumn) < 0.0)
    896                                 possible = false;
    897                             if (!coinModel.isInteger(iColumn))
    898                                 possible = false;
    899                             if (whichKnapsack[iColumn] >= 0)
    900                                 possible = false;
    901                         }
    902                     } else {
    903                         possible = false; // non linear
    904                     }
    905                     triple = coinModel.next(triple);
    906                 }
    907                 if (n - n1 > 1 && possible) {
    908                     double lower = coinModel.getRowLower(iRow);
    909                     double upper = coinModel.getRowUpper(iRow);
    910                     CoinModelLink triple = coinModel.firstInRow(iRow);
    911                     while (triple.column() >= 0) {
    912                         int iColumn = triple.column();
    913                         lower -= coinModel.columnLower(iColumn) * triple.value();
    914                         upper -= coinModel.columnLower(iColumn) * triple.value();
    915                         triple = coinModel.next(triple);
    916                     }
    917                     printf("%d is possible %g <=", iRow, lower);
    918                     // print
    919                     triple = coinModel.firstInRow(iRow);
    920                     while (triple.column() >= 0) {
    921                         int iColumn = triple.column();
    922                         if (coinModel.columnLower(iColumn) != coinModel.columnUpper(iColumn))
    923                             printf(" (%d,el %g up %g)", iColumn, triple.value(),
    924                                    coinModel.columnUpper(iColumn) - coinModel.columnLower(iColumn));
    925                         triple = coinModel.next(triple);
    926                     }
    927                     printf(" <= %g\n", upper);
    928                 }
    929             }
    930         }
    931     }
    932     numberKnapsack = 0;
    933     for (kRow = 0; kRow < numberRows; kRow++) {
    934         iRow = kRow;
    935         /* Just obvious one at first
    936            positive non unit coefficients
    937            all integer
    938            positive rowUpper
    939            for now - linear (but further down in code may use nonlinear)
    940            column bounds should be tight
    941         */
    942         //double lower = coinModel.getRowLower(iRow);
    943         double upper = coinModel.getRowUpper(iRow);
    944         if (upper < 1.0e10) {
    945             CoinModelLink triple = coinModel.firstInRow(iRow);
    946             bool possible = true;
    947             int n = 0;
    948             int n1 = 0;
    949             while (triple.column() >= 0) {
    950                 int iColumn = triple.column();
    951                 const char *  el = coinModel.getElementAsString(iRow, iColumn);
    952                 if (!strcmp("Numeric", el)) {
    953                     if (coinModel.columnLower(iColumn) == coinModel.columnUpper(iColumn)) {
    954                         triple = coinModel.next(triple);
    955                         continue; // fixed
    956                     }
    957                     double value = coinModel.getElement(iRow, iColumn);
    958                     if (value < 0.0) {
    959                         possible = false;
    960                     } else {
    961                         n++;
    962                         if (value == 1.0)
    963                             n1++;
    964                         if (coinModel.columnLower(iColumn) < 0.0)
    965                             possible = false;
    966                         if (!coinModel.isInteger(iColumn))
    967                             possible = false;
    968                         if (whichKnapsack[iColumn] >= 0)
    969                             possible = false;
    970                     }
    971                 } else {
    972                     possible = false; // non linear
    973                 }
    974                 triple = coinModel.next(triple);
    975             }
    976             if (n - n1 > 1 && possible) {
    977                 // try
    978                 CoinModelLink triple = coinModel.firstInRow(iRow);
    979                 while (triple.column() >= 0) {
    980                     int iColumn = triple.column();
    981                     if (coinModel.columnLower(iColumn) != coinModel.columnUpper(iColumn))
    982                         whichKnapsack[iColumn] = numberKnapsack;
    983                     triple = coinModel.next(triple);
    984                 }
    985                 knapsackRow[numberKnapsack++] = iRow;
    986             }
    987         }
    988     }
    989     if (logLevel > 0)
    990         printf("%d out of %d candidate rows are possible\n", numberKnapsack, numberRows);
    991     // Check whether we can get rid of nonlinearities
    992     /* mark rows
    993        -2 in knapsack and other variables
    994        -1 not involved
    995        n only in knapsack n
    996     */
    997     int * markRow = new int [numberRows];
    998     for (iRow = 0; iRow < numberRows; iRow++)
    999         markRow[iRow] = -1;
    1000     int canDo = 1; // OK and linear
    1001     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    1002         CoinModelLink triple = coinModel.firstInColumn(iColumn);
    1003         int iKnapsack = whichKnapsack[iColumn];
    1004         bool linear = true;
    1005         // See if quadratic objective
    1006         const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
    1007         if (strcmp(expr, "Numeric")) {
    1008             linear = false;
    1009         }
    1010         while (triple.row() >= 0) {
    1011             int iRow = triple.row();
    1012             if (iKnapsack >= 0) {
    1013                 if (markRow[iRow] == -1) {
    1014                     markRow[iRow] = iKnapsack;
    1015                 } else if (markRow[iRow] != iKnapsack) {
    1016                     markRow[iRow] = -2;
    1017                 }
    1018             }
    1019             const char * expr = coinModel.getElementAsString(iRow, iColumn);
    1020             if (strcmp(expr, "Numeric")) {
    1021                 linear = false;
    1022             }
    1023             triple = coinModel.next(triple);
    1024         }
    1025         if (!linear) {
    1026             if (whichKnapsack[iColumn] < 0) {
    1027                 canDo = 0;
    1028                 break;
    1029             } else {
    1030                 canDo = 2;
    1031             }
    1032         }
    1033     }
    1034     int * markKnapsack = NULL;
    1035     double * coefficient = NULL;
    1036     double * linear = NULL;
    1037     int * whichRow = NULL;
    1038     int * lookupRow = NULL;
    1039     badModel = (canDo == 0);
    1040     if (numberKnapsack && canDo) {
    1041         /* double check - OK if
    1042            no nonlinear
    1043            nonlinear only on columns in knapsack
    1044            nonlinear only on columns in knapsack * ONE other - same for all in knapsack
    1045            AND that is only row connected to knapsack
    1046            (theoretically could split knapsack if two other and small numbers)
    1047            also ONE could be ONE expression - not just a variable
    1048         */
    1049         int iKnapsack;
    1050         markKnapsack = new int [numberKnapsack];
    1051         coefficient = new double [numberKnapsack];
    1052         linear = new double [numberColumns];
    1053         for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++)
    1054             markKnapsack[iKnapsack] = -1;
    1055         if (canDo == 2) {
    1056             for (iRow = -1; iRow < numberRows; iRow++) {
    1057                 int numberOdd;
    1058                 CoinPackedMatrix * row = coinModel.quadraticRow(iRow, linear, numberOdd);
    1059                 if (row) {
    1060                     // see if valid
    1061                     const double * element = row->getElements();
    1062                     const int * column = row->getIndices();
    1063                     const CoinBigIndex * columnStart = row->getVectorStarts();
    1064                     const int * columnLength = row->getVectorLengths();
    1065                     int numberLook = row->getNumCols();
    1066                     for (int i = 0; i < numberLook; i++) {
    1067                         int iKnapsack = whichKnapsack[i];
    1068                         if (iKnapsack < 0) {
    1069                             // might be able to swap - but for now can't have knapsack in
    1070                             for (int j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    1071                                 int iColumn = column[j];
    1072                                 if (whichKnapsack[iColumn] >= 0) {
    1073                                     canDo = 0; // no good
    1074                                     badModel = true;
    1075                                     break;
    1076                                 }
    1077                             }
    1078                         } else {
    1079                             // OK if in same knapsack - or maybe just one
    1080                             int marked = markKnapsack[iKnapsack];
    1081                             for (int j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    1082                                 int iColumn = column[j];
    1083                                 if (whichKnapsack[iColumn] != iKnapsack && whichKnapsack[iColumn] >= 0) {
    1084                                     canDo = 0; // no good
    1085                                     badModel = true;
    1086                                     break;
    1087                                 } else if (marked == -1) {
    1088                                     markKnapsack[iKnapsack] = iColumn;
    1089                                     marked = iColumn;
    1090                                     coefficient[iKnapsack] = element[j];
    1091                                     coinModel.associateElement(coinModel.columnName(iColumn), 1.0);
    1092                                 } else if (marked != iColumn) {
    1093                                     badModel = true;
    1094                                     canDo = 0; // no good
    1095                                     break;
    1096                                 } else {
    1097                                     // could manage with different coefficients - but for now ...
    1098                                     assert(coefficient[iKnapsack] == element[j]);
    1099                                 }
    1100                             }
    1101                         }
    1102                     }
    1103                     delete row;
    1104                 }
    1105             }
    1106         }
    1107         if (canDo) {
    1108             // for any rows which are cuts
    1109             whichRow = new int [numberRows];
    1110             lookupRow = new int [numberRows];
    1111             bool someNonlinear = false;
    1112             double maxCoefficient = 1.0;
    1113             for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
    1114                 if (markKnapsack[iKnapsack] >= 0) {
    1115                     someNonlinear = true;
    1116                     int iColumn = markKnapsack[iKnapsack];
    1117                     maxCoefficient = CoinMax(maxCoefficient, fabs(coefficient[iKnapsack] * coinModel.columnUpper(iColumn)));
    1118                 }
    1119             }
    1120             if (someNonlinear) {
    1121                 // associate all columns to stop possible error messages
    1122                 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    1123                     coinModel.associateElement(coinModel.columnName(iColumn), 1.0);
    1124                 }
    1125             }
    1126             ClpSimplex tempModel;
    1127             tempModel.loadProblem(coinModel);
    1128             // Create final model - first without knapsacks
    1129             int nCol = 0;
    1130             int nRow = 0;
    1131             for (iRow = 0; iRow < numberRows; iRow++) {
    1132                 if (markRow[iRow] < 0) {
    1133                     lookupRow[iRow] = nRow;
    1134                     whichRow[nRow++] = iRow;
    1135                 } else {
    1136                     lookupRow[iRow] = -1;
    1137                 }
    1138             }
    1139             for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    1140                 if (whichKnapsack[iColumn] < 0)
    1141                     whichColumn[nCol++] = iColumn;
    1142             }
    1143             ClpSimplex finalModelX(&tempModel, nRow, whichRow, nCol, whichColumn, false, false, false);
    1144             OsiClpSolverInterface finalModelY(&finalModelX, true);
    1145             finalModel = finalModelY.clone();
    1146             finalModelY.releaseClp();
    1147             // Put back priorities
    1148             const int * priorities = model.priorities();
    1149             if (priorities) {
    1150                 finalModel->findIntegers(false);
    1151                 OsiObject ** objects = finalModel->objects();
    1152                 int numberObjects = finalModel->numberObjects();
    1153                 for (int iObj = 0; iObj < numberObjects; iObj++) {
    1154                     int iColumn = objects[iObj]->columnNumber();
    1155                     if (iColumn >= 0 && iColumn < nCol) {
    1156 #ifndef NDEBUG
    1157                         OsiSimpleInteger * obj =
    1158                             dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
    1159 #endif
    1160                         assert (obj);
    1161                         int iPriority = priorities[whichColumn[iColumn]];
    1162                         if (iPriority > 0)
    1163                             objects[iObj]->setPriority(iPriority);
    1164                     }
    1165                 }
    1166             }
    1167             for (iRow = 0; iRow < numberRows; iRow++) {
    1168                 whichRow[iRow] = iRow;
    1169             }
    1170             int numberOther = finalModel->getNumCols();
    1171             int nLargest = 0;
    1172             int nelLargest = 0;
    1173             int nTotal = 0;
    1174             for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
    1175                 iRow = knapsackRow[iKnapsack];
    1176                 int nCreate = maxTotal;
    1177                 int nelCreate = coinModel.expandKnapsack(iRow, nCreate, NULL, NULL, NULL, NULL);
    1178                 if (nelCreate < 0)
    1179                     badModel = true;
    1180                 nTotal += nCreate;
    1181                 nLargest = CoinMax(nLargest, nCreate);
    1182                 nelLargest = CoinMax(nelLargest, nelCreate);
    1183             }
    1184             if (nTotal > maxTotal)
    1185                 badModel = true;
    1186             if (!badModel) {
    1187                 // Now arrays for building
    1188                 nelLargest = CoinMax(nelLargest, nLargest) + 1;
    1189                 double * buildObj = new double [nLargest];
    1190                 double * buildElement = new double [nelLargest];
    1191                 int * buildStart = new int[nLargest+1];
    1192                 int * buildRow = new int[nelLargest];
    1193                 // alow for integers in knapsacks
    1194                 OsiObject ** object = new OsiObject * [numberKnapsack+nTotal];
    1195                 int nSOS = 0;
    1196                 int nObj = numberKnapsack;
    1197                 for (iKnapsack = 0; iKnapsack < numberKnapsack; iKnapsack++) {
    1198                     knapsackStart[iKnapsack] = finalModel->getNumCols();
    1199                     iRow = knapsackRow[iKnapsack];
    1200                     int nCreate = 10000;
    1201                     coinModel.expandKnapsack(iRow, nCreate, buildObj, buildStart, buildRow, buildElement);
    1202                     // Redo row numbers
    1203                     for (iColumn = 0; iColumn < nCreate; iColumn++) {
    1204                         for (int j = buildStart[iColumn]; j < buildStart[iColumn+1]; j++) {
    1205                             int jRow = buildRow[j];
    1206                             jRow = lookupRow[jRow];
    1207                             assert (jRow >= 0 && jRow < nRow);
    1208                             buildRow[j] = jRow;
    1209                         }
    1210                     }
    1211                     finalModel->addCols(nCreate, buildStart, buildRow, buildElement, NULL, NULL, buildObj);
    1212                     int numberFinal = finalModel->getNumCols();
    1213                     for (iColumn = numberOther; iColumn < numberFinal; iColumn++) {
    1214                         if (markKnapsack[iKnapsack] < 0) {
    1215                             finalModel->setColUpper(iColumn, maxCoefficient);
    1216                             finalModel->setInteger(iColumn);
    1217                         } else {
    1218                             finalModel->setColUpper(iColumn, maxCoefficient + 1.0);
    1219                             finalModel->setInteger(iColumn);
    1220                         }
    1221                         OsiSimpleInteger * sosObject = new OsiSimpleInteger(finalModel, iColumn);
    1222                         sosObject->setPriority(1000000);
    1223                         object[nObj++] = sosObject;
    1224                         buildRow[iColumn-numberOther] = iColumn;
    1225                         buildElement[iColumn-numberOther] = 1.0;
    1226                     }
    1227                     if (markKnapsack[iKnapsack] < 0) {
    1228                         // convexity row
    1229                         finalModel->addRow(numberFinal - numberOther, buildRow, buildElement, 1.0, 1.0);
    1230                     } else {
    1231                         int iColumn = markKnapsack[iKnapsack];
    1232                         int n = numberFinal - numberOther;
    1233                         buildRow[n] = iColumn;
    1234                         buildElement[n++] = -fabs(coefficient[iKnapsack]);
    1235                         // convexity row (sort of)
    1236                         finalModel->addRow(n, buildRow, buildElement, 0.0, 0.0);
    1237                         OsiSOS * sosObject = new OsiSOS(finalModel, n - 1, buildRow, NULL, 1);
    1238                         sosObject->setPriority(iKnapsack + SOSPriority);
    1239                         // Say not integral even if is (switch off heuristics)
    1240                         sosObject->setIntegerValued(false);
    1241                         object[nSOS++] = sosObject;
    1242                     }
    1243                     numberOther = numberFinal;
    1244                 }
    1245                 finalModel->addObjects(nObj, object);
    1246                 for (iKnapsack = 0; iKnapsack < nObj; iKnapsack++)
    1247                     delete object[iKnapsack];
    1248                 delete [] object;
    1249                 // Can we move any rows to cuts
    1250                 const int * cutMarker = coinModel.cutMarker();
    1251                 if (cutMarker && 0) {
    1252                     printf("AMPL CUTS OFF until global cuts fixed\n");
    1253                     cutMarker = NULL;
    1254                 }
    1255                 if (cutMarker) {
    1256                     // Row copy
    1257                     const CoinPackedMatrix * matrixByRow = finalModel->getMatrixByRow();
    1258                     const double * elementByRow = matrixByRow->getElements();
    1259                     const int * column = matrixByRow->getIndices();
    1260                     const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
    1261                     const int * rowLength = matrixByRow->getVectorLengths();
    1262 
    1263                     const double * rowLower = finalModel->getRowLower();
    1264                     const double * rowUpper = finalModel->getRowUpper();
    1265                     int nDelete = 0;
    1266                     for (iRow = 0; iRow < numberRows; iRow++) {
    1267                         if (cutMarker[iRow] && lookupRow[iRow] >= 0) {
    1268                             int jRow = lookupRow[iRow];
    1269                             whichRow[nDelete++] = jRow;
    1270                             int start = rowStart[jRow];
    1271                             stored.addCut(rowLower[jRow], rowUpper[jRow],
    1272                                           rowLength[jRow], column + start, elementByRow + start);
    1273                         }
    1274                     }
    1275                     finalModel->deleteRows(nDelete, whichRow);
    1276                 }
    1277                 knapsackStart[numberKnapsack] = finalModel->getNumCols();
    1278                 delete [] buildObj;
    1279                 delete [] buildElement;
    1280                 delete [] buildStart;
    1281                 delete [] buildRow;
    1282                 finalModel->writeMps("full");
    1283             }
    1284         }
    1285     }
    1286     delete [] whichKnapsack;
    1287     delete [] markRow;
    1288     delete [] markKnapsack;
    1289     delete [] coefficient;
    1290     delete [] linear;
    1291     delete [] whichRow;
    1292     delete [] lookupRow;
    1293     delete si;
    1294     si = NULL;
    1295     if (!badModel && finalModel) {
    1296         finalModel->setDblParam(OsiObjOffset, coinModel.objectiveOffset());
    1297         return finalModel;
    1298     } else {
    1299         delete finalModel;
    1300         printf("can't make knapsacks - did you set fixedPriority (extra1)\n");
    1301         return NULL;
    1302     }
    1303 }
    1304 #endif  //COIN_HAS_LINK
    1305 
    1306800
    1307801
Note: See TracChangeset for help on using the changeset viewer.