Changeset 1398 for branches/sandbox/Cbc/src
 Timestamp:
 Dec 10, 2009 11:46:08 PM (10 years ago)
 Location:
 branches/sandbox/Cbc/src
 Files:

 2 added
 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/sandbox/Cbc/src/CbcSolver.cpp
r1395 r1398 202 202 203 203 #include "CbcSolverAnalyze.hpp" 204 #include "CbcSolverExpandKnapsack.hpp" 204 205 205 206 #include "CbcSolver.hpp" … … 797 798 } 798 799 #endif 799 800 801 802 #ifdef COIN_HAS_LINK803 /* Returns OsiSolverInterface (User should delete)804 On entry numberKnapsack is maximum number of Total entries805 806 Expanding possibilities of x*y, where x*y are both integers, constructing807 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 model817 OsiSolverLink *si = new OsiSolverLink();818 OsiSolverInterface * finalModel = NULL;819 si>setDefaultMeshSize(0.001);820 // need some relative granularity821 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 priorities828 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 NDEBUG837 OsiSimpleInteger * obj =838 dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;839 #endif840 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 variables857 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 // analyze864 if (logLevel > 1) {865 for (iRow = 0; iRow < numberRows; iRow++) {866 /* Just obvious one at first867 positive non unit coefficients868 all integer869 positive rowUpper870 for now  linear (but further down in code may use nonlinear)871 column bounds should be tight872 */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; // fixed887 }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 linear904 }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 // print919 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 first936 positive non unit coefficients937 all integer938 positive rowUpper939 for now  linear (but further down in code may use nonlinear)940 column bounds should be tight941 */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; // fixed956 }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 linear973 }974 triple = coinModel.next(triple);975 }976 if (n  n1 > 1 && possible) {977 // try978 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 nonlinearities992 /* mark rows993 2 in knapsack and other variables994 1 not involved995 n only in knapsack n996 */997 int * markRow = new int [numberRows];998 for (iRow = 0; iRow < numberRows; iRow++)999 markRow[iRow] = 1;1000 int canDo = 1; // OK and linear1001 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 objective1006 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 if1042 no nonlinear1043 nonlinear only on columns in knapsack1044 nonlinear only on columns in knapsack * ONE other  same for all in knapsack1045 AND that is only row connected to knapsack1046 (theoretically could split knapsack if two other and small numbers)1047 also ONE could be ONE expression  not just a variable1048 */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 valid1061 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 in1070 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 good1074 badModel = true;1075 break;1076 }1077 }1078 } else {1079 // OK if in same knapsack  or maybe just one1080 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 good1085 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 good1095 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 cuts1109 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 messages1122 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 knapsacks1129 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 priorities1148 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 NDEBUG1157 OsiSimpleInteger * obj =1158 dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;1159 #endif1160 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 building1188 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 knapsacks1194 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 numbers1203 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[iColumnnumberOther] = iColumn;1225 buildElement[iColumnnumberOther] = 1.0;1226 }1227 if (markKnapsack[iKnapsack] < 0) {1228 // convexity row1229 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 cuts1250 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 copy1257 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_LINK1305 1306 800 1307 801
Note: See TracChangeset
for help on using the changeset viewer.