Changeset 1383 for branches/sandbox
 Timestamp:
 Dec 9, 2009 11:04:28 AM (10 years ago)
 Location:
 branches/sandbox/Cbc/src
 Files:

 2 added
 2 edited
Legend:
 Unmodified
 Added
 Removed

branches/sandbox/Cbc/src/CbcSolver.cpp
r1381 r1383 45 45 #include "OsiChooseVariable.hpp" 46 46 #include "OsiAuxInfo.hpp" 47 48 #include "CbcSolverHeuristics.hpp" 49 47 50 // Version 48 51 #ifndef CBCVERSION … … 1206 1209 } 1207 1210 #endif 1208 #if 11209 #include "ClpSimplexOther.hpp"1210 1211 1211 // Crunch down model1212 static void1213 crunchIt(ClpSimplex * model)1214 {1215 #if 01216 model>dual();1217 #else1218 int numberColumns = model>numberColumns();1219 int numberRows = model>numberRows();1220 // Use dual region1221 double * rhs = model>dualRowSolution();1222 int * whichRow = new int[3*numberRows];1223 int * whichColumn = new int[2*numberColumns];1224 int nBound;1225 ClpSimplex * small = static_cast<ClpSimplexOther *> (model)>crunch(rhs, whichRow, whichColumn,1226 nBound, false, false);1227 if (small) {1228 small>dual();1229 if (small>problemStatus() == 0) {1230 model>setProblemStatus(0);1231 static_cast<ClpSimplexOther *> (model)>afterCrunch(*small, whichRow, whichColumn, nBound);1232 } else if (small>problemStatus() != 3) {1233 model>setProblemStatus(1);1234 } else {1235 if (small>problemStatus() == 3) {1236 // may be problems1237 small>computeObjectiveValue();1238 model>setObjectiveValue(small>objectiveValue());1239 model>setProblemStatus(3);1240 } else {1241 model>setProblemStatus(3);1242 }1243 }1244 delete small;1245 } else {1246 model>setProblemStatus(1);1247 }1248 delete [] whichRow;1249 delete [] whichColumn;1250 #endif1251 }1252 /*1253 On input1254 doAction  0 just fix in original and return NULL1255 1 return fixed nonpresolved solver1256 2 as one but use presolve Inside this1257 3 use presolve and fix ones with large cost1258 ? do heuristics and set best solution1259 ? do BAB and just set best solution1260 10+ then use lastSolution and relax a few1261 2 cleanup afterwards if using 21262 On output  number fixed1263 */1264 static OsiClpSolverInterface *1265 fixVubs(CbcModel & model, int skipZero2,1266 int & doAction,1267 CoinMessageHandler * /*generalMessageHandler*/,1268 const double * lastSolution, double dextra[6],1269 int extra[5])1270 {1271 if (doAction == 11 && !lastSolution)1272 lastSolution = model.bestSolution();1273 assert (((doAction >= 0 && doAction <= 3) && !lastSolution)  (doAction == 11 && lastSolution));1274 double fractionIntFixed = dextra[3];1275 double fractionFixed = dextra[4];1276 double fixAbove = dextra[2];1277 double fixAboveValue = (dextra[5] > 0.0) ? dextra[5] : 1.0;1278 double time1 = CoinCpuTime();1279 int leaveIntFree = extra[1];1280 OsiSolverInterface * originalSolver = model.solver();1281 OsiClpSolverInterface * originalClpSolver = dynamic_cast< OsiClpSolverInterface*> (originalSolver);1282 ClpSimplex * originalLpSolver = originalClpSolver>getModelPtr();1283 int * originalColumns = NULL;1284 OsiClpSolverInterface * clpSolver;1285 ClpSimplex * lpSolver;1286 ClpPresolve pinfo;1287 assert(originalSolver>getObjSense() > 0);1288 if (doAction == 2  doAction == 3) {1289 double * saveLB = NULL;1290 double * saveUB = NULL;1291 int numberColumns = originalLpSolver>numberColumns();1292 if (fixAbove > 0.0) {1293 double time1 = CoinCpuTime();1294 originalClpSolver>initialSolve();1295 printf("first solve took %g seconds\n", CoinCpuTime()  time1);1296 double * columnLower = originalLpSolver>columnLower() ;1297 double * columnUpper = originalLpSolver>columnUpper() ;1298 const double * solution = originalLpSolver>primalColumnSolution();1299 saveLB = CoinCopyOfArray(columnLower, numberColumns);1300 saveUB = CoinCopyOfArray(columnUpper, numberColumns);1301 const double * objective = originalLpSolver>getObjCoefficients() ;1302 int iColumn;1303 int nFix = 0;1304 int nArt = 0;1305 for (iColumn = 0; iColumn < numberColumns; iColumn++) {1306 if (objective[iColumn] > fixAbove) {1307 if (solution[iColumn] < columnLower[iColumn] + 1.0e8) {1308 columnUpper[iColumn] = columnLower[iColumn];1309 nFix++;1310 } else {1311 nArt++;1312 }1313 } else if (objective[iColumn] < fixAbove) {1314 if (solution[iColumn] > columnUpper[iColumn]  1.0e8) {1315 columnLower[iColumn] = columnUpper[iColumn];1316 nFix++;1317 } else {1318 nArt++;1319 }1320 }1321 }1322 printf("%d artificials fixed, %d left as in solution\n", nFix, nArt);1323 lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e8, true, 10);1324 if (!lpSolver  doAction == 2) {1325 // take off fixing in original1326 memcpy(columnLower, saveLB, numberColumns*sizeof(double));1327 memcpy(columnUpper, saveUB, numberColumns*sizeof(double));1328 }1329 delete [] saveLB;1330 delete [] saveUB;1331 if (!lpSolver) {1332 // try again1333 pinfo.destroyPresolve();1334 lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e8, true, 10);1335 assert (lpSolver);1336 }1337 } else {1338 lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e8, true, 10);1339 assert (lpSolver);1340 }1341 clpSolver = new OsiClpSolverInterface(lpSolver, true);1342 assert(lpSolver == clpSolver>getModelPtr());1343 numberColumns = lpSolver>numberColumns();1344 originalColumns = CoinCopyOfArray(pinfo.originalColumns(), numberColumns);1345 doAction = 1;1346 } else {1347 OsiSolverInterface * solver = originalSolver>clone();1348 clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);1349 lpSolver = clpSolver>getModelPtr();1350 }1351 // Tighten bounds1352 lpSolver>tightenPrimalBounds(0.0, 11, true);1353 int numberColumns = clpSolver>getNumCols() ;1354 double * saveColumnLower = CoinCopyOfArray(lpSolver>columnLower(), numberColumns);1355 double * saveColumnUpper = CoinCopyOfArray(lpSolver>columnUpper(), numberColumns);1356 //char generalPrint[200];1357 const double *objective = lpSolver>getObjCoefficients() ;1358 double *columnLower = lpSolver>columnLower() ;1359 double *columnUpper = lpSolver>columnUpper() ;1360 int numberRows = clpSolver>getNumRows();1361 int iRow, iColumn;1362 1212 1363 // Row copy1364 CoinPackedMatrix matrixByRow(*clpSolver>getMatrixByRow());1365 const double * elementByRow = matrixByRow.getElements();1366 const int * column = matrixByRow.getIndices();1367 const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();1368 const int * rowLength = matrixByRow.getVectorLengths();1369 1370 // Column copy1371 CoinPackedMatrix matrixByCol(*clpSolver>getMatrixByCol());1372 //const double * element = matrixByCol.getElements();1373 const int * row = matrixByCol.getIndices();1374 const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();1375 const int * columnLength = matrixByCol.getVectorLengths();1376 1377 const double * rowLower = clpSolver>getRowLower();1378 const double * rowUpper = clpSolver>getRowUpper();1379 1380 // Get maximum size of VUB tree1381 // otherColumn is one fixed to 0 if this one zero1382 int nEl = matrixByCol.getNumElements();1383 CoinBigIndex * fixColumn = new CoinBigIndex [numberColumns+1];1384 int * otherColumn = new int [nEl];1385 int * fix = new int[numberColumns];1386 char * mark = new char [numberColumns];1387 memset(mark, 0, numberColumns);1388 int numberInteger = 0;1389 int numberOther = 0;1390 fixColumn[0] = 0;1391 double large = lpSolver>largeValue(); // treat bounds > this as infinite1392 #ifndef NDEBUG1393 double large2 = 1.0e10 * large;1394 #endif1395 double tolerance = lpSolver>primalTolerance();1396 int * check = new int[numberRows];1397 for (iRow = 0; iRow < numberRows; iRow++) {1398 check[iRow] = 2; // don't check1399 if (rowLower[iRow] < 1.0e6 && rowUpper[iRow] > 1.0e6)1400 continue;// unlikely1401 // possible row1402 int numberPositive = 0;1403 int iPositive = 1;1404 int numberNegative = 0;1405 int iNegative = 1;1406 CoinBigIndex rStart = rowStart[iRow];1407 CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];1408 CoinBigIndex j;1409 int kColumn;1410 for (j = rStart; j < rEnd; ++j) {1411 double value = elementByRow[j];1412 kColumn = column[j];1413 if (columnUpper[kColumn] > columnLower[kColumn]) {1414 if (value > 0.0) {1415 numberPositive++;1416 iPositive = kColumn;1417 } else {1418 numberNegative++;1419 iNegative = kColumn;1420 }1421 }1422 }1423 if (numberPositive == 1 && numberNegative == 1)1424 check[iRow] = 1; // try both1425 if (numberPositive == 1 && rowLower[iRow] > 1.0e20)1426 check[iRow] = iPositive;1427 else if (numberNegative == 1 && rowUpper[iRow] < 1.0e20)1428 check[iRow] = iNegative;1429 }1430 for (iColumn = 0; iColumn < numberColumns; iColumn++) {1431 fix[iColumn] = 1;1432 if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e8) {1433 if (clpSolver>isInteger(iColumn))1434 numberInteger++;1435 if (columnLower[iColumn] == 0.0) {1436 bool infeasible = false;1437 fix[iColumn] = 0;1438 // fake upper bound1439 double saveUpper = columnUpper[iColumn];1440 columnUpper[iColumn] = 0.0;1441 for (CoinBigIndex i = columnStart[iColumn];1442 i < columnStart[iColumn] + columnLength[iColumn]; i++) {1443 iRow = row[i];1444 if (check[iRow] != 1 && check[iRow] != iColumn)1445 continue; // unlikely1446 // possible row1447 int infiniteUpper = 0;1448 int infiniteLower = 0;1449 double maximumUp = 0.0;1450 double maximumDown = 0.0;1451 double newBound;1452 CoinBigIndex rStart = rowStart[iRow];1453 CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];1454 CoinBigIndex j;1455 int kColumn;1456 // Compute possible lower and upper ranges1457 for (j = rStart; j < rEnd; ++j) {1458 double value = elementByRow[j];1459 kColumn = column[j];1460 if (value > 0.0) {1461 if (columnUpper[kColumn] >= large) {1462 ++infiniteUpper;1463 } else {1464 maximumUp += columnUpper[kColumn] * value;1465 }1466 if (columnLower[kColumn] <= large) {1467 ++infiniteLower;1468 } else {1469 maximumDown += columnLower[kColumn] * value;1470 }1471 } else if (value < 0.0) {1472 if (columnUpper[kColumn] >= large) {1473 ++infiniteLower;1474 } else {1475 maximumDown += columnUpper[kColumn] * value;1476 }1477 if (columnLower[kColumn] <= large) {1478 ++infiniteUpper;1479 } else {1480 maximumUp += columnLower[kColumn] * value;1481 }1482 }1483 }1484 // Build in a margin of error1485 maximumUp += 1.0e8 * fabs(maximumUp);1486 maximumDown = 1.0e8 * fabs(maximumDown);1487 double maxUp = maximumUp + infiniteUpper * 1.0e31;1488 double maxDown = maximumDown  infiniteLower * 1.0e31;1489 if (maxUp <= rowUpper[iRow] + tolerance &&1490 maxDown >= rowLower[iRow]  tolerance) {1491 //printf("Redundant row in vubs %d\n",iRow);1492 } else {1493 if (maxUp < rowLower[iRow]  100.0*tolerance 1494 maxDown > rowUpper[iRow] + 100.0*tolerance) {1495 infeasible = true;1496 break;1497 }1498 double lower = rowLower[iRow];1499 double upper = rowUpper[iRow];1500 for (j = rStart; j < rEnd; ++j) {1501 double value = elementByRow[j];1502 kColumn = column[j];1503 double nowLower = columnLower[kColumn];1504 double nowUpper = columnUpper[kColumn];1505 if (value > 0.0) {1506 // positive value1507 if (lower > large) {1508 if (!infiniteUpper) {1509 assert(nowUpper < large2);1510 newBound = nowUpper +1511 (lower  maximumUp) / value;1512 // relax if original was large1513 if (fabs(maximumUp) > 1.0e8)1514 newBound = 1.0e12 * fabs(maximumUp);1515 } else if (infiniteUpper == 1 && nowUpper > large) {1516 newBound = (lower  maximumUp) / value;1517 // relax if original was large1518 if (fabs(maximumUp) > 1.0e8)1519 newBound = 1.0e12 * fabs(maximumUp);1520 } else {1521 newBound = COIN_DBL_MAX;1522 }1523 if (newBound > nowLower + 1.0e12 && newBound > large) {1524 // Tighten the lower bound1525 // check infeasible (relaxed)1526 if (nowUpper < newBound) {1527 if (nowUpper  newBound <1528 100.0*tolerance) {1529 infeasible = true;1530 break;1531 }1532 }1533 }1534 }1535 if (upper < large) {1536 if (!infiniteLower) {1537 assert(nowLower >  large2);1538 newBound = nowLower +1539 (upper  maximumDown) / value;1540 // relax if original was large1541 if (fabs(maximumDown) > 1.0e8)1542 newBound += 1.0e12 * fabs(maximumDown);1543 } else if (infiniteLower == 1 && nowLower < large) {1544 newBound = (upper  maximumDown) / value;1545 // relax if original was large1546 if (fabs(maximumDown) > 1.0e8)1547 newBound += 1.0e12 * fabs(maximumDown);1548 } else {1549 newBound = COIN_DBL_MAX;1550 }1551 if (newBound < nowUpper  1.0e12 && newBound < large) {1552 // Tighten the upper bound1553 // check infeasible (relaxed)1554 if (nowLower > newBound) {1555 if (newBound  nowLower <1556 100.0*tolerance) {1557 infeasible = true;1558 break;1559 } else {1560 newBound = nowLower;1561 }1562 }1563 if (!newBound  (clpSolver>isInteger(kColumn) && newBound < 0.999)) {1564 // fix to zero1565 if (!mark[kColumn]) {1566 otherColumn[numberOther++] = kColumn;1567 mark[kColumn] = 1;1568 if (check[iRow] == 1)1569 check[iRow] = iColumn;1570 else1571 assert(check[iRow] == iColumn);1572 }1573 }1574 }1575 }1576 } else {1577 // negative value1578 if (lower > large) {1579 if (!infiniteUpper) {1580 assert(nowLower < large2);1581 newBound = nowLower +1582 (lower  maximumUp) / value;1583 // relax if original was large1584 if (fabs(maximumUp) > 1.0e8)1585 newBound += 1.0e12 * fabs(maximumUp);1586 } else if (infiniteUpper == 1 && nowLower < large) {1587 newBound = (lower  maximumUp) / value;1588 // relax if original was large1589 if (fabs(maximumUp) > 1.0e8)1590 newBound += 1.0e12 * fabs(maximumUp);1591 } else {1592 newBound = COIN_DBL_MAX;1593 }1594 if (newBound < nowUpper  1.0e12 && newBound < large) {1595 // Tighten the upper bound1596 // check infeasible (relaxed)1597 if (nowLower > newBound) {1598 if (newBound  nowLower <1599 100.0*tolerance) {1600 infeasible = true;1601 break;1602 } else {1603 newBound = nowLower;1604 }1605 }1606 if (!newBound  (clpSolver>isInteger(kColumn) && newBound < 0.999)) {1607 // fix to zero1608 if (!mark[kColumn]) {1609 otherColumn[numberOther++] = kColumn;1610 mark[kColumn] = 1;1611 if (check[iRow] == 1)1612 check[iRow] = iColumn;1613 else1614 assert(check[iRow] == iColumn);1615 }1616 }1617 }1618 }1619 if (upper < large) {1620 if (!infiniteLower) {1621 assert(nowUpper < large2);1622 newBound = nowUpper +1623 (upper  maximumDown) / value;1624 // relax if original was large1625 if (fabs(maximumDown) > 1.0e8)1626 newBound = 1.0e12 * fabs(maximumDown);1627 } else if (infiniteLower == 1 && nowUpper > large) {1628 newBound = (upper  maximumDown) / value;1629 // relax if original was large1630 if (fabs(maximumDown) > 1.0e8)1631 newBound = 1.0e12 * fabs(maximumDown);1632 } else {1633 newBound = COIN_DBL_MAX;1634 }1635 if (newBound > nowLower + 1.0e12 && newBound > large) {1636 // Tighten the lower bound1637 // check infeasible (relaxed)1638 if (nowUpper < newBound) {1639 if (nowUpper  newBound <1640 100.0*tolerance) {1641 infeasible = true;1642 break;1643 }1644 }1645 }1646 }1647 }1648 }1649 }1650 }1651 for (int i = fixColumn[iColumn]; i < numberOther; i++)1652 mark[otherColumn[i]] = 0;1653 // reset bound unless infeasible1654 if (!infeasible  !clpSolver>isInteger(iColumn))1655 columnUpper[iColumn] = saveUpper;1656 else if (clpSolver>isInteger(iColumn))1657 columnLower[iColumn] = 1.0;1658 }1659 }1660 fixColumn[iColumn+1] = numberOther;1661 }1662 delete [] check;1663 delete [] mark;1664 // Now do reverse way1665 int * counts = new int [numberColumns];1666 CoinZeroN(counts, numberColumns);1667 for (iColumn = 0; iColumn < numberColumns; iColumn++) {1668 for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++)1669 counts[otherColumn[i]]++;1670 }1671 numberOther = 0;1672 CoinBigIndex * fixColumn2 = new CoinBigIndex [numberColumns+1];1673 int * otherColumn2 = new int [fixColumn[numberColumns]];1674 fixColumn2[0] = 0;1675 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {1676 numberOther += counts[iColumn];1677 counts[iColumn] = 0;1678 fixColumn2[iColumn+1] = numberOther;1679 }1680 // Create other way1681 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {1682 for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {1683 int jColumn = otherColumn[i];1684 CoinBigIndex put = fixColumn2[jColumn] + counts[jColumn];1685 counts[jColumn]++;1686 otherColumn2[put] = iColumn;1687 }1688 }1689 // get top layer i.e. those which are not fixed by any other1690 int kLayer = 0;1691 while (true) {1692 int numberLayered = 0;1693 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {1694 if (fix[iColumn] == kLayer) {1695 for (int i = fixColumn2[iColumn]; i < fixColumn2[iColumn+1]; i++) {1696 int jColumn = otherColumn2[i];1697 if (fix[jColumn] == kLayer) {1698 fix[iColumn] = kLayer + 100;1699 }1700 }1701 }1702 if (fix[iColumn] == kLayer) {1703 numberLayered++;1704 }1705 }1706 if (numberLayered) {1707 kLayer += 100;1708 } else {1709 break;1710 }1711 }1712 for (int iPass = 0; iPass < 2; iPass++) {1713 for (int jLayer = 0; jLayer < kLayer; jLayer++) {1714 int check[] = { 1, 0, 1, 2, 3, 4, 5, 10, 50, 100, 500, 1000, 5000, 10000, COIN_INT_MAX};1715 int nCheck = static_cast<int> (sizeof(check) / sizeof(int));1716 int countsI[20];1717 int countsC[20];1718 assert (nCheck <= 20);1719 memset(countsI, 0, nCheck*sizeof(int));1720 memset(countsC, 0, nCheck*sizeof(int));1721 check[nCheck1] = numberColumns;1722 int numberLayered = 0;1723 int numberInteger = 0;1724 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {1725 if (fix[iColumn] == jLayer) {1726 numberLayered++;1727 int nFix = fixColumn[iColumn+1]  fixColumn[iColumn];1728 if (iPass) {1729 // just integers1730 nFix = 0;1731 for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {1732 int jColumn = otherColumn[i];1733 if (clpSolver>isInteger(jColumn))1734 nFix++;1735 }1736 }1737 int iFix;1738 for (iFix = 0; iFix < nCheck; iFix++) {1739 if (nFix <= check[iFix])1740 break;1741 }1742 assert (iFix < nCheck);1743 if (clpSolver>isInteger(iColumn)) {1744 numberInteger++;1745 countsI[iFix]++;1746 } else {1747 countsC[iFix]++;1748 }1749 }1750 }1751 if (numberLayered) {1752 printf("%d (%d integer) at priority %d\n", numberLayered, numberInteger, 1 + (jLayer / 100));1753 char buffer[50];1754 for (int i = 1; i < nCheck; i++) {1755 if (countsI[i]  countsC[i]) {1756 if (i == 1)1757 sprintf(buffer, " == zero ");1758 else if (i < nCheck  1)1759 sprintf(buffer, "> %6d and <= %6d ", check[i1], check[i]);1760 else1761 sprintf(buffer, "> %6d ", check[i1]);1762 printf("%s %8d integers and %8d continuous\n", buffer, countsI[i], countsC[i]);1763 }1764 }1765 }1766 }1767 }1768 delete [] counts;1769 // Now do fixing1770 {1771 // switch off presolve and up weight1772 ClpSolve solveOptions;1773 //solveOptions.setPresolveType(ClpSolve::presolveOff,0);1774 solveOptions.setSolveType(ClpSolve::usePrimalorSprint);1775 //solveOptions.setSpecialOption(1,3,30); // sprint1776 int numberColumns = lpSolver>numberColumns();1777 int iColumn;1778 bool allSlack = true;1779 for (iColumn = 0; iColumn < numberColumns; iColumn++) {1780 if (lpSolver>getColumnStatus(iColumn) == ClpSimplex::basic) {1781 allSlack = false;1782 break;1783 }1784 }1785 if (allSlack)1786 solveOptions.setSpecialOption(1, 2, 50); // idiot1787 lpSolver>setInfeasibilityCost(1.0e11);1788 lpSolver>defaultFactorizationFrequency();1789 if (doAction != 11)1790 lpSolver>initialSolve(solveOptions);1791 double * columnLower = lpSolver>columnLower();1792 double * columnUpper = lpSolver>columnUpper();1793 double * fullSolution = lpSolver>primalColumnSolution();1794 const double * dj = lpSolver>dualColumnSolution();1795 int iPass = 0;1796 #define MAXPROB 21797 ClpSimplex models[MAXPROB];1798 int pass[MAXPROB];1799 int kPass = 1;1800 int kLayer = 0;1801 int skipZero = 0;1802 if (skipZero2 == 1)1803 skipZero2 = 40; //1;1804 /* 0 fixed to 0 by choice1805 1 lb of 1 by choice1806 2 fixed to 0 by another1807 3 as 2 but this go1808 1 free1809 */1810 char * state = new char [numberColumns];1811 for (iColumn = 0; iColumn < numberColumns; iColumn++)1812 state[iColumn] = 1;1813 while (true) {1814 double largest = 0.1;1815 double smallest = 1.1;1816 int iLargest = 1;1817 int iSmallest = 1;1818 int atZero = 0;1819 int atOne = 0;1820 int toZero = 0;1821 int toOne = 0;1822 int numberFree = 0;1823 int numberGreater = 0;1824 columnLower = lpSolver>columnLower();1825 columnUpper = lpSolver>columnUpper();1826 fullSolution = lpSolver>primalColumnSolution();1827 if (doAction == 11) {1828 {1829 double * columnLower = lpSolver>columnLower();1830 double * columnUpper = lpSolver>columnUpper();1831 // lpSolver>dual();1832 memcpy(columnLower, saveColumnLower, numberColumns*sizeof(double));1833 memcpy(columnUpper, saveColumnUpper, numberColumns*sizeof(double));1834 // lpSolver>dual();1835 int iColumn;1836 for (iColumn = 0; iColumn < numberColumns; iColumn++) {1837 if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e8) {1838 if (clpSolver>isInteger(iColumn)) {1839 double value = lastSolution[iColumn];1840 int iValue = static_cast<int> (value + 0.5);1841 assert (fabs(value  static_cast<double> (iValue)) < 1.0e3);1842 assert (iValue >= columnLower[iColumn] &&1843 iValue <= columnUpper[iColumn]);1844 columnLower[iColumn] = iValue;1845 columnUpper[iColumn] = iValue;1846 }1847 }1848 }1849 lpSolver>initialSolve(solveOptions);1850 memcpy(columnLower, saveColumnLower, numberColumns*sizeof(double));1851 memcpy(columnUpper, saveColumnUpper, numberColumns*sizeof(double));1852 }1853 for (iColumn = 0; iColumn < numberColumns; iColumn++) {1854 if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e8) {1855 if (clpSolver>isInteger(iColumn)) {1856 double value = lastSolution[iColumn];1857 int iValue = static_cast<int> (value + 0.5);1858 assert (fabs(value  static_cast<double> (iValue)) < 1.0e3);1859 assert (iValue >= columnLower[iColumn] &&1860 iValue <= columnUpper[iColumn]);1861 if (!fix[iColumn]) {1862 if (iValue == 0) {1863 state[iColumn] = 0;1864 assert (!columnLower[iColumn]);1865 columnUpper[iColumn] = 0.0;1866 } else if (iValue == 1) {1867 state[iColumn] = 1;1868 columnLower[iColumn] = 1.0;1869 } else {1870 // leave fixed1871 columnLower[iColumn] = iValue;1872 columnUpper[iColumn] = iValue;1873 }1874 } else if (iValue == 0) {1875 state[iColumn] = 10;1876 columnUpper[iColumn] = 0.0;1877 } else {1878 // leave fixed1879 columnLower[iColumn] = iValue;1880 columnUpper[iColumn] = iValue;1881 }1882 }1883 }1884 }1885 int jLayer = 0;1886 int nFixed = 1;1887 int nTotalFixed = 0;1888 while (nFixed) {1889 nFixed = 0;1890 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {1891 if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {1892 for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {1893 int jColumn = otherColumn[i];1894 if (columnUpper[jColumn]) {1895 bool canFix = true;1896 for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {1897 int kColumn = otherColumn2[k];1898 if (state[kColumn] == 1) {1899 canFix = false;1900 break;1901 }1902 }1903 if (canFix) {1904 columnUpper[jColumn] = 0.0;1905 nFixed++;1906 }1907 }1908 }1909 }1910 }1911 nTotalFixed += nFixed;1912 jLayer += 100;1913 }1914 printf("This fixes %d variables in lower priorities\n", nTotalFixed);1915 break;1916 }1917 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {1918 if (!clpSolver>isInteger(iColumn)  fix[iColumn] > kLayer)1919 continue;1920 // skip if fixes nothing1921 if (fixColumn[iColumn+1]  fixColumn[iColumn] <= skipZero2)1922 continue;1923 double value = fullSolution[iColumn];1924 if (value > 1.00001) {1925 numberGreater++;1926 continue;1927 }1928 double lower = columnLower[iColumn];1929 double upper = columnUpper[iColumn];1930 if (lower == upper) {1931 if (lower)1932 atOne++;1933 else1934 atZero++;1935 continue;1936 }1937 if (value < 1.0e7) {1938 toZero++;1939 columnUpper[iColumn] = 0.0;1940 state[iColumn] = 10;1941 continue;1942 }1943 if (value > 1.0  1.0e7) {1944 toOne++;1945 columnLower[iColumn] = 1.0;1946 state[iColumn] = 1;1947 continue;1948 }1949 numberFree++;1950 // skip if fixes nothing1951 if (fixColumn[iColumn+1]  fixColumn[iColumn] <= skipZero)1952 continue;1953 if (value < smallest) {1954 smallest = value;1955 iSmallest = iColumn;1956 }1957 if (value > largest) {1958 largest = value;1959 iLargest = iColumn;1960 }1961 }1962 if (toZero  toOne)1963 printf("%d at 0 fixed and %d at one fixed\n", toZero, toOne);1964 printf("%d variables free, %d fixed to 0, %d to 1  smallest %g, largest %g\n",1965 numberFree, atZero, atOne, smallest, largest);1966 if (numberGreater && !iPass)1967 printf("%d variables have value > 1.0\n", numberGreater);1968 //skipZero2=0; // leave 0 fixing1969 int jLayer = 0;1970 int nFixed = 1;1971 int nTotalFixed = 0;1972 while (nFixed) {1973 nFixed = 0;1974 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {1975 if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {1976 for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {1977 int jColumn = otherColumn[i];1978 if (columnUpper[jColumn]) {1979 bool canFix = true;1980 for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {1981 int kColumn = otherColumn2[k];1982 if (state[kColumn] == 1) {1983 canFix = false;1984 break;1985 }1986 }1987 if (canFix) {1988 columnUpper[jColumn] = 0.0;1989 nFixed++;1990 }1991 }1992 }1993 }1994 }1995 nTotalFixed += nFixed;1996 jLayer += 100;1997 }1998 printf("This fixes %d variables in lower priorities\n", nTotalFixed);1999 if (iLargest < 0  numberFree <= leaveIntFree)2000 break;2001 double movement;2002 int way;2003 if (smallest <= 1.0  largest && smallest < 0.2 && largest < fixAboveValue) {2004 columnUpper[iSmallest] = 0.0;2005 state[iSmallest] = 0;2006 movement = smallest;2007 way = 1;2008 } else {2009 columnLower[iLargest] = 1.0;2010 state[iLargest] = 1;2011 movement = 1.0  largest;2012 way = 1;2013 }2014 double saveObj = lpSolver>objectiveValue();2015 iPass++;2016 kPass = iPass % MAXPROB;2017 models[kPass] = *lpSolver;2018 if (way == 1) {2019 // fix others2020 for (int i = fixColumn[iSmallest]; i < fixColumn[iSmallest+1]; i++) {2021 int jColumn = otherColumn[i];2022 if (state[jColumn] == 1) {2023 columnUpper[jColumn] = 0.0;2024 state[jColumn] = 3;2025 }2026 }2027 }2028 pass[kPass] = iPass;2029 double maxCostUp = COIN_DBL_MAX;2030 objective = lpSolver>getObjCoefficients() ;2031 if (way == 1)2032 maxCostUp = (1.0  movement) * objective[iSmallest];2033 lpSolver>setDualObjectiveLimit(saveObj + maxCostUp);2034 crunchIt(lpSolver);2035 double moveObj = lpSolver>objectiveValue()  saveObj;2036 printf("movement %s was %g costing %g\n",2037 (way == 1) ? "down" : "up", movement, moveObj);2038 if (way == 1 && (moveObj >= maxCostUp  lpSolver>status())) {2039 // go up2040 columnLower = models[kPass].columnLower();2041 columnUpper = models[kPass].columnUpper();2042 columnLower[iSmallest] = 1.0;2043 columnUpper[iSmallest] = saveColumnUpper[iSmallest];2044 *lpSolver = models[kPass];2045 columnLower = lpSolver>columnLower();2046 columnUpper = lpSolver>columnUpper();2047 fullSolution = lpSolver>primalColumnSolution();2048 dj = lpSolver>dualColumnSolution();2049 columnLower[iSmallest] = 1.0;2050 columnUpper[iSmallest] = saveColumnUpper[iSmallest];2051 state[iSmallest] = 1;2052 // unfix others2053 for (int i = fixColumn[iSmallest]; i < fixColumn[iSmallest+1]; i++) {2054 int jColumn = otherColumn[i];2055 if (state[jColumn] == 3) {2056 columnUpper[jColumn] = saveColumnUpper[jColumn];2057 state[jColumn] = 1;2058 }2059 }2060 crunchIt(lpSolver);2061 }2062 models[kPass] = *lpSolver;2063 }2064 lpSolver>dual();2065 printf("Fixing took %g seconds\n", CoinCpuTime()  time1);2066 columnLower = lpSolver>columnLower();2067 columnUpper = lpSolver>columnUpper();2068 fullSolution = lpSolver>primalColumnSolution();2069 dj = lpSolver>dualColumnSolution();2070 int * sort = new int[numberColumns];2071 double * dsort = new double[numberColumns];2072 int chunk = 20;2073 int iRelax = 0;2074 //double fractionFixed=6.0/8.0;2075 // relax while lots fixed2076 while (true) {2077 if (skipZero2 > 10 && doAction < 10)2078 break;2079 iRelax++;2080 int n = 0;2081 double sum0 = 0.0;2082 double sum00 = 0.0;2083 double sum1 = 0.0;2084 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {2085 if (!clpSolver>isInteger(iColumn)  fix[iColumn] > kLayer)2086 continue;2087 // skip if fixes nothing2088 if (fixColumn[iColumn+1]  fixColumn[iColumn] == 0 && doAction < 10)2089 continue;2090 double djValue = dj[iColumn];2091 if (state[iColumn] == 1) {2092 assert (columnLower[iColumn]);2093 assert (fullSolution[iColumn] > 0.1);2094 if (djValue > 0.0) {2095 //printf("YY dj of %d at %g is %g\n",iColumn,value,djValue);2096 sum1 += djValue;2097 sort[n] = iColumn;2098 dsort[n++] = djValue;2099 } else {2100 //printf("dj of %d at %g is %g\n",iColumn,value,djValue);2101 }2102 } else if (state[iColumn] == 0  state[iColumn] == 10) {2103 assert (fullSolution[iColumn] < 0.1);2104 assert (!columnUpper[iColumn]);2105 double otherValue = 0.0;2106 int nn = 0;2107 for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {2108 int jColumn = otherColumn[i];2109 if (columnUpper[jColumn] == 0.0) {2110 if (dj[jColumn] < 1.0e5) {2111 nn++;2112 otherValue += dj[jColumn]; // really need to look at rest2113 }2114 }2115 }2116 if (djValue < 1.0e2  otherValue < 1.0e2) {2117 //printf("XX dj of %d at %g is %g  %d out of %d contribute %g\n",iColumn,value,djValue,2118 // nn,fixColumn[iColumn+1]fixColumn[iColumn],otherValue);2119 if (djValue < 1.0e8) {2120 sum0 = djValue;2121 sum00 = otherValue;2122 sort[n] = iColumn;2123 if (djValue < 1.0e2)2124 dsort[n++] = djValue + otherValue;2125 else2126 dsort[n++] = djValue + 0.001 * otherValue;2127 }2128 } else {2129 //printf("dj of %d at %g is %g  no contribution from %d\n",iColumn,value,djValue,2130 // fixColumn[iColumn+1]fixColumn[iColumn]);2131 }2132 }2133 }2134 CoinSort_2(dsort, dsort + n, sort);2135 double * originalColumnLower = saveColumnLower;2136 double * originalColumnUpper = saveColumnUpper;2137 double * lo = CoinCopyOfArray(columnLower, numberColumns);2138 double * up = CoinCopyOfArray(columnUpper, numberColumns);2139 for (int k = 0; k < CoinMin(chunk, n); k++) {2140 iColumn = sort[k];2141 state[iColumn] = 2;2142 }2143 memcpy(columnLower, originalColumnLower, numberColumns*sizeof(double));2144 memcpy(columnUpper, originalColumnUpper, numberColumns*sizeof(double));2145 int nFixed = 0;2146 int nFixed0 = 0;2147 int nFixed1 = 0;2148 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {2149 if (state[iColumn] == 0  state[iColumn] == 10) {2150 columnUpper[iColumn] = 0.0;2151 assert (lo[iColumn] == 0.0);2152 nFixed++;2153 nFixed0++;2154 for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {2155 int jColumn = otherColumn[i];2156 if (columnUpper[jColumn]) {2157 bool canFix = true;2158 for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {2159 int kColumn = otherColumn2[k];2160 if (state[kColumn] == 1  state[kColumn] == 2) {2161 canFix = false;2162 break;2163 }2164 }2165 if (canFix) {2166 columnUpper[jColumn] = 0.0;2167 assert (lo[jColumn] == 0.0);2168 nFixed++;2169 }2170 }2171 }2172 } else if (state[iColumn] == 1) {2173 columnLower[iColumn] = 1.0;2174 nFixed1++;2175 }2176 }2177 printf("%d fixed %d orig 0 %d 1\n", nFixed, nFixed0, nFixed1);2178 int jLayer = 0;2179 nFixed = 1;2180 int nTotalFixed = 0;2181 while (nFixed) {2182 nFixed = 0;2183 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {2184 if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {2185 for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {2186 int jColumn = otherColumn[i];2187 if (columnUpper[jColumn]) {2188 bool canFix = true;2189 for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {2190 int kColumn = otherColumn2[k];2191 if (state[kColumn] == 1  state[kColumn] == 2) {2192 canFix = false;2193 break;2194 }2195 }2196 if (canFix) {2197 columnUpper[jColumn] = 0.0;2198 assert (lo[jColumn] == 0.0);2199 nFixed++;2200 }2201 }2202 }2203 }2204 }2205 nTotalFixed += nFixed;2206 jLayer += 100;2207 }2208 nFixed = 0;2209 int nFixedI = 0;2210 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {2211 if (columnLower[iColumn] == columnUpper[iColumn]) {2212 if (clpSolver>isInteger(iColumn))2213 nFixedI++;2214 nFixed++;2215 }2216 }2217 printf("This fixes %d variables in lower priorities  total %d (%d integer)  all target %d, int target %d\n",2218 nTotalFixed, nFixed, nFixedI, static_cast<int>(fractionFixed*numberColumns), static_cast<int> (fractionIntFixed*numberInteger));2219 int nBad = 0;2220 int nRelax = 0;2221 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {2222 if (lo[iColumn] < columnLower[iColumn] 2223 up[iColumn] > columnUpper[iColumn]) {2224 printf("bad %d old %g %g, new %g %g\n", iColumn, lo[iColumn], up[iColumn],2225 columnLower[iColumn], columnUpper[iColumn]);2226 nBad++;2227 }2228 if (lo[iColumn] > columnLower[iColumn] 2229 up[iColumn] < columnUpper[iColumn]) {2230 nRelax++;2231 }2232 }2233 printf("%d relaxed\n", nRelax);2234 if (iRelax > 20 && nRelax == chunk)2235 nRelax = 0;2236 if (iRelax > 50)2237 nRelax = 0;2238 assert (!nBad);2239 delete [] lo;2240 delete [] up;2241 lpSolver>primal(1);2242 if (nFixed < fractionFixed*numberColumns  nFixedI < fractionIntFixed*numberInteger  !nRelax)2243 break;2244 }2245 delete [] state;2246 delete [] sort;2247 delete [] dsort;2248 }2249 delete [] fix;2250 delete [] fixColumn;2251 delete [] otherColumn;2252 delete [] otherColumn2;2253 delete [] fixColumn2;2254 // See if was presolved2255 if (originalColumns) {2256 columnLower = lpSolver>columnLower();2257 columnUpper = lpSolver>columnUpper();2258 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {2259 saveColumnLower[iColumn] = columnLower[iColumn];2260 saveColumnUpper[iColumn] = columnUpper[iColumn];2261 }2262 pinfo.postsolve(true);2263 columnLower = originalLpSolver>columnLower();2264 columnUpper = originalLpSolver>columnUpper();2265 double * newColumnLower = lpSolver>columnLower();2266 double * newColumnUpper = lpSolver>columnUpper();2267 for (iColumn = 0; iColumn < numberColumns; iColumn++) {2268 int jColumn = originalColumns[iColumn];2269 columnLower[jColumn] = CoinMax(columnLower[jColumn], newColumnLower[iColumn]);2270 columnUpper[jColumn] = CoinMin(columnUpper[jColumn], newColumnUpper[iColumn]);2271 }2272 numberColumns = originalLpSolver>numberColumns();2273 delete [] originalColumns;2274 }2275 delete [] saveColumnLower;2276 delete [] saveColumnUpper;2277 if (!originalColumns) {2278 // Basis2279 memcpy(originalLpSolver>statusArray(), lpSolver>statusArray(), numberRows + numberColumns);2280 memcpy(originalLpSolver>primalColumnSolution(), lpSolver>primalColumnSolution(), numberColumns*sizeof(double));2281 memcpy(originalLpSolver>primalRowSolution(), lpSolver>primalRowSolution(), numberRows*sizeof(double));2282 // Fix in solver2283 columnLower = lpSolver>columnLower();2284 columnUpper = lpSolver>columnUpper();2285 }2286 double * originalColumnLower = originalLpSolver>columnLower();2287 double * originalColumnUpper = originalLpSolver>columnUpper();2288 // number fixed2289 doAction = 0;2290 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {2291 originalColumnLower[iColumn] = columnLower[iColumn];2292 originalColumnUpper[iColumn] = columnUpper[iColumn];2293 if (columnLower[iColumn] == columnUpper[iColumn])2294 doAction++;2295 }2296 printf("%d fixed by vub preprocessing\n", doAction);2297 if (originalColumns) {2298 originalLpSolver>initialSolve();2299 }2300 delete clpSolver;2301 return NULL;2302 }2303 #endif2304 1213 #ifdef COIN_HAS_LINK 2305 1214 /* Returns OsiSolverInterface (User should delete) … … 3157 2066 static CbcOrClpParam parameters[CBCMAXPARAMETERS]; 3158 2067 static int numberParameters = 0 ; 2068 2069 3159 2070 void CbcMain0 (CbcModel & model) 3160 2071 { … … 3308 2219 3  for miplib test so skip some 3309 2220 */ 3310 #if NEW_STYLE_SOLVER==03311 static int doHeuristics(CbcModel * model, int type)3312 #else3313 int3314 CbcSolver::doHeuristics(CbcModel * model, int type)3315 #endif3316 {3317 #if NEW_STYLE_SOLVER==03318 CbcOrClpParam * parameters_ = parameters;3319 int numberParameters_ = numberParameters;3320 bool noPrinting_ = noPrinting;3321 #endif3322 char generalPrint[10000];3323 CoinMessages generalMessages = model>messages();3324 CoinMessageHandler * generalMessageHandler = model>messageHandler();3325 //generalMessageHandler>setPrefix(false);3326 bool anyToDo = false;3327 int logLevel = parameters_[whichParam(CLP_PARAM_INT_LOGLEVEL, numberParameters_, parameters_)].intValue();3328 int useFpump = parameters_[whichParam(CBC_PARAM_STR_FPUMP, numberParameters_, parameters_)].currentOptionAsInteger();3329 int useRounding = parameters_[whichParam(CBC_PARAM_STR_ROUNDING, numberParameters_, parameters_)].currentOptionAsInteger();3330 int useGreedy = parameters_[whichParam(CBC_PARAM_STR_GREEDY, numberParameters_, parameters_)].currentOptionAsInteger();3331 int useCombine = parameters_[whichParam(CBC_PARAM_STR_COMBINE, numberParameters_, parameters_)].currentOptionAsInteger();3332 int useCrossover = parameters_[whichParam(CBC_PARAM_STR_CROSSOVER2, numberParameters_, parameters_)].currentOptionAsInteger();3333 //int usePivotC = parameters_[whichParam(CBC_PARAM_STR_PIVOTANDCOMPLEMENT, numberParameters_, parameters_)].currentOptionAsInteger();3334 int usePivotF = parameters_[whichParam(CBC_PARAM_STR_PIVOTANDFIX, numberParameters_, parameters_)].currentOptionAsInteger();3335 int useRand = parameters_[whichParam(CBC_PARAM_STR_RANDROUND, numberParameters_, parameters_)].currentOptionAsInteger();3336 int useRINS = parameters_[whichParam(CBC_PARAM_STR_RINS, numberParameters_, parameters_)].currentOptionAsInteger();3337 int useRENS = parameters_[whichParam(CBC_PARAM_STR_RENS, numberParameters_, parameters_)].currentOptionAsInteger();3338 int useDINS = parameters_[whichParam(CBC_PARAM_STR_DINS, numberParameters_, parameters_)].currentOptionAsInteger();3339 int useDIVING2 = parameters_[whichParam(CBC_PARAM_STR_DIVINGS, numberParameters_, parameters_)].currentOptionAsInteger();3340 int useNaive = parameters_[whichParam(CBC_PARAM_STR_NAIVE, numberParameters_, parameters_)].currentOptionAsInteger();3341 int kType = (type < 10) ? type : 1;3342 assert (kType == 1  kType == 2);3343 // FPump done first as it only works if no solution3344 if (useFpump >= kType && useFpump <= kType + 1) {3345 anyToDo = true;3346 CbcHeuristicFPump heuristic4(*model);3347 double dextra3 = parameters_[whichParam(CBC_PARAM_DBL_SMALLBAB, numberParameters_, parameters_)].doubleValue();3348 heuristic4.setFractionSmall(dextra3);3349 double dextra1 = parameters_[whichParam(CBC_PARAM_DBL_ARTIFICIALCOST, numberParameters_, parameters_)].doubleValue();3350 if (dextra1)3351 heuristic4.setArtificialCost(dextra1);3352 heuristic4.setMaximumPasses(parameters_[whichParam(CBC_PARAM_INT_FPUMPITS, numberParameters_, parameters_)].intValue());3353 if (parameters_[whichParam(CBC_PARAM_INT_FPUMPITS, numberParameters_, parameters_)].intValue() == 21)3354 heuristic4.setIterationRatio(1.0);3355 int pumpTune = parameters_[whichParam(CBC_PARAM_INT_FPUMPTUNE, numberParameters_, parameters_)].intValue();3356 int pumpTune2 = parameters_[whichParam(CBC_PARAM_INT_FPUMPTUNE2, numberParameters_, parameters_)].intValue();3357 if (pumpTune > 0) {3358 bool printStuff = (pumpTune != initialPumpTune  logLevel > 1  pumpTune2 > 0)3359 && !noPrinting_;3360 if (printStuff) {3361 generalMessageHandler>message(CBC_GENERAL, generalMessages)3362 << "Options for feasibility pump  "3363 << CoinMessageEol;3364 }3365 /*3366 >=10000000 for using obj3367 >=1000000 use as accumulate switch3368 >=1000 use index+1 as number of large loops3369 >=100 use dextra1 as cutoff3370 %100 == 10,20 etc for experimentation3371 1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds3372 4 and static continuous, 5 as 3 but no internal integers3373 6 as 3 but all slack basis!3374 */3375 double value = model>solver()>getObjSense() * model>solver()>getObjValue();3376 int w = pumpTune / 10;3377 int i = w % 10;3378 w /= 10;3379 int c = w % 10;3380 w /= 10;3381 int r = w;3382 int accumulate = r / 1000;3383 r = 1000 * accumulate;3384 if (accumulate >= 10) {3385 int which = accumulate / 10;3386 accumulate = 10 * which;3387 which;3388 // weights and factors3389 double weight[] = {0.01, 0.01, 0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};3390 double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};3391 heuristic4.setInitialWeight(weight[which]);3392 heuristic4.setWeightFactor(factor[which]);3393 if (printStuff) {3394 sprintf(generalPrint, "Initial weight for objective %g, decay factor %g",3395 weight[which], factor[which]);3396 generalMessageHandler>message(CBC_GENERAL, generalMessages)3397 << generalPrint3398 << CoinMessageEol;3399 }3400 3401 }3402 // fake cutoff3403 if (c) {3404 double cutoff;3405 model>solver()>getDblParam(OsiDualObjectiveLimit, cutoff);3406 cutoff = CoinMin(cutoff, value + 0.05 * fabs(value) * c);3407 double fakeCutoff = parameters_[whichParam(CBC_PARAM_DBL_FAKECUTOFF, numberParameters_, parameters_)].doubleValue();3408 if (fakeCutoff)3409 cutoff = fakeCutoff;3410 heuristic4.setFakeCutoff(cutoff);3411 if (printStuff) {3412 sprintf(generalPrint, "Fake cutoff of %g", cutoff);3413 generalMessageHandler>message(CBC_GENERAL, generalMessages)3414 << generalPrint3415 << CoinMessageEol;3416 }3417 }3418 int offRandomEtc = 0;3419 if (pumpTune2) {3420 if ((pumpTune2 / 1000) != 0) {3421 offRandomEtc = 1000000 * (pumpTune2 / 1000);3422 if (printStuff) {3423 generalMessageHandler>message(CBC_GENERAL, generalMessages)3424 << "Feasibility pump may run twice"3425 << CoinMessageEol;3426 }3427 pumpTune2 = pumpTune2 % 1000;3428 }3429 if ((pumpTune2 / 100) != 0) {3430 offRandomEtc += 100 * (pumpTune2 / 100);3431 if (printStuff) {3432 generalMessageHandler>message(CBC_GENERAL, generalMessages)3433 << "Not using randomized objective"3434 << CoinMessageEol;3435 }3436 }3437 int maxAllowed = pumpTune2 % 100;3438 if (maxAllowed) {3439 offRandomEtc += 1000 * maxAllowed;3440 if (printStuff) {3441 sprintf(generalPrint, "Fixing if same for %d passes",3442 maxAllowed);3443 generalMessageHandler>message(CBC_GENERAL, generalMessages)3444 << generalPrint3445 << CoinMessageEol;3446 }3447 }3448 }3449 if (accumulate) {3450 heuristic4.setAccumulate(accumulate);3451 if (printStuff) {3452 if (accumulate) {3453 sprintf(generalPrint, "Accumulate of %d", accumulate);3454 generalMessageHandler>message(CBC_GENERAL, generalMessages)3455 << generalPrint3456 << CoinMessageEol;3457 }3458 }3459 }3460 if (r) {3461 // also set increment3462 //double increment = (0.01*i+0.005)*(fabs(value)+1.0e12);3463 double increment = 0.0;3464 double fakeIncrement = parameters_[whichParam(CBC_PARAM_DBL_FAKEINCREMENT, numberParameters_, parameters_)].doubleValue();3465 if (fakeIncrement)3466 increment = fakeIncrement;3467 heuristic4.setAbsoluteIncrement(increment);3468 heuristic4.setMaximumRetries(r + 1);3469 if (printStuff) {3470 if (increment) {3471 sprintf(generalPrint, "Increment of %g", increment);3472 generalMessageHandler>message(CBC_GENERAL, generalMessages)3473 << generalPrint3474 << CoinMessageEol;3475 }3476 sprintf(generalPrint, "%d retries", r + 1);3477 generalMessageHandler>message(CBC_GENERAL, generalMessages)3478 << generalPrint3479 << CoinMessageEol;3480 }3481 }3482 if (i + offRandomEtc) {3483 heuristic4.setFeasibilityPumpOptions(i*10 + offRandomEtc);3484 if (printStuff) {3485 sprintf(generalPrint, "Feasibility pump options of %d",3486 i*10 + offRandomEtc);3487 generalMessageHandler>message(CBC_GENERAL, generalMessages)3488 << generalPrint3489 << CoinMessageEol;3490 }3491 }3492 pumpTune = pumpTune % 100;3493 if (pumpTune == 6)3494 pumpTune = 13;3495 heuristic4.setWhen((pumpTune % 10) + 10);3496 if (printStuff) {3497 sprintf(generalPrint, "Tuning (fixing) %d", pumpTune % 10);3498 generalMessageHandler>message(CBC_GENERAL, generalMessages)3499 << generalPrint3500 << CoinMessageEol;3501 }3502 }3503 heuristic4.setHeuristicName("feasibility pump");3504 //#define ROLF3505 #ifdef ROLF3506 CbcHeuristicFPump pump(*model);3507 pump.setMaximumTime(60);3508 pump.setMaximumPasses(100);3509 pump.setMaximumRetries(1);3510 pump.setFixOnReducedCosts(0);3511 pump.setHeuristicName("Feasibility pump");3512 pump.setFractionSmall(1.0);3513 pump.setWhen(13);3514 model>addHeuristic(&pump);3515 #else3516 model>addHeuristic(&heuristic4);3517 #endif3518 }3519 if (useRounding >= type && useRounding >= kType && useRounding <= kType + 1) {3520 CbcRounding heuristic1(*model);3521 heuristic1.setHeuristicName("rounding");3522 model>addHeuristic(&heuristic1) ;3523 anyToDo = true;3524 }3525 if (useGreedy >= type && useGreedy >= kType && useGreedy <= kType + 1) {3526 CbcHeuristicGreedyCover heuristic3(*model);3527 heuristic3.setHeuristicName("greedy cover");3528 CbcHeuristicGreedyEquality heuristic3a(*model);3529 heuristic3a.setHeuristicName("greedy equality");3530 model>addHeuristic(&heuristic3);3531 model>addHeuristic(&heuristic3a);3532 anyToDo = true;3533 }3534 if (useRENS >= kType && useRENS <= kType + 1) {3535 #if 13536 CbcHeuristicRENS heuristic6(*model);3537 heuristic6.setHeuristicName("RENS");3538 heuristic6.setFractionSmall(0.4);3539 heuristic6.setFeasibilityPumpOptions(1008003);3540 int nodes [] = { 2, 50, 50, 50, 200, 1000, 10000};3541 heuristic6.setNumberNodes(nodes[useRENS]);3542 #else3543 CbcHeuristicVND heuristic6(*model);3544 heuristic6.setHeuristicName("VND");3545 heuristic5.setFractionSmall(0.5);3546 heuristic5.setDecayFactor(5.0);3547 #endif3548 model>addHeuristic(&heuristic6) ;3549 anyToDo = true;3550 }3551 if (useNaive >= kType && useNaive <= kType + 1) {3552 CbcHeuristicNaive heuristic5b(*model);3553 heuristic5b.setHeuristicName("Naive");3554 heuristic5b.setFractionSmall(0.4);3555 heuristic5b.setNumberNodes(50);3556 model>addHeuristic(&heuristic5b) ;3557 anyToDo = true;3558 }3559 int useDIVING = 0;3560 {3561 int useD;3562 useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGV, numberParameters_, parameters_)].currentOptionAsInteger();3563 useDIVING = 1 * ((useD >= kType) ? 1 : 0);3564 useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGG, numberParameters_, parameters_)].currentOptionAsInteger();3565 useDIVING = 2 * ((useD >= kType) ? 1 : 0);3566 useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGF, numberParameters_, parameters_)].currentOptionAsInteger();3567 useDIVING = 4 * ((useD >= kType) ? 1 : 0);3568 useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGC, numberParameters_, parameters_)].currentOptionAsInteger();3569 useDIVING = 8 * ((useD >= kType) ? 1 : 0);3570 useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGL, numberParameters_, parameters_)].currentOptionAsInteger();3571 useDIVING = 16 * ((useD >= kType) ? 1 : 0);3572 useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGP, numberParameters_, parameters_)].currentOptionAsInteger();3573 useDIVING = 32 * ((useD >= kType) ? 1 : 0);3574 }3575 if (useDIVING2 >= kType && useDIVING2 <= kType + 1) {3576 int diveOptions = parameters_[whichParam(CBC_PARAM_INT_DIVEOPT, numberParameters_, parameters_)].intValue();3577 if (diveOptions < 0  diveOptions > 10)3578 diveOptions = 2;3579 CbcHeuristicJustOne heuristicJustOne(*model);3580 heuristicJustOne.setHeuristicName("DiveAny");3581 heuristicJustOne.setWhen(diveOptions);3582 // add in others3583 CbcHeuristicDiveCoefficient heuristicDC(*model);3584 heuristicDC.setHeuristicName("DiveCoefficient");3585 heuristicJustOne.addHeuristic(&heuristicDC, 1.0) ;3586 CbcHeuristicDiveFractional heuristicDF(*model);3587 heuristicDF.setHeuristicName("DiveFractional");3588 heuristicJustOne.addHeuristic(&heuristicDF, 1.0) ;3589 CbcHeuristicDiveGuided heuristicDG(*model);3590 heuristicDG.setHeuristicName("DiveGuided");3591 heuristicJustOne.addHeuristic(&heuristicDG, 1.0) ;3592 CbcHeuristicDiveLineSearch heuristicDL(*model);3593 heuristicDL.setHeuristicName("DiveLineSearch");3594 heuristicJustOne.addHeuristic(&heuristicDL, 1.0) ;3595 CbcHeuristicDivePseudoCost heuristicDP(*model);3596 heuristicDP.setHeuristicName("DivePseudoCost");3597 heuristicJustOne.addHeuristic(&heuristicDP, 1.0) ;3598 CbcHeuristicDiveVectorLength heuristicDV(*model);3599 heuristicDV.setHeuristicName("DiveVectorLength");3600 heuristicJustOne.addHeuristic(&heuristicDV, 1.0) ;3601 // Now normalize probabilities3602 heuristicJustOne.normalizeProbabilities();3603 model>addHeuristic(&heuristicJustOne) ;3604 }3605 3606 if (useDIVING > 0) {3607 int diveOptions2 = parameters_[whichParam(CBC_PARAM_INT_DIVEOPT, numberParameters_, parameters_)].intValue();3608 int diveOptions;3609 if (diveOptions2 > 99) {3610 // switch on various active set stuff3611 diveOptions = diveOptions2 % 100;3612 diveOptions2 = diveOptions;3613 } else {3614 diveOptions = diveOptions2;3615 diveOptions2 = 0;3616 }3617 if (diveOptions < 0  diveOptions > 9)3618 diveOptions = 2;3619 if ((useDIVING&1) != 0) {3620 CbcHeuristicDiveVectorLength heuristicDV(*model);3621 heuristicDV.setHeuristicName("DiveVectorLength");3622 heuristicDV.setWhen(diveOptions);3623 model>addHeuristic(&heuristicDV) ;3624 }3625 if ((useDIVING&2) != 0) {3626 CbcHeuristicDiveGuided heuristicDG(*model);3627 heuristicDG.setHeuristicName("DiveGuided");3628 heuristicDG.setWhen(diveOptions);3629 model>addHeuristic(&heuristicDG) ;3630 }3631 if ((useDIVING&4) != 0) {3632 CbcHeuristicDiveFractional heuristicDF(*model);3633 heuristicDF.setHeuristicName("DiveFractional");3634 heuristicDF.setWhen(diveOptions);3635 model>addHeuristic(&heuristicDF) ;3636 }3637 if ((useDIVING&8) != 0) {3638 CbcHeuristicDiveCoefficient heuristicDC(*model);3639 heuristicDC.setHeuristicName("DiveCoefficient");3640 heuristicDC.setWhen(diveOptions);3641 model>addHeuristic(&heuristicDC) ;3642 }3643 if ((useDIVING&16) != 0) {3644 CbcHeuristicDiveLineSearch heuristicDL(*model);3645 heuristicDL.setHeuristicName("DiveLineSearch");3646 heuristicDL.setWhen(diveOptions);3647 model>addHeuristic(&heuristicDL) ;3648 }3649 if ((useDIVING&32) != 0) {3650 CbcHeuristicDivePseudoCost heuristicDP(*model);3651 heuristicDP.setHeuristicName("DivePseudoCost");3652 heuristicDP.setWhen(diveOptions + diveOptions2);3653 model>addHeuristic(&heuristicDP) ;3654 }3655 anyToDo = true;3656 }3657 #if 03658 if (usePivotC >= type && usePivotC <= kType + 1) {3659 CbcHeuristicPivotAndComplement heuristic(*model);3660 heuristic.setHeuristicName("pivot and complement");3661 heuristic.setFractionSmall(10.0); // normally 0.53662 model>addHeuristic(&heuristic);3663 anyToDo = true;3664 }3665 #endif3666 if (usePivotF >= type && usePivotF <= kType + 1) {3667 CbcHeuristicPivotAndFix heuristic(*model);3668 heuristic.setHeuristicName("pivot and fix");3669 heuristic.setFractionSmall(10.0); // normally 0.53670 model>addHeuristic(&heuristic);3671 anyToDo = true;3672 }3673 if (useRand >= type && useRand <= kType + 1) {3674 CbcHeuristicRandRound heuristic(*model);3675 heuristic.setHeuristicName("randomized rounding");3676 heuristic.setFractionSmall(10.0); // normally 0.53677 model>addHeuristic(&heuristic);3678 anyToDo = true;3679 }3680 if (useDINS >= kType && useDINS <= kType + 1) {3681 CbcHeuristicDINS heuristic5a(*model);3682 heuristic5a.setHeuristicName("DINS");3683 heuristic5a.setFractionSmall(0.6);3684 if (useDINS < 4)3685 heuristic5a.setDecayFactor(5.0);3686 else3687 heuristic5a.setDecayFactor(1.5);3688 heuristic5a.setNumberNodes(1000);3689 model>addHeuristic(&heuristic5a) ;3690 anyToDo = true;3691 }3692 if (useRINS >= kType && useRINS <= kType + 1) {3693 CbcHeuristicRINS heuristic5(*model);3694 heuristic5.setHeuristicName("RINS");3695 if (useRINS < 4) {3696 heuristic5.setFractionSmall(0.5);3697 heuristic5.setDecayFactor(5.0);3698 } else {3699 heuristic5.setFractionSmall(0.6);3700 heuristic5.setDecayFactor(1.5);3701 }3702 model>addHeuristic(&heuristic5) ;3703 anyToDo = true;3704 }3705 if (useCombine >= kType && useCombine <= kType + 1) {3706 CbcHeuristicLocal heuristic2(*model);3707 heuristic2.setHeuristicName("combine solutions");3708 heuristic2.setFractionSmall(0.5);3709 heuristic2.setSearchType(1);3710 model>addHeuristic(&heuristic2);3711 anyToDo = true;3712 }3713 if (useCrossover >= kType && useCrossover <= kType + 1) {3714 CbcHeuristicCrossover heuristic2a(*model);3715 heuristic2a.setHeuristicName("crossover");3716 heuristic2a.setFractionSmall(0.3);3717 // just fix at lower3718 heuristic2a.setWhen(11);3719 model>addHeuristic(&heuristic2a);3720 model>setMaximumSavedSolutions(5);3721 anyToDo = true;3722 }3723 int heurSwitches = parameters_[whichParam(CBC_PARAM_INT_HOPTIONS, numberParameters_, parameters_)].intValue() % 100;3724 if (heurSwitches) {3725 for (int iHeur = 0; iHeur < model>numberHeuristics(); iHeur++) {3726 CbcHeuristic * heuristic = model>heuristic(iHeur);3727 heuristic>setSwitches(heurSwitches);3728 }3729 }3730 if (type == 2 && anyToDo) {3731 // Do heuristics3732 #if 13733 // clean copy3734 CbcModel model2(*model);3735 // But get rid of heuristics in model3736 model>doHeuristicsAtRoot(2);3737 if (logLevel <= 1)3738 model2.solver()>setHintParam(OsiDoReducePrint, true, OsiHintTry);3739 OsiBabSolver defaultC;3740 //solver_>setAuxiliaryInfo(&defaultC);3741 model2.passInSolverCharacteristics(&defaultC);3742 // Save bounds3743 int numberColumns = model2.solver()>getNumCols();3744 model2.createContinuousSolver();3745 bool cleanModel = !model2.numberIntegers() && !model2.numberObjects();3746 model2.findIntegers(false);3747 int heurOptions = (parameters_[whichParam(CBC_PARAM_INT_HOPTIONS, numberParameters_, parameters_)].intValue() / 100) % 100;3748 if (heurOptions == 0  heurOptions == 2) {3749 model2.doHeuristicsAtRoot(1);3750 } else if (heurOptions == 1  heurOptions == 3) {3751 model2.setMaximumNodes(1);3752 CbcStrategyDefault strategy(0, 5, 5);3753 strategy.setupPreProcessing(1, 0);3754 model2.setStrategy(strategy);3755 model2.branchAndBound();3756 }3757 if (cleanModel)3758 model2.zapIntegerInformation(false);3759 if (model2.bestSolution()) {3760 double value = model2.getMinimizationObjValue();3761 model>setCutoff(value);3762 model>setBestSolution(model2.bestSolution(), numberColumns, value);3763 model>setSolutionCount(1);3764 model>setNumberHeuristicSolutions(1);3765 }3766 #else3767 if (logLevel <= 1)3768 model>solver()>setHintParam(OsiDoReducePrint, true, OsiHintTry);3769 OsiBabSolver defaultC;3770 //solver_>setAuxiliaryInfo(&defaultC);3771 model>passInSolverCharacteristics(&defaultC);3772 // Save bounds3773 int numberColumns = model>solver()>getNumCols();3774 model>createContinuousSolver();3775 bool cleanModel = !model>numberIntegers() && !model>numberObjects();3776 model>findIntegers(false);3777 model>doHeuristicsAtRoot(1);3778 if (cleanModel)3779 model>zapIntegerInformation(false);3780 #endif3781 return 0;3782 } else {3783 return 0;3784 }3785 }3786 3787 2221 int CbcClpUnitTest (const CbcModel & saveModel, 3788 2222 std::string& dirMiplib, int testSwitch, … … 5716 4150 } 5717 4151 // Actually do heuristics 5718 doHeuristics(&model_, 2 );4152 doHeuristics(&model_, 2, parameters_, numberParameters_, noPrinting_, initialPumpTune); 5719 4153 if (model_.bestSolution()) { 5720 4154 model_.setProblemStatus(1); … … 6759 5193 } 6760 5194 // Set up heuristics 6761 doHeuristics(babModel_, (!miplib) ? 1 : 10 );5195 doHeuristics(babModel_, (!miplib) ? 1 : 10, parameters_, numberParameters_, noPrinting_, initialPumpTune); 6762 5196 if (!miplib) { 6763 5197 if (parameters_[whichParam(CBC_PARAM_STR_LOCALTREE, numberParameters_, parameters_)].currentOptionAsInteger()) { 
branches/sandbox/Cbc/src/CbcSolver.hpp
r1361 r1383 109 109 (out model later) 110 110 */ 111 int doHeuristics(CbcModel * model, int type);111 //int doHeuristics(CbcModel * model, int type); 112 112 /** Updates model_ from babModel_ according to returnMode 113 113 returnMode 
Note: See TracChangeset
for help on using the changeset viewer.