Changeset 1591 for trunk/Cbc/src/CbcHeuristicGreedy.cpp
 Timestamp:
 Feb 7, 2011 7:48:46 AM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Cbc/src/CbcHeuristicGreedy.cpp
r1589 r1591 990 990 const int * columnLength = matrix_.getVectorLengths(); 991 991 int * sosRow = new int [numberColumns]; 992 int nonSOS=0; 992 993 // If bit set then use current 993 994 if ((algorithm_&1)!=0) { … … 1005 1006 // SOS 1006 1007 rhs[iRow]=1.0; 1007 } else if (rowLower[iRow] > 0.0 ) {1008 } else if (rowLower[iRow] > 0.0 && rowUpper[iRow] < 1.0e10) { 1008 1009 good = false; 1009 1010 } else if (rowUpper[iRow] < 0.0) { 1010 1011 good = false; 1012 } else if (rowUpper[iRow] < 1.0e10) { 1013 rhs[iRow]=rowUpper[iRow]; 1011 1014 } else { 1012 rhs[iRow]=row Upper[iRow];1015 rhs[iRow]=rowLower[iRow]; 1013 1016 } 1014 1017 } 1015 1018 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 1019 if (!columnLength[iColumn]) 1020 continue; 1016 1021 if (columnLower[iColumn] < 0.0  columnUpper[iColumn] > 1.0) 1017 1022 good = false; … … 1019 1024 int nSOS=0; 1020 1025 int iSOS=1; 1026 if (!solver>isInteger(iColumn)) 1027 good = false; 1021 1028 for (j = columnStart[iColumn]; 1022 1029 j < columnStart[iColumn] + columnLength[iColumn]; j++) { … … 1031 1038 } 1032 1039 } 1033 if (nSOS>1 !solver>isBinary(iColumn))1040 if (nSOS>1) 1034 1041 good = false; 1042 else if (!nSOS) 1043 nonSOS++; 1035 1044 sosRow[iColumn] = iSOS; 1036 1045 } 1037 1046 if (!good) { 1047 delete [] sosRow; 1038 1048 delete [] rhs; 1039 delete [] sosRow;1040 1049 setWhen(0); // switch off 1041 1050 return 0; … … 1046 1055 const double * solution = solver>getColSolution(); 1047 1056 const double * objective = solver>getObjCoefficients(); 1057 const double * rowLower = solver>getRowLower(); 1058 const double * rowUpper = solver>getRowUpper(); 1048 1059 double integerTolerance = model_>getDblParam(CbcModel::CbcIntegerTolerance); 1049 1060 double primalTolerance; … … 1052 1063 numRuns_++; 1053 1064 assert (numberRows == matrix_.getNumRows()); 1065 // set up linked list for sets 1066 int * firstGub = new int [numberRows]; 1067 int * nextGub = new int [numberColumns]; 1054 1068 int iRow, iColumn; 1055 1069 double direction = solver>getObjSense(); 1056 1070 double * slackCost = new double [numberRows]; 1057 1071 double * modifiedCost = CoinCopyOfArray(objective,numberColumns); 1058 for (int iRow = 0;iRow < numberRows; iRow++) 1072 for (int iRow = 0;iRow < numberRows; iRow++) { 1059 1073 slackCost[iRow]=1.0e30; 1074 firstGub[iRow]=1; 1075 } 1060 1076 // Take off cost of gub slack 1061 1077 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 1078 nextGub[iColumn]=1; 1062 1079 int iRow = sosRow[iColumn]; 1063 1080 if (columnLength[iColumn] == 1&&iRow>=0) { … … 1075 1092 sos[iRow]=1; 1076 1093 rhs[iRow]=1.0; 1094 } else if (rhs[iRow] != rowUpper[iRow]) { 1095 // G row 1096 sos[iRow]=1; 1077 1097 } 1078 1098 if( slackCost[iRow] == 1.0e30) { … … 1089 1109 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1090 1110 int iRow = row[j]; 1091 if (sos[iRow] ) {1111 if (sos[iRow]>0) { 1092 1112 cost = slackCost[iRow]; 1113 if (firstGub[iRow]<0) { 1114 firstGub[iRow]=iColumn; 1115 } else { 1116 int jColumn = firstGub[iRow]; 1117 while (nextGub[jColumn]>=0) 1118 jColumn=nextGub[jColumn]; 1119 nextGub[jColumn]=iColumn; 1120 } 1093 1121 } 1094 1122 } … … 1105 1133 double * newSolution = new double [numberColumns]; 1106 1134 double * rowActivity = new double[numberRows]; 1107 memset(rowActivity, 0, numberRows*sizeof(double)); 1135 double * contribution = new double [numberColumns]; 1136 int * which = new int [numberColumns]; 1137 double * newSolution0 = new double [numberColumns]; 1108 1138 if ((algorithm_&(24))==0) { 1109 1139 // get solution as small as possible 1110 1140 for (iColumn = 0; iColumn < numberColumns; iColumn++) 1111 newSolution [iColumn] = columnLower[iColumn];1141 newSolution0[iColumn] = columnLower[iColumn]; 1112 1142 } else { 1113 1143 // Get rounded down solution 1114 1144 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1115 1145 double value = solution[iColumn]; 1116 1146 // Round down integer 1117 1147 if (fabs(floor(value + 0.5)  value) < integerTolerance) { … … 1120 1150 value = CoinMax(floor(value), columnLower[iColumn]); 1121 1151 } 1122 1123 1124 1125 newSolution[iColumn] = value;1152 // make sure clean 1153 value = CoinMin(value, columnUpper[iColumn]); 1154 value = CoinMax(value, columnLower[iColumn]); 1155 newSolution0[iColumn] = value; 1126 1156 } 1127 1157 } 1128 // get row activity 1129 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1130 CoinBigIndex j; 1131 double value = newSolution[iColumn]; 1132 for (j = columnStart[iColumn]; 1133 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1134 int iRow = row[j]; 1135 rowActivity[iRow] += value * element[j]; 1158 double * rowWeight = new double [numberRows]; 1159 for (int i=0;i<numberRows;i++) 1160 rowWeight[i]=1.0; 1161 double costBias = 0.0; 1162 int nPass = ((algorithm_&4)!=0) ? 1 : 10; 1163 for (int iPass=0;iPass<nPass;iPass++) { 1164 newSolutionValue = offset+offset2; 1165 memcpy(newSolution,newSolution0,numberColumns*sizeof(double)); 1166 // get row activity 1167 memset(rowActivity, 0, numberRows*sizeof(double)); 1168 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1169 CoinBigIndex j; 1170 double value = newSolution[iColumn]; 1171 for (j = columnStart[iColumn]; 1172 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1173 int iRow = row[j]; 1174 rowActivity[iRow] += value * element[j]; 1175 } 1136 1176 } 1137 } 1138 double * contribution = new double [numberColumns]; 1139 int * which = new int [numberColumns]; 1140 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1141 CoinBigIndex j; 1142 double value = newSolution[iColumn]; 1143 double cost = modifiedCost[iColumn]; 1144 double forSort = 0.0; 1145 bool hasSlack=false; 1146 bool willFit=true; 1147 newSolutionValue += value * cost; 1148 for (j = columnStart[iColumn]; 1149 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1150 int iRow = row[j]; 1151 if (sos[iRow] == 2) 1152 hasSlack = true; 1153 forSort += element[j]; 1154 double gap = rhs[iRow]  rowActivity[iRow]+1.0e8; 1155 if (gap<element[j]) 1156 willFit = false; 1177 if (!rowWeight) { 1178 rowWeight = CoinCopyOfArray(rowActivity,numberRows); 1157 1179 } 1158 bool isSlack = hasSlack && (columnLength[iColumn]==1); 1159 if ((algorithm_&4)!=0) 1160 forSort=1.0; 1161 // Use smallest cost if will fit 1162 if (willFit /*&& hasSlack*/ && 1163 value == 0.0 && columnUpper[iColumn]) { 1164 if (hasSlack) { 1165 if (cost>0.0) { 1166 forSort = 2.0e30; 1167 } else if (cost==0.0) { 1168 if (!isSlack) 1169 forSort = 1.0e29; 1180 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1181 CoinBigIndex j; 1182 double value = newSolution[iColumn]; 1183 double cost = modifiedCost[iColumn]; 1184 double forSort = 1.0e24; 1185 bool hasSlack=false; 1186 bool willFit=true; 1187 bool gRow=false; 1188 newSolutionValue += value * cost; 1189 cost += 1.0e12; 1190 for (j = columnStart[iColumn]; 1191 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1192 int iRow = row[j]; 1193 int type = sos[iRow]; 1194 double gap = rhs[iRow]  rowActivity[iRow]+1.0e8; 1195 switch (type) { 1196 case 1: 1197 // G row 1198 gRow = true; 1199 #if 0 1200 if (rhs[iRow]>rowWeight[iRow](algorithm_&(24))!=0) 1201 forSort += element[j]; 1170 1202 else 1171 forSort = 1.0e28; 1203 forSort += 0.1*element[j]; 1204 #else 1205 forSort += rowWeight[iRow]*element[j]; 1206 #endif 1207 break; 1208 case 0: 1209 // L row 1210 if (gap<element[j]) { 1211 willFit = false; 1212 } else { 1213 forSort += element[j]; 1214 } 1215 break; 1216 case 1: 1217 // SOS without slack 1218 if (gap<element[j]) { 1219 willFit = false; 1220 } 1221 break; 1222 case 2: 1223 // SOS with slack 1224 hasSlack = true; 1225 if (gap<element[j]) { 1226 willFit = false; 1227 } 1228 break; 1229 } 1230 } 1231 bool isSlack = hasSlack && (columnLength[iColumn]==1); 1232 if (forSort<1.0e24) 1233 forSort = 1.0e12; 1234 if ((algorithm_&4)!=0 && forSort > 1.0e24) 1235 forSort=1.0; 1236 // Use smallest cost if will fit 1237 if (willFit && (hasSlackgRow) && 1238 value == 0.0 && columnUpper[iColumn]) { 1239 if (hasSlack && !gRow) { 1240 if (cost>1.0e12) { 1241 forSort = 2.0e30; 1242 } else if (cost==1.0e12) { 1243 if (!isSlack) 1244 forSort = 1.0e29; 1245 else 1246 forSort = 1.0e28; 1247 } else { 1248 forSort = cost/forSort; 1249 } 1172 1250 } else { 1173 forSort = cost/forSort; 1251 if (!gRowtrue) 1252 forSort = (cost+costBias)/forSort; 1253 else 1254 forSort = 1.0e12/forSort; 1174 1255 } 1175 1256 } else { 1176 forSort = cost/forSort; 1257 // put at end 1258 forSort = 1.0e30; 1177 1259 } 1178 } else { 1179 // put at end 1180 forSort = 1.0e30; 1260 which[iColumn]=iColumn; 1261 contribution[iColumn]= forSort; 1181 1262 } 1182 which[iColumn]=iColumn; 1183 contribution[iColumn]= forSort; 1184 } 1185 CoinSort_2(contribution,contribution+numberColumns,which); 1186 // Go through columns 1187 for (int jColumn = 0; jColumn < numberColumns; jColumn++) { 1188 int iColumn = which[jColumn]; 1189 double value = newSolution[iColumn]; 1190 if (value) 1191 continue; 1192 bool possible = true; 1193 CoinBigIndex j; 1194 for (j = columnStart[iColumn]; 1195 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1196 int iRow = row[j]; 1197 if (sos[iRow]&&rowActivity[iRow]) { 1198 possible = false; 1199 } else { 1200 double gap = rhs[iRow]  rowActivity[iRow]+1.0e8; 1201 if (gap<element[j]) 1263 CoinSort_2(contribution,contribution+numberColumns,which); 1264 // Go through columns 1265 int nAdded=0; 1266 int nSlacks=0; 1267 for (int jColumn = 0; jColumn < numberColumns; jColumn++) { 1268 if (contribution[jColumn]>=1.0e30) 1269 break; 1270 int iColumn = which[jColumn]; 1271 double value = newSolution[iColumn]; 1272 if (value) 1273 continue; 1274 bool possible = true; 1275 CoinBigIndex j; 1276 for (j = columnStart[iColumn]; 1277 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1278 int iRow = row[j]; 1279 if (sos[iRow]>0&&rowActivity[iRow]) { 1202 1280 possible = false; 1281 } else { 1282 double gap = rhs[iRow]  rowActivity[iRow]+1.0e8; 1283 if (gap<element[j]&&sos[iRow]>=0) 1284 possible = false; 1285 } 1286 } 1287 if (possible) { 1288 //#define REPORT 1 1289 #ifdef REPORT 1290 if ((nAdded%1000)==0) { 1291 double gap=0.0; 1292 for (int i=0;i<numberRows;i++) { 1293 if (rowUpper[i]>1.0e20) 1294 gap += CoinMax(rowLower[i]rowActivity[i],0.0); 1295 } 1296 if (gap) 1297 printf("after %d added gap %g  %d slacks\n", 1298 nAdded,gap,nSlacks); 1299 } 1300 #endif 1301 nAdded++; 1302 if (columnLength[iColumn]==1) 1303 nSlacks++; 1304 // Increase chosen column 1305 newSolution[iColumn] = 1.0; 1306 double cost = modifiedCost[iColumn]; 1307 newSolutionValue += cost; 1308 for (CoinBigIndex j = columnStart[iColumn]; 1309 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1310 int iRow = row[j]; 1311 rowActivity[iRow] += element[j]; 1312 } 1203 1313 } 1204 1314 } 1205 if (possible) { 1206 // Increase chosen column 1207 newSolution[iColumn] = 1.0; 1208 double cost = modifiedCost[iColumn]; 1209 newSolutionValue += cost; 1210 for (CoinBigIndex j = columnStart[iColumn]; 1211 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1212 int iRow = row[j]; 1213 rowActivity[iRow] += element[j]; 1214 } 1315 #ifdef REPORT 1316 { 1317 double under=0.0; 1318 double over=0.0; 1319 double gap = 0.0; 1320 int nUnder=0; 1321 int nOver=0; 1322 int nGap=0; 1323 for (iRow = 0; iRow < numberRows; iRow++) { 1324 if (rowActivity[iRow] < rowLower[iRow]  10.0*primalTolerance) { 1325 double value = rowLower[iRow]rowActivity[iRow]; 1326 #if REPORT>1 1327 printf("below on %d is %g  activity %g lower %g\n", 1328 iRow,value,rowActivity[iRow],rowLower[iRow]); 1329 #endif 1330 under += value; 1331 nUnder++; 1332 } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) { 1333 double value = rowActivity[iRow]rowUpper[iRow]; 1334 #if REPORT>1 1335 printf("above on %d is %g  activity %g upper %g\n", 1336 iRow,value,rowActivity[iRow],rowUpper[iRow]); 1337 #endif 1338 over += value; 1339 nOver++; 1340 } else { 1341 double value = rowActivity[iRow]rowLower[iRow]; 1342 if (value && value < 1.0e20) { 1343 #if REPORT>1 1344 printf("gap on %d is %g  activity %g lower %g\n", 1345 iRow,value,rowActivity[iRow],rowLower[iRow]); 1346 #endif 1347 gap += value; 1348 nGap++; 1349 } 1350 } 1351 } 1352 printf("final under %g (%d)  over %g (%d)  free %g (%d)  %d added  solvalue %g\n", 1353 under,nUnder,over,nOver,gap,nGap,nAdded,newSolutionValue); 1215 1354 } 1216 } 1355 #endif 1356 double gap = 0.0; 1357 double over = 0.0; 1358 int nL=0; 1359 int nG=0; 1360 int iUnder=1; 1361 int nUnder=0; 1362 for (iRow = 0; iRow < numberRows; iRow++) { 1363 if (rowLower[iRow]<1.0e20) 1364 nL++; 1365 if (rowUpper[iRow]>1.0e20) 1366 nG++; 1367 if (rowActivity[iRow] < rowLower[iRow]  10.0*primalTolerance) { 1368 gap += rowLower[iRow]rowActivity[iRow]; 1369 nUnder++; 1370 iUnder=iRow; 1371 rowWeight[iRow] *= 1.1; 1372 } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) { 1373 gap += rowActivity[iRow]rowUpper[iRow]; 1374 } else { 1375 over += rowActivity[iRow]rowLower[iRow]; 1376 //rowWeight[iRow] *= 0.9; 1377 } 1378 } 1379 if (nG&&!nL) { 1380 // can we fix 1381 // get list of columns which can go down without making 1382 // things much worse 1383 int nPossible=0; 1384 int nEasyDown=0; 1385 int nSlackDown=0; 1386 for (int iColumn=0;iColumn<numberColumns;iColumn++) { 1387 if (newSolution[iColumn]&& 1388 columnUpper[iColumn]>columnLower[iColumn]) { 1389 bool canGoDown=true; 1390 bool under = false; 1391 int iSos=1; 1392 for (CoinBigIndex j = columnStart[iColumn]; 1393 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1394 int iRow = row[j]; 1395 if (sos[iRow]<0) { 1396 double over = rowActivity[iRow]rowLower[iRow]; 1397 if (over>=0.0&&element[j]>over+1.0e12) { 1398 canGoDown=false; 1399 break; 1400 } else if (over<0.0) { 1401 under = true; 1402 } 1403 } else { 1404 iSos=iRow; 1405 } 1406 } 1407 if (canGoDown) { 1408 if (!under) { 1409 if (iSos>=0) { 1410 // find cheapest 1411 double cheapest=modifiedCost[iColumn]; 1412 int iCheapest = 1; 1413 int jColumn = firstGub[iSos]; 1414 assert (jColumn>=0); 1415 while (jColumn>=0) { 1416 if (modifiedCost[jColumn]<cheapest) { 1417 cheapest=modifiedCost[jColumn]; 1418 iCheapest=jColumn; 1419 } 1420 jColumn = nextGub[jColumn]; 1421 } 1422 if (iCheapest>=0) { 1423 // Decrease column 1424 newSolution[iColumn] = 0.0; 1425 newSolutionValue = modifiedCost[iColumn]; 1426 for (CoinBigIndex j = columnStart[iColumn]; 1427 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1428 int iRow = row[j]; 1429 rowActivity[iRow] = element[j]; 1430 } 1431 // Increase chosen column 1432 newSolution[iCheapest] = 1.0; 1433 newSolutionValue += modifiedCost[iCheapest]; 1434 for (CoinBigIndex j = columnStart[iCheapest]; 1435 j < columnStart[iCheapest] + columnLength[iCheapest]; j++) { 1436 int iRow = row[j]; 1437 rowActivity[iRow] += element[j]; 1438 } 1439 nEasyDown++; 1440 if (columnLength[iColumn]>1) { 1441 //printf("%d is easy down\n",iColumn); 1442 } else { 1443 nSlackDown++; 1444 } 1445 } 1446 } else if (modifiedCost[iColumn]>0.0) { 1447 // easy down 1448 // Decrease column 1449 newSolution[iColumn] = 0.0; 1450 newSolutionValue = modifiedCost[iColumn]; 1451 for (CoinBigIndex j = columnStart[iColumn]; 1452 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1453 int iRow = row[j]; 1454 rowActivity[iRow] = element[j]; 1455 } 1456 nEasyDown++; 1457 } 1458 } else { 1459 which[nPossible++]=iColumn; 1460 } 1461 } 1462 } 1463 } 1464 #ifdef REPORT 1465 printf("%d possible down, %d easy down of which %d are slacks\n", 1466 nPossible,nEasyDown,nSlackDown); 1467 #endif 1468 double * needed = new double [numberRows]; 1469 for (int i=0;i<numberRows;i++) { 1470 double value = rowLower[i]  rowActivity[i]; 1471 if (value<1.0e8) 1472 value=0.0; 1473 needed[i]=value; 1474 } 1475 if (gap && /*nUnder==1 &&*/ nonSOS) { 1476 double * weight = new double [numberColumns]; 1477 int * sort = new int [numberColumns]; 1478 // look at ones not in set 1479 int nPossible=0; 1480 for (int iColumn=0;iColumn<numberColumns;iColumn++) { 1481 if (!newSolution[iColumn]&& 1482 columnUpper[iColumn]>columnLower[iColumn]) { 1483 int iSos=1; 1484 double value=0.0; 1485 for (CoinBigIndex j = columnStart[iColumn]; 1486 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1487 int iRow = row[j]; 1488 if (sos[iRow]<0) { 1489 if (needed[iRow]) 1490 value += CoinMin(element[j]/needed[iRow],1.0); 1491 } else { 1492 iSos=iRow; 1493 } 1494 } 1495 if (value && iSos<0) { 1496 weight[nPossible]=value; 1497 sort[nPossible++]=iColumn; 1498 } 1499 } 1500 } 1501 CoinSort_2(weight,weight+nPossible,sort); 1502 for (int i=0;i<nPossible;i++) { 1503 int iColumn = sort[i]; 1504 double helps=0.0; 1505 for (CoinBigIndex j = columnStart[iColumn]; 1506 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1507 int iRow = row[j]; 1508 if (needed[iRow]) 1509 helps += CoinMin(needed[iRow],element[j]); 1510 } 1511 if (helps) { 1512 newSolution[iColumn] = 1.0; 1513 newSolutionValue += modifiedCost[iColumn]; 1514 for (CoinBigIndex j = columnStart[iColumn]; 1515 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1516 int iRow = row[j]; 1517 rowActivity[iRow] += element[j]; 1518 if (needed[iRow]) { 1519 needed[iRow] = element[j]; 1520 if (needed[iRow]<1.0e8) 1521 needed[iRow]=0.0; 1522 } 1523 } 1524 gap = helps; 1525 #ifdef REPORT 1526 { 1527 double gap2 = 0.0; 1528 for (iRow = 0; iRow < numberRows; iRow++) { 1529 if (rowActivity[iRow] < rowLower[iRow]  10.0*primalTolerance) { 1530 gap2 += rowLower[iRow]rowActivity[iRow]; 1531 } 1532 } 1533 printf("estimated gap (nonsos) %g  computed %g\n", 1534 gap,gap2); 1535 } 1536 #endif 1537 if (gap<1.0e12) 1538 break; 1539 } 1540 } 1541 delete [] weight; 1542 delete [] sort; 1543 } 1544 if (gap&&nPossible/*&&nUnder==1*/&&true&&model_>bestSolution()) { 1545 double * weight = new double [numberColumns]; 1546 int * sort = new int [numberColumns]; 1547 // look at ones in sets 1548 const double * goodSolution = model_>bestSolution(); 1549 int nPossible=0; 1550 double largestWeight=0.0; 1551 for (int iColumn=0;iColumn<numberColumns;iColumn++) { 1552 if (!newSolution[iColumn]&&goodSolution[iColumn]&& 1553 columnUpper[iColumn]>columnLower[iColumn]) { 1554 int iSos=1; 1555 double value=0.0; 1556 for (CoinBigIndex j = columnStart[iColumn]; 1557 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1558 int iRow = row[j]; 1559 if (sos[iRow]<0) { 1560 if (needed[iRow]) 1561 value += CoinMin(element[j]/needed[iRow],1.0); 1562 } else { 1563 iSos=iRow; 1564 } 1565 } 1566 if (value&&iSos>=0) { 1567 // see if value bigger than current 1568 int jColumn = firstGub[iSos]; 1569 assert (jColumn>=0); 1570 while (jColumn>=0) { 1571 if (newSolution[jColumn]) 1572 break; 1573 jColumn = nextGub[jColumn]; 1574 } 1575 assert (jColumn>=0); 1576 double value2=0.0; 1577 for (CoinBigIndex j = columnStart[jColumn]; 1578 j < columnStart[jColumn] + columnLength[jColumn]; j++) { 1579 int iRow = row[j]; 1580 if (needed[iRow]) 1581 value2 += CoinMin(element[j]/needed[iRow],1.0); 1582 } 1583 if (value>value2) { 1584 weight[nPossible]=(valuevalue2); 1585 largestWeight = CoinMax(largestWeight,(valuevalue2)); 1586 sort[nPossible++]=iColumn; 1587 } 1588 } 1589 } 1590 } 1591 if (nPossible) { 1592 double * temp = new double [numberRows]; 1593 int * which2 = new int [numberRows]; 1594 memset(temp,0,numberRows*sizeof(double)); 1595 // modify so ones just more than gap best 1596 if (largestWeight>gap&&nUnder==1) { 1597 double offset = 4*largestWeight; 1598 for (int i=0;i<nPossible;i++) { 1599 double value = weight[i]; 1600 if (value>gap1.0e12) 1601 weight[i] = (offset(valuegap)); 1602 } 1603 } 1604 CoinSort_2(weight,weight+nPossible,sort); 1605 for (int i=0;i<nPossible;i++) { 1606 int iColumn = sort[i]; 1607 int n=0; 1608 // find jColumn 1609 int iSos=1; 1610 for (CoinBigIndex j = columnStart[iColumn]; 1611 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 1612 int iRow = row[j]; 1613 temp[iRow]=element[j]; 1614 which2[n++]=iRow; 1615 if (sos[iRow]>=0) { 1616 iSos=iRow; 1617 } 1618 } 1619 int jColumn = firstGub[iSos]; 1620 assert (jColumn>=0); 1621 while (jColumn>=0) { 1622 if (newSolution[jColumn]) 1623 break; 1624 jColumn = nextGub[jColumn]; 1625 } 1626 assert (jColumn>=0); 1627 for (CoinBigIndex j = columnStart[jColumn]; 1628 j < columnStart[jColumn] + columnLength[jColumn]; j++) { 1629 int iRow = row[j]; 1630 if (!temp[iRow]) 1631 which2[n++]=iRow; 1632 temp[iRow] = element[j]; 1633 } 1634 double helps = 0.0; 1635 for (int i=0;i<n;i++) { 1636 int iRow = which2[i]; 1637 double newValue = rowActivity[iRow]+temp[iRow]; 1638 if (temp[iRow]>1.0e8) { 1639 if (rowActivity[iRow]<rowLower[iRow]1.0e8) { 1640 helps += CoinMin(temp[iRow], 1641 rowLower[iRow]rowActivity[iRow]); 1642 } 1643 } else if (temp[iRow]<1.0e8) { 1644 if (newValue<rowLower[iRow]1.0e12) { 1645 helps = CoinMin(temp[iRow], 1646 1.0*(rowLower[iRow]newValue)); 1647 } 1648 } 1649 } 1650 if (helps>0.0) { 1651 newSolution[iColumn]=1.0; 1652 newSolution[jColumn]=0.0; 1653 newSolutionValue += modifiedCost[iColumn]modifiedCost[jColumn]; 1654 for (int i=0;i<n;i++) { 1655 int iRow = which2[i]; 1656 double newValue = rowActivity[iRow]+temp[iRow]; 1657 rowActivity[iRow] = newValue; 1658 temp[iRow]=0.0; 1659 } 1660 gap = helps; 1661 #ifdef REPORT 1662 { 1663 double gap2 = 0.0; 1664 for (iRow = 0; iRow < numberRows; iRow++) { 1665 if (rowActivity[iRow] < rowLower[iRow]  10.0*primalTolerance) { 1666 gap2 += rowLower[iRow]rowActivity[iRow]; 1667 } 1668 } 1669 printf("estimated gap %g  computed %g\n", 1670 gap,gap2); 1671 } 1672 #endif 1673 if (gap<1.0e8) 1674 break; 1675 } else { 1676 for (int i=0;i<n;i++) 1677 temp[which2[i]]=0.0; 1678 } 1679 } 1680 delete [] which2; 1681 delete [] temp; 1682 } 1683 delete [] weight; 1684 delete [] sort; 1685 } 1686 delete [] needed; 1687 } 1688 #ifdef REPORT 1689 { 1690 double gap=0.0; 1691 double over = 0.0; 1692 for (iRow = 0; iRow < numberRows; iRow++) { 1693 if (rowActivity[iRow] < rowLower[iRow]  10.0*primalTolerance) { 1694 double value = rowLower[iRow]rowActivity[iRow]; 1695 #if REPORT>1 1696 printf("below on %d is %g  activity %g lower %g\n", 1697 iRow,value,rowActivity[iRow],rowLower[iRow]); 1698 #endif 1699 gap += value; 1700 } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) { 1701 double value = rowActivity[iRow]rowUpper[iRow]; 1702 #if REPORT>1 1703 printf("above on %d is %g  activity %g upper %g\n", 1704 iRow,value,rowActivity[iRow],rowUpper[iRow]); 1705 #endif 1706 gap += value; 1707 } else { 1708 double value = rowActivity[iRow]rowLower[iRow]; 1709 if (value) { 1710 #if REPORT>1 1711 printf("over on %d is %g  activity %g lower %g\n", 1712 iRow,value,rowActivity[iRow],rowLower[iRow]); 1713 #endif 1714 over += value; 1715 } 1716 } 1717 } 1718 printf("modified final gap %g  over %g  %d added  solvalue %g\n", 1719 gap,over,nAdded,newSolutionValue); 1720 } 1721 #endif 1722 if (!gap) { 1723 break; 1724 } else { 1725 if (iPass==0) { 1726 costBias = 10.0*newSolutionValue/static_cast<double>(nAdded); 1727 } else { 1728 costBias *= 10.0; 1729 } 1730 } 1731 } 1732 delete [] newSolution0; 1733 delete [] rowWeight; 1217 1734 delete [] sos; 1735 delete [] firstGub; 1736 delete [] nextGub; 1218 1737 if (newSolutionValue < solutionValue) { 1219 1738 // check feasible 1220 const double * rowLower = solver>getRowLower();1221 const double * rowUpper = solver>getRowUpper();1222 1739 memset(rowActivity, 0, numberRows*sizeof(double)); 1223 1740 for (iColumn = 0; iColumn < numberColumns; iColumn++) { … … 1293 1810 setWhen(0); 1294 1811 } 1295 // Only works if coefficients positive and all rows L or SOS1812 // Only works if coefficients positive and all rows L/G or SOS 1296 1813 OsiSolverInterface * solver = model_>solver(); 1297 1814 const double * columnUpper = solver>getColUpper(); … … 1312 1829 // SOS 1313 1830 originalRhs_[iRow]=1.0; 1314 } else if (rowLower[iRow] > 0.0 ) {1831 } else if (rowLower[iRow] > 0.0 && rowUpper[iRow] < 1.0e10) { 1315 1832 good = false; 1316 1833 } else if (rowUpper[iRow] < 0.0) { 1317 1834 good = false; 1835 } else if (rowUpper[iRow] < 1.0e10) { 1836 originalRhs_[iRow]=rowUpper[iRow]; 1318 1837 } else { 1319 originalRhs_[iRow]=row Upper[iRow];1838 originalRhs_[iRow]=rowLower[iRow]; 1320 1839 } 1321 1840 } 1322 1841 int numberColumns = solver>getNumCols(); 1323 1842 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 1843 if (!columnLength[iColumn]) 1844 continue; 1324 1845 if (columnLower[iColumn] < 0.0  columnUpper[iColumn] > 1.0) 1325 1846 good = false; 1326 1847 CoinBigIndex j; 1327 1848 int nSOS=0; 1849 if (!solver>isInteger(iColumn)) 1850 good = false; 1328 1851 for (j = columnStart[iColumn]; 1329 1852 j < columnStart[iColumn] + columnLength[iColumn]; j++) { … … 1337 1860 } 1338 1861 } 1339 if (nSOS >1!solver>isBinary(iColumn))1862 if (nSOS > 1) 1340 1863 good = false; 1341 1864 }
Note: See TracChangeset
for help on using the changeset viewer.