1574 | | } else if (!numberObjects_) { |
1575 | | // nothing to do |
1576 | | solverCharacteristics_ = NULL; |
1577 | | bestObjective_ = solver_->getObjValue() * solver_->getObjSense(); |
1578 | | int numberColumns = solver_->getNumCols(); |
1579 | | delete [] bestSolution_; |
1580 | | bestSolution_ = new double[numberColumns]; |
1581 | | CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_); |
1582 | | return ; |
1583 | | } |
1584 | | // Convert to Osi if wanted |
1585 | | bool useOsiBranching = false; |
1586 | | //OsiBranchingInformation * persistentInfo = NULL; |
1587 | | if (branchingMethod_ && branchingMethod_->chooseMethod()) { |
1588 | | useOsiBranching = true; |
1589 | | //persistentInfo = new OsiBranchingInformation(solver_); |
1590 | | if (numberOriginalObjects) { |
1591 | | for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) { |
1592 | | CbcObject * obj = |
1593 | | dynamic_cast <CbcObject *>(object_[iObject]) ; |
1594 | | if (obj) { |
1595 | | CbcSimpleInteger * obj2 = |
1596 | | dynamic_cast <CbcSimpleInteger *>(obj) ; |
1597 | | if (obj2) { |
1598 | | // back to Osi land |
1599 | | object_[iObject] = obj2->osiObject(); |
1600 | | delete obj; |
1601 | | } else { |
1602 | | OsiSimpleInteger * obj3 = |
1603 | | dynamic_cast <OsiSimpleInteger *>(obj) ; |
1604 | | if (!obj3) { |
1605 | | OsiSOS * obj4 = |
1606 | | dynamic_cast <OsiSOS *>(obj) ; |
1607 | | if (!obj4) { |
1608 | | CbcSOS * obj5 = |
1609 | | dynamic_cast <CbcSOS *>(obj) ; |
1610 | | if (obj5) { |
1611 | | // back to Osi land |
1612 | | object_[iObject] = obj5->osiObject(solver_); |
1613 | | } else { |
1614 | | printf("Code up CbcObject type in Osi land\n"); |
1615 | | abort(); |
1616 | | } |
1617 | | } |
1618 | | } |
1619 | | } |
1620 | | } |
1621 | | } |
1622 | | // and add to solver |
1623 | | //if (!solver_->numberObjects()) { |
1624 | | solver_->addObjects(numberObjects_, object_); |
1625 | | //} else { |
1626 | | //if (solver_->numberObjects()!=numberOriginalObjects) { |
1627 | | //printf("should have trapped that solver has objects before\n"); |
1628 | | //abort(); |
1629 | | //} |
1630 | | //} |
1631 | | } else { |
1632 | | // do from solver |
1633 | | deleteObjects(false); |
1634 | | solver_->findIntegersAndSOS(false); |
1635 | | numberObjects_ = solver_->numberObjects(); |
1636 | | object_ = solver_->objects(); |
1637 | | ownObjects_ = false; |
1638 | | } |
1639 | | branchingMethod_->chooseMethod()->setSolver(solver_); |
1640 | | } |
1641 | | // take off heuristics if have to |
1642 | | { |
1643 | | int numberOdd = 0; |
1644 | | int numberSOS = 0; |
1645 | | for (int i = 0; i < numberObjects_; i++) { |
1646 | | if (!object_[i]->canDoHeuristics()) |
1647 | | numberOdd++; |
1648 | | CbcSOS * obj = |
1649 | | dynamic_cast <CbcSOS *>(object_[i]) ; |
1650 | | if (obj) |
1651 | | numberSOS++; |
1652 | | } |
1653 | | if (numberOdd) { |
1654 | | if (numberHeuristics_) { |
1655 | | int k = 0; |
1656 | | for (int i = 0; i < numberHeuristics_; i++) { |
1657 | | if (!heuristic_[i]->canDealWithOdd()) |
1658 | | delete heuristic_[i]; |
1659 | | else |
1660 | | heuristic_[k++] = heuristic_[i]; |
1661 | | } |
1662 | | if (!k) { |
1663 | | delete [] heuristic_; |
1664 | | heuristic_ = NULL; |
1665 | | } |
1666 | | numberHeuristics_ = k; |
1667 | | handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol ; |
1668 | | } |
1669 | | } else if (numberSOS) { |
1670 | | specialOptions_ |= 128; // say can do SOS in dynamic mode |
1671 | | // switch off fast nodes for now |
1672 | | fastNodeDepth_ = -1; |
1673 | | } |
1674 | | if (numberThreads_ > 0) { |
1675 | | // switch off fast nodes for now |
1676 | | fastNodeDepth_ = -1; |
1677 | | } |
1678 | | } |
1679 | | // Save objective (just so user can access it) |
1680 | | originalContinuousObjective_ = solver_->getObjValue(); |
1681 | | bestPossibleObjective_ = originalContinuousObjective_; |
1682 | | sumChangeObjective1_ = 0.0; |
1683 | | sumChangeObjective2_ = 0.0; |
1684 | | /* |
1685 | | OsiRowCutDebugger knows an optimal answer for a subset of MIP problems. |
1686 | | Assuming it recognises the problem, when called upon it will check a cut to |
1687 | | see if it cuts off the optimal answer. |
1688 | | */ |
1689 | | // If debugger exists set specialOptions_ bit |
1690 | | if (solver_->getRowCutDebuggerAlways()) { |
1691 | | specialOptions_ |= 1; |
1692 | | } |
1693 | | |
1694 | | # ifdef CBC_DEBUG |
1695 | | if ((specialOptions_&1) == 0) |
1696 | | solver_->activateRowCutDebugger(problemName.c_str()) ; |
1697 | | if (solver_->getRowCutDebuggerAlways()) |
1698 | | specialOptions_ |= 1; |
1699 | | # endif |
1700 | | |
1701 | | /* |
1702 | | Begin setup to process a feasible root node. |
1703 | | */ |
1704 | | bestObjective_ = CoinMin(bestObjective_, 1.0e50) ; |
1705 | | if (!bestSolution_) { |
1706 | | numberSolutions_ = 0 ; |
1707 | | numberHeuristicSolutions_ = 0 ; |
1708 | | } |
1709 | | stateOfSearch_ = 0; |
1710 | | // Everything is minimization |
1711 | | { |
1712 | | // needed to sync cutoffs |
1713 | | double value ; |
1714 | | solver_->getDblParam(OsiDualObjectiveLimit, value) ; |
1715 | | dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense(); |
1716 | | } |
1717 | | double cutoff = getCutoff() ; |
1718 | | double direction = solver_->getObjSense() ; |
1719 | | dblParam_[CbcOptimizationDirection] = direction; |
1720 | | if (cutoff < 1.0e20 && direction < 0.0) |
1721 | | messageHandler()->message(CBC_CUTOFF_WARNING1, |
1722 | | messages()) |
1723 | | << cutoff << -cutoff << CoinMessageEol ; |
1724 | | if (cutoff > bestObjective_) |
1725 | | cutoff = bestObjective_ ; |
1726 | | setCutoff(cutoff) ; |
1727 | | /* |
1728 | | We probably already have a current solution, but just in case ... |
1729 | | */ |
1730 | | int numberColumns = getNumCols() ; |
1731 | | if (!currentSolution_) |
1732 | | currentSolution_ = new double[numberColumns] ; |
1733 | | testSolution_ = currentSolution_; |
1734 | | /* |
1735 | | Create a copy of the solver, thus capturing the original (root node) |
1736 | | constraint system (aka the continuous system). |
1737 | | */ |
1738 | | continuousSolver_ = solver_->clone() ; |
1739 | | |
1740 | | numberRowsAtContinuous_ = getNumRows() ; |
1741 | | solver_->saveBaseModel(); |
1742 | | /* |
1743 | | Check the objective to see if we can deduce a nontrivial increment. If |
1744 | | it's better than the current value for CbcCutoffIncrement, it'll be |
1745 | | installed. |
1746 | | */ |
1747 | | if (solverCharacteristics_->reducedCostsAccurate()) |
1748 | | analyzeObjective() ; |
1749 | | { |
1750 | | // may be able to change cutoff now |
1751 | | double cutoff = getCutoff(); |
1752 | | double increment = getDblParam(CbcModel::CbcCutoffIncrement) ; |
1753 | | if (cutoff > bestObjective_ - increment) { |
1754 | | cutoff = bestObjective_ - increment ; |
1755 | | setCutoff(cutoff) ; |
1756 | | } |
1757 | | } |
1758 | | #ifdef COIN_HAS_CLP |
1759 | | // Possible save of pivot method |
1760 | | ClpDualRowPivot * savePivotMethod = NULL; |
1761 | | { |
1762 | | // pass tolerance and increment to solver |
1763 | | OsiClpSolverInterface * clpSolver |
1764 | | = dynamic_cast<OsiClpSolverInterface *> (solver_); |
1765 | | if (clpSolver) |
1766 | | clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement()); |
1767 | | } |
1768 | | #endif |
1769 | | /* |
1770 | | Set up for cut generation. addedCuts_ holds the cuts which are relevant for |
1771 | | the active subproblem. whichGenerator will be used to record the generator |
1772 | | that produced a given cut. |
1773 | | */ |
1774 | | maximumWhich_ = 1000 ; |
1775 | | delete [] whichGenerator_; |
1776 | | whichGenerator_ = new int[maximumWhich_] ; |
1777 | | memset(whichGenerator_, 0, maximumWhich_*sizeof(int)); |
1778 | | maximumNumberCuts_ = 0 ; |
1779 | | currentNumberCuts_ = 0 ; |
1780 | | delete [] addedCuts_ ; |
1781 | | addedCuts_ = NULL ; |
1782 | | OsiObject ** saveObjects = NULL; |
1783 | | maximumRows_ = numberRowsAtContinuous_; |
1784 | | currentDepth_ = 0; |
1785 | | workingBasis_.resize(maximumRows_, numberColumns); |
1786 | | /* |
1787 | | Set up an empty heap and associated data structures to hold the live set |
1788 | | (problems which require further exploration). |
1789 | | */ |
1790 | | CbcCompareDefault * compareActual |
1791 | | = dynamic_cast<CbcCompareDefault *> (nodeCompare_); |
1792 | | if (compareActual) { |
1793 | | compareActual->setBestPossible(direction*solver_->getObjValue()); |
1794 | | compareActual->setCutoff(getCutoff()); |
1795 | | if (false && !numberThreads_ && !parentModel_) { |
1796 | | printf("CbcTreeArray ? threads ? parentArray\n"); |
1797 | | // Setup new style tree |
1798 | | delete tree_; |
1799 | | tree_ = new CbcTreeArray(); |
1800 | | } |
1801 | | } |
1802 | | tree_->setComparison(*nodeCompare_) ; |
1803 | | /* |
1804 | | Used to record the path from a node to the root of the search tree, so that |
1805 | | we can then traverse from the root to the node when restoring a subproblem. |
1806 | | */ |
1807 | | maximumDepth_ = 10 ; |
1808 | | delete [] walkback_ ; |
1809 | | walkback_ = new CbcNodeInfo * [maximumDepth_] ; |
1810 | | lastDepth_ = 0; |
1811 | | delete [] lastNodeInfo_ ; |
1812 | | lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ; |
1813 | | delete [] lastNumberCuts_ ; |
1814 | | lastNumberCuts_ = new int [maximumDepth_] ; |
1815 | | maximumCuts_ = 100; |
1816 | | lastNumberCuts2_ = 0; |
1817 | | delete [] lastCut_; |
1818 | | lastCut_ = new const OsiRowCut * [maximumCuts_]; |
1819 | | /* |
1820 | | Used to generate bound edits for CbcPartialNodeInfo. |
1821 | | */ |
1822 | | double * lowerBefore = new double [numberColumns] ; |
1823 | | double * upperBefore = new double [numberColumns] ; |
1824 | | /* |
1825 | | |
1826 | | Generate cuts at the root node and reoptimise. solveWithCuts does the heavy |
1827 | | lifting. It will iterate a generate/reoptimise loop (including reduced cost |
1828 | | fixing) until no cuts are generated, the change in objective falls off, or |
1829 | | the limit on the number of rounds of cut generation is exceeded. |
1830 | | |
1831 | | At the end of all this, any cuts will be recorded in cuts and also |
1832 | | installed in the solver's constraint system. We'll have reoptimised, and |
1833 | | removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been |
1834 | | adjusted accordingly). |
1835 | | |
1836 | | Tell cut generators they can be a bit more aggressive at root node |
1837 | | |
1838 | | TODO: Why don't we make a copy of the solution after solveWithCuts? |
1839 | | TODO: If numberUnsatisfied == 0, don't we have a solution? |
1840 | | */ |
1841 | | phase_ = 1; |
1842 | | int iCutGenerator; |
1843 | | for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) { |
1844 | | CglCutGenerator * generator = generator_[iCutGenerator]->generator(); |
1845 | | generator->setAggressiveness(generator->getAggressiveness() + 100); |
1846 | | } |
1847 | | OsiCuts cuts ; |
1848 | | int anyAction = -1 ; |
1849 | | numberOldActiveCuts_ = 0 ; |
1850 | | numberNewCuts_ = 0 ; |
1851 | | // Array to mark solution |
1852 | | delete [] usedInSolution_; |
1853 | | usedInSolution_ = new int[numberColumns]; |
1854 | | CoinZeroN(usedInSolution_, numberColumns); |
1855 | | /* |
1856 | | For printing totals and for CbcNode (numberNodes_) |
1857 | | */ |
1858 | | numberIterations_ = 0 ; |
1859 | | numberSolves_ = 0 ; |
1860 | | numberNodes_ = 0 ; |
1861 | | numberNodes2_ = 0 ; |
1862 | | maximumStatistics_ = 0; |
1863 | | maximumDepthActual_ = 0; |
1864 | | numberDJFixed_ = 0.0; |
1865 | | // Do heuristics |
1866 | | doHeuristicsAtRoot(); |
1867 | | if ( intParam_[CbcMaxNumNode] < 0) |
1868 | | eventHappened_ = true; // stop as fast as possible |
1869 | | stoppedOnGap_ = false ; |
1870 | | // See if can stop on gap |
1871 | | bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense(); |
1872 | | double testGap = CoinMax(dblParam_[CbcAllowableGap], |
1873 | | CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_)) |
1874 | | * dblParam_[CbcAllowableFractionGap]); |
1875 | | if (bestObjective_ - bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) { |
1876 | | if (bestPossibleObjective_ < getCutoff()) |
1877 | | stoppedOnGap_ = true ; |
1878 | | feasible = false; |
1879 | | //eventHappened_=true; // stop as fast as possible |
1880 | | } |
1881 | | statistics_ = NULL; |
1882 | | // Do on switch |
1883 | | if (doStatistics > 0 && doStatistics <= 100) { |
1884 | | maximumStatistics_ = 10000; |
1885 | | statistics_ = new CbcStatistics * [maximumStatistics_]; |
1886 | | memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *)); |
1887 | | } |
1888 | | // See if we can add integers |
1889 | | if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_) { |
1890 | | int numberColumns = continuousSolver_->getNumCols(); |
1891 | | int numberRows = continuousSolver_->getNumRows(); |
1892 | | int * del = new int [CoinMax(numberColumns, numberRows)]; |
1893 | | int * original = new int [numberColumns]; |
1894 | | char * possibleRow = new char [numberRows]; |
1895 | | { |
1896 | | const CoinPackedMatrix * rowCopy = continuousSolver_->getMatrixByRow(); |
1897 | | const int * column = rowCopy->getIndices(); |
1898 | | const int * rowLength = rowCopy->getVectorLengths(); |
1899 | | const CoinBigIndex * rowStart = rowCopy->getVectorStarts(); |
1900 | | const double * rowLower = continuousSolver_->getRowLower(); |
1901 | | const double * rowUpper = continuousSolver_->getRowUpper(); |
1902 | | const double * element = rowCopy->getElements(); |
1903 | | for (int i = 0; i < numberRows; i++) { |
1904 | | int nLeft = 0; |
1905 | | bool possible = false; |
1906 | | if (rowLower[i] < -1.0e20) { |
1907 | | double value = rowUpper[i]; |
1908 | | if (fabs(value - floor(value + 0.5)) < 1.0e-8) |
1909 | | possible = true; |
1910 | | } else if (rowUpper[i] > 1.0e20) { |
1911 | | double value = rowLower[i]; |
1912 | | if (fabs(value - floor(value + 0.5)) < 1.0e-8) |
1913 | | possible = true; |
1914 | | } else { |
1915 | | double value = rowUpper[i]; |
1916 | | if (rowLower[i] == rowUpper[i] && |
1917 | | fabs(value - floor(value + 0.5)) < 1.0e-8) |
1918 | | possible = true; |
1919 | | } |
1920 | | for (CoinBigIndex j = rowStart[i]; |
1921 | | j < rowStart[i] + rowLength[i]; j++) { |
1922 | | int iColumn = column[j]; |
1923 | | if (continuousSolver_->isInteger(iColumn)) { |
1924 | | if (fabs(element[j]) != 1.0) |
1925 | | possible = false; |
1926 | | } else { |
1927 | | nLeft++; |
1928 | | } |
1929 | | } |
1930 | | if (possible || !nLeft) |
1931 | | possibleRow[i] = 1; |
1932 | | else |
1933 | | possibleRow[i] = 0; |
1934 | | } |
1935 | | } |
1936 | | int nDel = 0; |
1937 | | for (int i = 0; i < numberColumns; i++) { |
1938 | | original[i] = i; |
1939 | | if (continuousSolver_->isInteger(i)) |
1940 | | del[nDel++] = i; |
1941 | | } |
1942 | | int nExtra = 0; |
1943 | | OsiSolverInterface * copy1 = continuousSolver_->clone(); |
1944 | | int nPass = 0; |
1945 | | while (nDel && nPass < 10) { |
1946 | | nPass++; |
1947 | | OsiSolverInterface * copy2 = copy1->clone(); |
1948 | | int nLeft = 0; |
1949 | | for (int i = 0; i < nDel; i++) |
1950 | | original[del[i]] = -1; |
1951 | | for (int i = 0; i < numberColumns; i++) { |
1952 | | int kOrig = original[i]; |
1953 | | if (kOrig >= 0) |
1954 | | original[nLeft++] = kOrig; |
1955 | | } |
1956 | | assert (nLeft == numberColumns - nDel); |
1957 | | copy2->deleteCols(nDel, del); |
1958 | | numberColumns = copy2->getNumCols(); |
1959 | | const CoinPackedMatrix * rowCopy = copy2->getMatrixByRow(); |
1960 | | numberRows = rowCopy->getNumRows(); |
1961 | | const int * column = rowCopy->getIndices(); |
1962 | | const int * rowLength = rowCopy->getVectorLengths(); |
1963 | | const CoinBigIndex * rowStart = rowCopy->getVectorStarts(); |
1964 | | const double * rowLower = copy2->getRowLower(); |
1965 | | const double * rowUpper = copy2->getRowUpper(); |
1966 | | const double * element = rowCopy->getElements(); |
1967 | | const CoinPackedMatrix * columnCopy = copy2->getMatrixByCol(); |
1968 | | const int * columnLength = columnCopy->getVectorLengths(); |
1969 | | nDel = 0; |
1970 | | // Could do gcd stuff on ones with costs |
1971 | | for (int i = 0; i < numberRows; i++) { |
1972 | | if (!rowLength[i]) { |
1973 | | del[nDel++] = i; |
1974 | | possibleRow[i] = 1; |
1975 | | } else if (possibleRow[i]) { |
1976 | | if (rowLength[i] == 1) { |
1977 | | int k = rowStart[i]; |
1978 | | int iColumn = column[k]; |
1979 | | if (!copy2->isInteger(iColumn)) { |
1980 | | double mult = 1.0 / fabs(element[k]); |
1981 | | if (rowLower[i] < -1.0e20) { |
1982 | | double value = rowUpper[i] * mult; |
1983 | | if (fabs(value - floor(value + 0.5)) < 1.0e-8) { |
1984 | | del[nDel++] = i; |
1985 | | if (columnLength[iColumn] == 1) { |
1986 | | copy2->setInteger(iColumn); |
1987 | | int kOrig = original[iColumn]; |
1988 | | setOptionalInteger(kOrig); |
1989 | | } |
1990 | | } |
1991 | | } else if (rowUpper[i] > 1.0e20) { |
1992 | | double value = rowLower[i] * mult; |
1993 | | if (fabs(value - floor(value + 0.5)) < 1.0e-8) { |
1994 | | del[nDel++] = i; |
1995 | | if (columnLength[iColumn] == 1) { |
1996 | | copy2->setInteger(iColumn); |
1997 | | int kOrig = original[iColumn]; |
1998 | | setOptionalInteger(kOrig); |
1999 | | } |
2000 | | } |
2001 | | } else { |
2002 | | double value = rowUpper[i] * mult; |
2003 | | if (rowLower[i] == rowUpper[i] && |
2004 | | fabs(value - floor(value + 0.5)) < 1.0e-8) { |
2005 | | del[nDel++] = i; |
2006 | | copy2->setInteger(iColumn); |
2007 | | int kOrig = original[iColumn]; |
2008 | | setOptionalInteger(kOrig); |
2009 | | } |
2010 | | } |
2011 | | } |
2012 | | } else { |
2013 | | // only if all singletons |
2014 | | bool possible = false; |
2015 | | if (rowLower[i] < -1.0e20) { |
2016 | | double value = rowUpper[i]; |
2017 | | if (fabs(value - floor(value + 0.5)) < 1.0e-8) |
2018 | | possible = true; |
2019 | | } else if (rowUpper[i] > 1.0e20) { |
2020 | | double value = rowLower[i]; |
2021 | | if (fabs(value - floor(value + 0.5)) < 1.0e-8) |
2022 | | possible = true; |
2023 | | } else { |
2024 | | double value = rowUpper[i]; |
2025 | | if (rowLower[i] == rowUpper[i] && |
2026 | | fabs(value - floor(value + 0.5)) < 1.0e-8) |
2027 | | possible = true; |
2028 | | } |
2029 | | if (possible) { |
2030 | | for (CoinBigIndex j = rowStart[i]; |
2031 | | j < rowStart[i] + rowLength[i]; j++) { |
2032 | | int iColumn = column[j]; |
2033 | | if (columnLength[iColumn] != 1 || fabs(element[j]) != 1.0) { |
2034 | | possible = false; |
2035 | | break; |
2036 | | } |
2037 | | } |
2038 | | if (possible) { |
2039 | | for (CoinBigIndex j = rowStart[i]; |
2040 | | j < rowStart[i] + rowLength[i]; j++) { |
2041 | | int iColumn = column[j]; |
2042 | | if (!copy2->isInteger(iColumn)) { |
2043 | | copy2->setInteger(iColumn); |
2044 | | int kOrig = original[iColumn]; |
2045 | | setOptionalInteger(kOrig); |
2046 | | } |
2047 | | } |
2048 | | del[nDel++] = i; |
2049 | | } |
2050 | | } |
2051 | | } |
2052 | | } |
2053 | | } |
2054 | | if (nDel) { |
2055 | | copy2->deleteRows(nDel, del); |
2056 | | } |
2057 | | if (nDel != numberRows) { |
2058 | | nDel = 0; |
2059 | | for (int i = 0; i < numberColumns; i++) { |
2060 | | if (copy2->isInteger(i)) { |
2061 | | del[nDel++] = i; |
2062 | | nExtra++; |
2063 | | } |
2064 | | } |
2065 | | } else { |
2066 | | nDel = 0; |
2067 | | } |
2068 | | delete copy1; |
2069 | | copy1 = copy2->clone(); |
2070 | | delete copy2; |
2071 | | } |
2072 | | // See if what's left is a network |
2073 | | bool couldBeNetwork = false; |
2074 | | if (copy1->getNumRows() && copy1->getNumCols()) { |
2075 | | #ifdef COIN_HAS_CLP |
2076 | | OsiClpSolverInterface * clpSolver |
2077 | | = dynamic_cast<OsiClpSolverInterface *> (copy1); |
2078 | | if (false && clpSolver) { |
2079 | | numberRows = clpSolver->getNumRows(); |
2080 | | char * rotate = new char[numberRows]; |
2081 | | int n = clpSolver->getModelPtr()->findNetwork(rotate, 1.0); |
2082 | | delete [] rotate; |
2083 | | #ifdef CLP_INVESTIGATE |
2084 | | printf("INTA network %d rows out of %d\n", n, numberRows); |
2085 | | #endif |
2086 | | if (CoinAbs(n) == numberRows) { |
2087 | | couldBeNetwork = true; |
2088 | | for (int i = 0; i < numberRows; i++) { |
2089 | | if (!possibleRow[i]) { |
2090 | | couldBeNetwork = false; |
2091 | | #ifdef CLP_INVESTIGATE |
2092 | | printf("but row %d is bad\n", i); |
2093 | | #endif |
2094 | | break; |
2095 | | } |
2096 | | } |
2097 | | } |
2098 | | } else |
2099 | | #endif |
2100 | | { |
2101 | | numberColumns = copy1->getNumCols(); |
2102 | | numberRows = copy1->getNumRows(); |
2103 | | const double * rowLower = copy1->getRowLower(); |
2104 | | const double * rowUpper = copy1->getRowUpper(); |
2105 | | couldBeNetwork = true; |
2106 | | for (int i = 0; i < numberRows; i++) { |
2107 | | if (rowLower[i] > -1.0e20 && fabs(rowLower[i] - floor(rowLower[i] + 0.5)) > 1.0e-12) { |
2108 | | couldBeNetwork = false; |
2109 | | break; |
2110 | | } |
2111 | | if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i] - floor(rowUpper[i] + 0.5)) > 1.0e-12) { |
2112 | | couldBeNetwork = false; |
2113 | | break; |
2114 | | } |
2115 | | } |
2116 | | if (couldBeNetwork) { |
2117 | | const CoinPackedMatrix * matrixByCol = copy1->getMatrixByCol(); |
2118 | | const double * element = matrixByCol->getElements(); |
2119 | | //const int * row = matrixByCol->getIndices(); |
2120 | | const CoinBigIndex * columnStart = matrixByCol->getVectorStarts(); |
2121 | | const int * columnLength = matrixByCol->getVectorLengths(); |
2122 | | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
2123 | | CoinBigIndex start = columnStart[iColumn]; |
2124 | | CoinBigIndex end = start + columnLength[iColumn]; |
2125 | | if (end > start + 2) { |
2126 | | couldBeNetwork = false; |
2127 | | break; |
2128 | | } |
2129 | | int type = 0; |
2130 | | for (CoinBigIndex j = start; j < end; j++) { |
2131 | | double value = element[j]; |
2132 | | if (fabs(value) != 1.0) { |
2133 | | couldBeNetwork = false; |
2134 | | break; |
2135 | | } else if (value == 1.0) { |
2136 | | if ((type&1) == 0) |
2137 | | type |= 1; |
2138 | | else |
2139 | | type = 7; |
2140 | | } else if (value == -1.0) { |
2141 | | if ((type&2) == 0) |
2142 | | type |= 2; |
2143 | | else |
2144 | | type = 7; |
2145 | | } |
2146 | | } |
2147 | | if (type > 3) { |
2148 | | couldBeNetwork = false; |
2149 | | break; |
2150 | | } |
2151 | | } |
2152 | | } |
2153 | | } |
2154 | | } |
2155 | | if (couldBeNetwork) { |
2156 | | for (int i = 0; i < numberColumns; i++) |
2157 | | setOptionalInteger(original[i]); |
2158 | | } |
2159 | | if (nExtra || couldBeNetwork) { |
2160 | | numberColumns = copy1->getNumCols(); |
2161 | | numberRows = copy1->getNumRows(); |
2162 | | if (!numberColumns || !numberRows) { |
2163 | | int numberColumns = solver_->getNumCols(); |
2164 | | for (int i = 0; i < numberColumns; i++) |
2165 | | assert(solver_->isInteger(i)); |
2166 | | } |
2167 | | #ifdef CLP_INVESTIGATE |
2168 | | if (couldBeNetwork || nExtra) |
2169 | | printf("INTA %d extra integers, %d left%s\n", nExtra, |
2170 | | numberColumns, |
2171 | | couldBeNetwork ? ", all network" : ""); |
2172 | | #endif |
2173 | | findIntegers(true, 2); |
2174 | | convertToDynamic(); |
2175 | | } |
2176 | | #ifdef CLP_INVESTIGATE |
2177 | | if (!couldBeNetwork && copy1->getNumCols() && |
2178 | | copy1->getNumRows()) { |
2179 | | printf("INTA %d rows and %d columns remain\n", |
2180 | | copy1->getNumRows(), copy1->getNumCols()); |
2181 | | if (copy1->getNumCols() < 200) { |
2182 | | copy1->writeMps("moreint"); |
2183 | | printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n", |
2184 | | copy1->getNumRows(), copy1->getNumCols()); |
2185 | | } |
2186 | | } |
2187 | | #endif |
2188 | | delete copy1; |
2189 | | delete [] del; |
2190 | | delete [] original; |
2191 | | delete [] possibleRow; |
2192 | | // double check increment |
2193 | | analyzeObjective(); |
2194 | | } |
2195 | | { |
2196 | | int iObject ; |
2197 | | int preferredWay ; |
2198 | | int numberUnsatisfied = 0 ; |
2199 | | delete [] currentSolution_; |
2200 | | currentSolution_ = new double [numberColumns]; |
2201 | | testSolution_ = currentSolution_; |
2202 | | memcpy(currentSolution_, solver_->getColSolution(), |
2203 | | numberColumns*sizeof(double)) ; |
2204 | | // point to useful information |
2205 | | OsiBranchingInformation usefulInfo = usefulInformation(); |
2206 | | |
2207 | | for (iObject = 0 ; iObject < numberObjects_ ; iObject++) { |
2208 | | double infeasibility = |
2209 | | object_[iObject]->infeasibility(&usefulInfo, preferredWay) ; |
2210 | | if (infeasibility ) numberUnsatisfied++ ; |
2211 | | } |
2212 | | // replace solverType |
2213 | | if (solverCharacteristics_->tryCuts()) { |
2214 | | |
2215 | | if (numberUnsatisfied) { |
2216 | | // User event |
2217 | | if (!eventHappened_ && feasible) |
2218 | | feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_, |
2219 | | NULL); |
2220 | | else |
2221 | | feasible = false; |
2222 | | } else if (solverCharacteristics_->solutionAddsCuts() || |
2223 | | solverCharacteristics_->alwaysTryCutsAtRootNode()) { |
2224 | | // may generate cuts and turn the solution |
2225 | | //to an infeasible one |
2226 | | feasible = solveWithCuts(cuts, 1, |
2227 | | NULL); |
2228 | | } |
2229 | | } |
2230 | | // check extra info on feasibility |
2231 | | if (!solverCharacteristics_->mipFeasible()) |
2232 | | feasible = false; |
2233 | | } |
2234 | | // make cut generators less aggressive |
2235 | | for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) { |
2236 | | CglCutGenerator * generator = generator_[iCutGenerator]->generator(); |
2237 | | generator->setAggressiveness(generator->getAggressiveness() - 100); |
2238 | | } |
2239 | | currentNumberCuts_ = numberNewCuts_ ; |
2240 | | // See if can stop on gap |
2241 | | bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense(); |
2242 | | testGap = CoinMax(dblParam_[CbcAllowableGap], |
2243 | | CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_)) |
2244 | | * dblParam_[CbcAllowableFractionGap]); |
2245 | | if (bestObjective_ - bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) { |
2246 | | if (bestPossibleObjective_ < getCutoff()) |
2247 | | stoppedOnGap_ = true ; |
2248 | | feasible = false; |
2249 | | } |
2250 | | // User event |
2251 | | if (eventHappened_) |
2252 | | feasible = false; |
2253 | | #if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX) |
2254 | | if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) { |
2255 | | // Use Cplex to do search! |
2256 | | double time1 = CoinCpuTime(); |
2257 | | OsiClpSolverInterface * clpSolver |
2258 | | = dynamic_cast<OsiClpSolverInterface *> (solver_); |
2259 | | OsiCpxSolverInterface cpxSolver; |
2260 | | double direction = clpSolver->getObjSense(); |
2261 | | cpxSolver.setObjSense(direction); |
2262 | | // load up cplex |
2263 | | const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol(); |
2264 | | const double * rowLower = continuousSolver_->getRowLower(); |
2265 | | const double * rowUpper = continuousSolver_->getRowUpper(); |
2266 | | const double * columnLower = continuousSolver_->getColLower(); |
2267 | | const double * columnUpper = continuousSolver_->getColUpper(); |
2268 | | const double * objective = continuousSolver_->getObjCoefficients(); |
2269 | | cpxSolver.loadProblem(*matrix, columnLower, columnUpper, |
2270 | | objective, rowLower, rowUpper); |
2271 | | double * setSol = new double [numberIntegers_]; |
2272 | | int * setVar = new int [numberIntegers_]; |
2273 | | // cplex doesn't know about objective offset |
2274 | | double offset = clpSolver->getModelPtr()->objectiveOffset(); |
2275 | | for (int i = 0; i < numberIntegers_; i++) { |
2276 | | int iColumn = integerVariable_[i]; |
2277 | | cpxSolver.setInteger(iColumn); |
2278 | | if (bestSolution_) { |
2279 | | setSol[i] = bestSolution_[iColumn]; |
2280 | | setVar[i] = iColumn; |
2281 | | } |
2282 | | } |
2283 | | CPXENVptr env = cpxSolver.getEnvironmentPtr(); |
2284 | | CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL); |
2285 | | cpxSolver.switchToMIP(); |
2286 | | if (bestSolution_) { |
2287 | | CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol); |
2288 | | } |
2289 | | if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) { |
2290 | | // add cuts |
2291 | | const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow(); |
2292 | | const double * rhs = clpSolver->getRightHandSide(); |
2293 | | const char * rowSense = clpSolver->getRowSense(); |
2294 | | const double * elementByRow = matrix->getElements(); |
2295 | | const int * column = matrix->getIndices(); |
2296 | | const CoinBigIndex * rowStart = matrix->getVectorStarts(); |
2297 | | const int * rowLength = matrix->getVectorLengths(); |
2298 | | int nStart = continuousSolver_->getNumRows(); |
2299 | | int nRows = clpSolver->getNumRows(); |
2300 | | int size = rowStart[nRows-1] + rowLength[nRows-1] - |
2301 | | rowStart[nStart]; |
2302 | | int nAdd = 0; |
2303 | | double * rmatval = new double [size]; |
2304 | | int * rmatind = new int [size]; |
2305 | | int * rmatbeg = new int [nRows-nStart+1]; |
2306 | | size = 0; |
2307 | | rmatbeg[0] = 0; |
2308 | | for (int i = nStart; i < nRows; i++) { |
2309 | | for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) { |
2310 | | rmatind[size] = column[k]; |
2311 | | rmatval[size++] = elementByRow[k]; |
2312 | | } |
2313 | | nAdd++; |
2314 | | rmatbeg[nAdd] = size; |
2315 | | } |
2316 | | CPXaddlazyconstraints(env, lpPtr, nAdd, size, |
2317 | | rhs, rowSense, rmatbeg, |
2318 | | rmatind, rmatval, NULL); |
2319 | | CPXsetintparam( env, CPX_PARAM_REDUCE, |
2320 | | // CPX_PREREDUCE_NOPRIMALORDUAL (0) |
2321 | | CPX_PREREDUCE_PRIMALONLY); |
2322 | | } |
2323 | | if (getCutoff() < 1.0e50) { |
2324 | | double useCutoff = getCutoff() + offset; |
2325 | | if (bestObjective_ < 1.0e50) |
2326 | | useCutoff = bestObjective_ + offset + 1.0e-7; |
2327 | | cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff* |
2328 | | direction); |
2329 | | if ( direction > 0.0 ) |
2330 | | CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min |
2331 | | else |
2332 | | CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max |
2333 | | } |
2334 | | CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]); |
2335 | | delete [] setSol; |
2336 | | delete [] setVar; |
2337 | | char printBuffer[200]; |
2338 | | if (offset) { |
2339 | | sprintf(printBuffer, "Add %g to all Cplex messages for true objective", |
2340 | | -offset); |
2341 | | messageHandler()->message(CBC_GENERAL, messages()) |
2342 | | << printBuffer << CoinMessageEol ; |
2343 | | cpxSolver.setDblParam(OsiObjOffset, offset); |
2344 | | } |
2345 | | cpxSolver.branchAndBound(); |
2346 | | double timeTaken = CoinCpuTime() - time1; |
2347 | | sprintf(printBuffer, "Cplex took %g seconds", |
2348 | | timeTaken); |
2349 | | messageHandler()->message(CBC_GENERAL, messages()) |
2350 | | << printBuffer << CoinMessageEol ; |
2351 | | numberExtraNodes_ = CPXgetnodecnt(env, lpPtr); |
2352 | | numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr); |
2353 | | double value = cpxSolver.getObjValue() * direction; |
2354 | | if (cpxSolver.isProvenOptimal() && value <= getCutoff()) { |
2355 | | feasible = true; |
2356 | | clpSolver->setWarmStart(NULL); |
2357 | | // try and do solution |
2358 | | double * newSolution = |
2359 | | CoinCopyOfArray(cpxSolver.getColSolution(), |
2360 | | getNumCols()); |
2361 | | setBestSolution(CBC_STRONGSOL, value, newSolution) ; |
2362 | | delete [] newSolution; |
2363 | | } |
2364 | | feasible = false; |
2365 | | } |
2366 | | #endif |
2367 | | if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) { |
2368 | | fastNodeDepth_ = -1; |
2369 | | CbcObject * obj = |
2370 | | new CbcFollowOn(this); |
2371 | | obj->setPriority(1); |
2372 | | addObjects(1, &obj); |
2373 | | } |
2374 | | int saveNumberSolves = numberSolves_; |
2375 | | int saveNumberIterations = numberIterations_; |
2376 | | if (fastNodeDepth_ >= 0 &&/*!parentModel_*/(specialOptions_&2048) == 0) { |
2377 | | // add in a general depth object doClp |
2378 | | int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100); |
2379 | | CbcObject * obj = |
2380 | | new CbcGeneralDepth(this, type); |
2381 | | addObjects(1, &obj); |
2382 | | // mark as done |
2383 | | fastNodeDepth_ += 1000000; |
2384 | | delete obj; |
2385 | | // fake number of objects |
2386 | | numberObjects_--; |
2387 | | if (parallelMode() < -1) { |
2388 | | // But make sure position is correct |
2389 | | OsiObject * obj2 = object_[numberObjects_]; |
2390 | | obj = dynamic_cast<CbcObject *> (obj2); |
2391 | | assert (obj); |
2392 | | obj->setPosition(numberObjects_); |
2393 | | } |
2394 | | } |
2395 | | #ifdef COIN_HAS_CLP |
2396 | | #ifdef NO_CRUNCH |
2397 | | if (true) { |
2398 | | OsiClpSolverInterface * clpSolver |
2399 | | = dynamic_cast<OsiClpSolverInterface *> (solver_); |
2400 | | if (clpSolver && !parentModel_) { |
2401 | | ClpSimplex * clpSimplex = clpSolver->getModelPtr(); |
2402 | | clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072); |
2403 | | //clpSimplex->startPermanentArrays(); |
2404 | | clpSimplex->setPersistenceFlag(2); |
2405 | | } |
2406 | | } |
2407 | | #endif |
2408 | | #endif |
2409 | | // Save copy of solver |
2410 | | OsiSolverInterface * saveSolver = NULL; |
2411 | | if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0) |
2412 | | saveSolver = solver_->clone(); |
2413 | | double checkCutoffForRestart = 1.0e100; |
2414 | | if (saveSolver && (specialOptions_&32768) != 0) { |
2415 | | // See if worth trying reduction |
2416 | | checkCutoffForRestart = getCutoff(); |
2417 | | bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() && |
2418 | | (checkCutoffForRestart < 1.0e20); |
| 1229 | } else if (numberObjects_ && object_) { |
| 1230 | numberOriginalObjects = numberObjects_; |
| 1231 | // redo sequence |
| 1232 | numberIntegers_ = 0; |
2420 | | if (tryNewSearch) { |
2421 | | #ifdef CLP_INVESTIGATE |
2422 | | printf("after %d nodes, cutoff %g - looking\n", |
2423 | | numberNodes_, getCutoff()); |
2424 | | #endif |
2425 | | saveSolver->resolve(); |
2426 | | double direction = saveSolver->getObjSense() ; |
2427 | | double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ; |
2428 | | double tolerance; |
2429 | | saveSolver->getDblParam(OsiDualTolerance, tolerance) ; |
2430 | | if (gap <= 0.0) |
2431 | | gap = tolerance; |
2432 | | gap += 100.0 * tolerance; |
2433 | | double integerTolerance = getDblParam(CbcIntegerTolerance) ; |
2434 | | |
2435 | | const double *lower = saveSolver->getColLower() ; |
2436 | | const double *upper = saveSolver->getColUpper() ; |
2437 | | const double *solution = saveSolver->getColSolution() ; |
2438 | | const double *reducedCost = saveSolver->getReducedCost() ; |
2439 | | |
2440 | | int numberFixed = 0 ; |
2441 | | int numberFixed2 = 0; |
2442 | | for (int i = 0 ; i < numberIntegers_ ; i++) { |
2443 | | int iColumn = integerVariable_[i] ; |
2444 | | double djValue = direction * reducedCost[iColumn] ; |
2445 | | if (upper[iColumn] - lower[iColumn] > integerTolerance) { |
2446 | | if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) { |
2447 | | saveSolver->setColUpper(iColumn, lower[iColumn]) ; |
2448 | | numberFixed++ ; |
2449 | | } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) { |
2450 | | saveSolver->setColLower(iColumn, upper[iColumn]) ; |
2451 | | numberFixed++ ; |
2452 | | } |
2453 | | } else { |
2454 | | numberFixed2++; |
2455 | | } |
2456 | | } |
2457 | | #ifdef COIN_DEVELOP |
2458 | | if ((specialOptions_&1) != 0) { |
2459 | | const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ; |
2460 | | if (debugger) { |
2461 | | printf("Contains optimal\n") ; |
2462 | | OsiSolverInterface * temp = saveSolver->clone(); |
2463 | | const double * solution = debugger->optimalSolution(); |
2464 | | const double *lower = temp->getColLower() ; |
2465 | | const double *upper = temp->getColUpper() ; |
2466 | | int n = temp->getNumCols(); |
2467 | | for (int i = 0; i < n; i++) { |
2468 | | if (temp->isInteger(i)) { |
2469 | | double value = floor(solution[i] + 0.5); |
2470 | | assert (value >= lower[i] && value <= upper[i]); |
2471 | | temp->setColLower(i, value); |
2472 | | temp->setColUpper(i, value); |
2473 | | } |
2474 | | } |
2475 | | temp->writeMps("reduced_fix"); |
2476 | | delete temp; |
2477 | | saveSolver->writeMps("reduced"); |
2478 | | } else { |
2479 | | abort(); |
2480 | | } |
2481 | | } |
2482 | | printf("Restart could fix %d integers (%d already fixed)\n", |
2483 | | numberFixed + numberFixed2, numberFixed2); |
2484 | | #endif |
2485 | | numberFixed += numberFixed2; |
2486 | | if (numberFixed*20 < numberColumns) |
2487 | | tryNewSearch = false; |
2488 | | } |
2489 | | if (tryNewSearch) { |
2490 | | // back to solver without cuts? |
2491 | | OsiSolverInterface * solver2 = continuousSolver_->clone(); |
2492 | | const double *lower = saveSolver->getColLower() ; |
2493 | | const double *upper = saveSolver->getColUpper() ; |
2494 | | for (int i = 0 ; i < numberIntegers_ ; i++) { |
2495 | | int iColumn = integerVariable_[i] ; |
2496 | | solver2->setColLower(iColumn, lower[iColumn]); |
2497 | | solver2->setColUpper(iColumn, upper[iColumn]); |
2498 | | } |
2499 | | // swap |
2500 | | delete saveSolver; |
2501 | | saveSolver = solver2; |
2502 | | double * newSolution = new double[numberColumns]; |
2503 | | double objectiveValue = checkCutoffForRestart; |
2504 | | CbcSerendipity heuristic(*this); |
2505 | | if (bestSolution_) |
2506 | | heuristic.setInputSolution(bestSolution_, bestObjective_); |
2507 | | heuristic.setFractionSmall(0.5); |
2508 | | heuristic.setFeasibilityPumpOptions(1008013); |
2509 | | // Use numberNodes to say how many are original rows |
2510 | | heuristic.setNumberNodes(continuousSolver_->getNumRows()); |
2511 | | #ifdef COIN_DEVELOP |
2512 | | if (continuousSolver_->getNumRows() < |
2513 | | saveSolver->getNumRows()) |
2514 | | printf("%d rows added ZZZZZ\n", |
2515 | | solver_->getNumRows() - continuousSolver_->getNumRows()); |
2516 | | #endif |
2517 | | int returnCode = heuristic.smallBranchAndBound(saveSolver, |
2518 | | -1, newSolution, |
2519 | | objectiveValue, |
2520 | | checkCutoffForRestart, "Reduce"); |
2521 | | if (returnCode < 0) { |
2522 | | #ifdef COIN_DEVELOP |
2523 | | printf("Restart - not small enough to do search after fixing\n"); |
2524 | | #endif |
2525 | | delete [] newSolution; |
2526 | | } else { |
2527 | | if ((returnCode&1) != 0) { |
2528 | | // increment number of solutions so other heuristics can test |
2529 | | numberSolutions_++; |
2530 | | numberHeuristicSolutions_++; |
2531 | | lastHeuristic_ = NULL; |
2532 | | setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ; |
2533 | | } |
2534 | | delete [] newSolution; |
2535 | | feasible = false; // stop search |
2536 | | } |
2537 | | } |
2538 | | } |
2539 | | /* |
2540 | | We've taken the continuous relaxation as far as we can. Time to branch. |
2541 | | The first order of business is to actually create a node. chooseBranch |
2542 | | currently uses strong branching to evaluate branch object candidates, |
2543 | | unless forced back to simple branching. If chooseBranch concludes that a |
2544 | | branching candidate is monotone (anyAction == -1) or infeasible (anyAction |
2545 | | == -2) when forced to integer values, it returns here immediately. |
2546 | | |
2547 | | Monotone variables trigger a call to resolve(). If the problem remains |
2548 | | feasible, try again to choose a branching variable. At the end of the loop, |
2549 | | resolved == true indicates that some variables were fixed. |
2550 | | |
2551 | | Loss of feasibility will result in the deletion of newNode. |
2552 | | */ |
2553 | | |
2554 | | bool resolved = false ; |
2555 | | CbcNode *newNode = NULL ; |
2556 | | numberFixedAtRoot_ = 0; |
2557 | | numberFixedNow_ = 0; |
2558 | | int numberIterationsAtContinuous = numberIterations_; |
2559 | | //solverCharacteristics_->setSolver(solver_); |
2560 | | if (feasible) { |
2561 | | //#define HOTSTART -1 |
2562 | | #if HOTSTART<0 |
2563 | | if (bestSolution_ && !parentModel_ && !hotstartSolution_ && |
2564 | | (moreSpecialOptions_&1024) != 0) { |
2565 | | // Set priorities so only branch on ones we need to |
2566 | | // use djs and see if only few branches needed |
2567 | | #ifndef NDEBUG |
2568 | | double integerTolerance = getIntegerTolerance() ; |
2569 | | #endif |
2570 | | bool possible = true; |
2571 | | const double * saveLower = continuousSolver_->getColLower(); |
2572 | | const double * saveUpper = continuousSolver_->getColUpper(); |
2573 | | for (int i = 0; i < numberObjects_; i++) { |
2574 | | const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]); |
2575 | | if (thisOne) { |
2576 | | int iColumn = thisOne->columnNumber(); |
2577 | | if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) { |
2578 | | possible = false; |
2579 | | break; |
2580 | | } |
2581 | | } else { |
2582 | | possible = false; |
2583 | | break; |
2584 | | } |
2585 | | } |
2586 | | if (possible) { |
2587 | | OsiSolverInterface * solver = continuousSolver_->clone(); |
2588 | | int numberColumns = solver->getNumCols(); |
2589 | | for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) { |
2590 | | double value = bestSolution_[iColumn] ; |
2591 | | value = CoinMax(value, saveLower[iColumn]) ; |
2592 | | value = CoinMin(value, saveUpper[iColumn]) ; |
2593 | | value = floor(value + 0.5); |
2594 | | if (solver->isInteger(iColumn)) { |
2595 | | solver->setColLower(iColumn, value); |
2596 | | solver->setColUpper(iColumn, value); |
2597 | | } |
2598 | | } |
2599 | | solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry); |
2600 | | // objlim and all slack |
2601 | | double direction = solver->getObjSense(); |
2602 | | solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction); |
2603 | | CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart()); |
2604 | | solver->setWarmStart(basis); |
2605 | | delete basis; |
2606 | | bool changed = true; |
2607 | | hotstartPriorities_ = new int [numberColumns]; |
2608 | | for (int iColumn = 0; iColumn < numberColumns; iColumn++) |
2609 | | hotstartPriorities_[iColumn] = 1; |
2610 | | while (changed) { |
2611 | | changed = false; |
2612 | | solver->resolve(); |
2613 | | if (!solver->isProvenOptimal()) { |
2614 | | possible = false; |
2615 | | break; |
2616 | | } |
2617 | | const double * dj = solver->getReducedCost(); |
2618 | | const double * colLower = solver->getColLower(); |
2619 | | const double * colUpper = solver->getColUpper(); |
2620 | | const double * solution = solver->getColSolution(); |
2621 | | int nAtLbNatural = 0; |
2622 | | int nAtUbNatural = 0; |
2623 | | int nZeroDj = 0; |
2624 | | int nForced = 0; |
2625 | | for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) { |
2626 | | double value = solution[iColumn] ; |
2627 | | value = CoinMax(value, saveLower[iColumn]) ; |
2628 | | value = CoinMin(value, saveUpper[iColumn]) ; |
2629 | | if (solver->isInteger(iColumn)) { |
2630 | | assert(fabs(value - solution[iColumn]) <= integerTolerance) ; |
2631 | | if (hotstartPriorities_[iColumn] == 1) { |
2632 | | if (dj[iColumn] < -1.0e-6) { |
2633 | | // negative dj |
2634 | | if (saveUpper[iColumn] == colUpper[iColumn]) { |
2635 | | nAtUbNatural++; |
2636 | | hotstartPriorities_[iColumn] = 2; |
2637 | | solver->setColLower(iColumn, saveLower[iColumn]); |
2638 | | solver->setColUpper(iColumn, saveUpper[iColumn]); |
2639 | | } else { |
2640 | | nForced++; |
2641 | | } |
2642 | | } else if (dj[iColumn] > 1.0e-6) { |
2643 | | // positive dj |
2644 | | if (saveLower[iColumn] == colLower[iColumn]) { |
2645 | | nAtLbNatural++; |
2646 | | hotstartPriorities_[iColumn] = 2; |
2647 | | solver->setColLower(iColumn, saveLower[iColumn]); |
2648 | | solver->setColUpper(iColumn, saveUpper[iColumn]); |
2649 | | } else { |
2650 | | nForced++; |
2651 | | } |
2652 | | } else { |
2653 | | // zero dj |
2654 | | nZeroDj++; |
2655 | | } |
2656 | | } |
2657 | | } |
2658 | | } |
2659 | | #ifdef CLP_INVESTIGATE |
2660 | | printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n", |
2661 | | nForced, nAtLbNatural, nAtUbNatural, nZeroDj); |
2662 | | #endif |
2663 | | if (nAtLbNatural || nAtUbNatural) { |
2664 | | changed = true; |
2665 | | } else { |
2666 | | if (nForced + nZeroDj > 50 || |
2667 | | (nForced + nZeroDj)*4 > numberIntegers_) |
2668 | | possible = false; |
2669 | | } |
2670 | | } |
2671 | | delete solver; |
2672 | | } |
2673 | | if (possible) { |
2674 | | setHotstartSolution(bestSolution_); |
2675 | | if (!saveCompare) { |
2676 | | // create depth first comparison |
2677 | | saveCompare = nodeCompare_; |
2678 | | // depth first |
2679 | | nodeCompare_ = new CbcCompareDepth(); |
2680 | | tree_->setComparison(*nodeCompare_) ; |
2681 | | } |
2682 | | } else { |
2683 | | delete [] hotstartPriorities_; |
2684 | | hotstartPriorities_ = NULL; |
2685 | | } |
2686 | | } |
2687 | | #endif |
2688 | | #if HOTSTART>0 |
2689 | | if (hotstartSolution_ && !hotstartPriorities_) { |
2690 | | // Set up hot start |
2691 | | OsiSolverInterface * solver = solver_->clone(); |
2692 | | double direction = solver_->getObjSense() ; |
2693 | | int numberColumns = solver->getNumCols(); |
2694 | | double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns); |
2695 | | double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns); |
2696 | | // move solution |
2697 | | solver->setColSolution(hotstartSolution_); |
2698 | | // point to useful information |
2699 | | const double * saveSolution = testSolution_; |
2700 | | testSolution_ = solver->getColSolution(); |
2701 | | OsiBranchingInformation usefulInfo = usefulInformation(); |
2702 | | testSolution_ = saveSolution; |
2703 | | /* |
2704 | | Run through the objects and use feasibleRegion() to set variable bounds |
2705 | | so as to fix the variables specified in the objects at their value in this |
2706 | | solution. Since the object list contains (at least) one object for every |
2707 | | integer variable, this has the effect of fixing all integer variables. |
2708 | | */ |
2709 | | for (int i = 0; i < numberObjects_; i++) |
2710 | | object_[i]->feasibleRegion(solver, &usefulInfo); |
2711 | | solver->resolve(); |
2712 | | assert (solver->isProvenOptimal()); |
2713 | | double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ; |
2714 | | const double * dj = solver->getReducedCost(); |
2715 | | const double * colLower = solver->getColLower(); |
2716 | | const double * colUpper = solver->getColUpper(); |
2717 | | const double * solution = solver->getColSolution(); |
2718 | | int nAtLbNatural = 0; |
2719 | | int nAtUbNatural = 0; |
2720 | | int nAtLbNaturalZero = 0; |
2721 | | int nAtUbNaturalZero = 0; |
2722 | | int nAtLbFixed = 0; |
2723 | | int nAtUbFixed = 0; |
2724 | | int nAtOther = 0; |
2725 | | int nAtOtherNatural = 0; |
2726 | | int nNotNeeded = 0; |
2727 | | delete [] hotstartSolution_; |
2728 | | hotstartSolution_ = new double [numberColumns]; |
2729 | | delete [] hotstartPriorities_; |
2730 | | hotstartPriorities_ = new int [numberColumns]; |
2731 | | int * order = (int *) saveUpper; |
2732 | | int nFix = 0; |
2733 | | double bestRatio = COIN_DBL_MAX; |
2734 | | for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) { |
2735 | | double value = solution[iColumn] ; |
2736 | | value = CoinMax(value, saveLower[iColumn]) ; |
2737 | | value = CoinMin(value, saveUpper[iColumn]) ; |
2738 | | double sortValue = COIN_DBL_MAX; |
2739 | | if (solver->isInteger(iColumn)) { |
2740 | | assert(fabs(value - solution[iColumn]) <= 1.0e-5) ; |
2741 | | double value2 = floor(value + 0.5); |
2742 | | if (dj[iColumn] < -1.0e-6) { |
2743 | | // negative dj |
2744 | | //assert (value2==colUpper[iColumn]); |
2745 | | if (saveUpper[iColumn] == colUpper[iColumn]) { |
2746 | | nAtUbNatural++; |
2747 | | sortValue = 0.0; |
2748 | | double value = -dj[iColumn]; |
2749 | | if (value > gap) |
2750 | | nFix++; |
2751 | | else if (gap < value*bestRatio) |
2752 | | bestRatio = gap / value; |
2753 | | if (saveLower[iColumn] != colLower[iColumn]) { |
2754 | | nNotNeeded++; |
2755 | | sortValue = 1.0e20; |
2756 | | } |
2757 | | } else if (saveLower[iColumn] == colUpper[iColumn]) { |
2758 | | nAtLbFixed++; |
2759 | | sortValue = dj[iColumn]; |
2760 | | } else { |
2761 | | nAtOther++; |
2762 | | sortValue = 0.0; |
2763 | | if (saveLower[iColumn] != colLower[iColumn] && |
2764 | | saveUpper[iColumn] != colUpper[iColumn]) { |
2765 | | nNotNeeded++; |
2766 | | sortValue = 1.0e20; |
2767 | | } |
2768 | | } |
2769 | | } else if (dj[iColumn] > 1.0e-6) { |
2770 | | // positive dj |
2771 | | //assert (value2==colLower[iColumn]); |
2772 | | if (saveLower[iColumn] == colLower[iColumn]) { |
2773 | | nAtLbNatural++; |
2774 | | sortValue = 0.0; |
2775 | | double value = dj[iColumn]; |
2776 | | if (value > gap) |
2777 | | nFix++; |
2778 | | else if (gap < value*bestRatio) |
2779 | | bestRatio = gap / value; |
2780 | | if (saveUpper[iColumn] != colUpper[iColumn]) { |
2781 | | nNotNeeded++; |
2782 | | sortValue = 1.0e20; |
2783 | | } |
2784 | | } else if (saveUpper[iColumn] == colLower[iColumn]) { |
2785 | | nAtUbFixed++; |
2786 | | sortValue = -dj[iColumn]; |
2787 | | } else { |
2788 | | nAtOther++; |
2789 | | sortValue = 0.0; |
2790 | | if (saveLower[iColumn] != colLower[iColumn] && |
2791 | | saveUpper[iColumn] != colUpper[iColumn]) { |
2792 | | nNotNeeded++; |
2793 | | sortValue = 1.0e20; |
2794 | | } |
2795 | | } |
2796 | | } else { |
2797 | | // zero dj |
2798 | | if (value2 == saveUpper[iColumn]) { |
2799 | | nAtUbNaturalZero++; |
2800 | | sortValue = 0.0; |
2801 | | if (saveLower[iColumn] != colLower[iColumn]) { |
2802 | | nNotNeeded++; |
2803 | | sortValue = 1.0e20; |
2804 | | } |
2805 | | } else if (value2 == saveLower[iColumn]) { |
2806 | | nAtLbNaturalZero++; |
2807 | | sortValue = 0.0; |
2808 | | } else { |
2809 | | nAtOtherNatural++; |
2810 | | sortValue = 0.0; |
2811 | | if (saveLower[iColumn] != colLower[iColumn] && |
2812 | | saveUpper[iColumn] != colUpper[iColumn]) { |
2813 | | nNotNeeded++; |
2814 | | sortValue = 1.0e20; |
2815 | | } |
2816 | | } |
2817 | | } |
2818 | | #if HOTSTART==3 |
2819 | | sortValue = -fabs(dj[iColumn]); |
2820 | | #endif |
2821 | | } |
2822 | | hotstartSolution_[iColumn] = value ; |
2823 | | saveLower[iColumn] = sortValue; |
2824 | | order[iColumn] = iColumn; |
2825 | | } |
2826 | | printf("** can fix %d columns - best ratio for others is %g on gap of %g\n", |
2827 | | nFix, bestRatio, gap); |
2828 | | int nNeg = 0; |
2829 | | CoinSort_2(saveLower, saveLower + numberColumns, order); |
2830 | | for (int i = 0; i < numberColumns; i++) { |
2831 | | if (saveLower[i] < 0.0) { |
2832 | | nNeg++; |
2833 | | #if HOTSTART==2||HOTSTART==3 |
2834 | | // swap sign ? |
2835 | | saveLower[i] = -saveLower[i]; |
2836 | | #endif |
2837 | | } |
2838 | | } |
2839 | | CoinSort_2(saveLower, saveLower + nNeg, order); |
2840 | | for (int i = 0; i < numberColumns; i++) { |
2841 | | #if HOTSTART==1 |
2842 | | hotstartPriorities_[order[i]] = 100; |
2843 | | #else |
2844 | | hotstartPriorities_[order[i]] = -(i + 1); |
2845 | | #endif |
2846 | | } |
2847 | | printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n", |
2848 | | nAtLbNatural, |
2849 | | nAtUbNatural, |
2850 | | nAtLbNaturalZero, |
2851 | | nAtUbNaturalZero, |
2852 | | nAtLbFixed, |
2853 | | nAtUbFixed, |
2854 | | nAtOther, |
2855 | | nAtOtherNatural, nNotNeeded, nNeg); |
2856 | | delete [] saveLower; |
2857 | | delete [] saveUpper; |
2858 | | if (!saveCompare) { |
2859 | | // create depth first comparison |
2860 | | saveCompare = nodeCompare_; |
2861 | | // depth first |
2862 | | nodeCompare_ = new CbcCompareDepth(); |
2863 | | tree_->setComparison(*nodeCompare_) ; |
2864 | | } |
2865 | | } |
2866 | | #endif |
2867 | | if (probingInfo_) { |
2868 | | int number01 = probingInfo_->numberIntegers(); |
2869 | | //const fixEntry * entry = probingInfo_->fixEntries(); |
2870 | | const int * toZero = probingInfo_->toZero(); |
2871 | | //const int * toOne = probingInfo_->toOne(); |
2872 | | //const int * integerVariable = probingInfo_->integerVariable(); |
2873 | | if (toZero[number01]) { |
2874 | | if (probingInfo_->packDown()) { |
2875 | | #ifdef CLP_INVESTIGATE |
2876 | | printf("%d implications on %d 0-1\n", toZero[number01], number01); |
2877 | | #endif |
2878 | | CglImplication implication(probingInfo_); |
2879 | | addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200); |
2880 | | } else { |
2881 | | delete probingInfo_; |
2882 | | probingInfo_ = NULL; |
2883 | | } |
2884 | | } else { |
2885 | | delete probingInfo_; |
2886 | | |
2887 | | probingInfo_ = NULL; |
2888 | | } |
2889 | | } |
2890 | | newNode = new CbcNode ; |
2891 | | // Set objective value (not so obvious if NLP etc) |
2892 | | setObjectiveValue(newNode, NULL); |
2893 | | anyAction = -1 ; |
2894 | | // To make depth available we may need a fake node |
2895 | | CbcNode fakeNode; |
2896 | | if (!currentNode_) { |
2897 | | // Not true if sub trees assert (!numberNodes_); |
2898 | | currentNode_ = &fakeNode; |
2899 | | } |
2900 | | phase_ = 3; |
2901 | | // only allow 1000 passes |
2902 | | int numberPassesLeft = 1000; |
2903 | | // This is first crude step |
2904 | | if (numberAnalyzeIterations_) { |
2905 | | delete [] analyzeResults_; |
2906 | | analyzeResults_ = new double [4*numberIntegers_]; |
2907 | | numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_); |
2908 | | if (numberFixedAtRoot_ > 0) { |
2909 | | printf("%d fixed by analysis\n", numberFixedAtRoot_); |
2910 | | setPointers(solver_); |
2911 | | numberFixedNow_ = numberFixedAtRoot_; |
2912 | | } else if (numberFixedAtRoot_ < 0) { |
2913 | | printf("analysis found to be infeasible\n"); |
2914 | | anyAction = -2; |
2915 | | delete newNode ; |
2916 | | newNode = NULL ; |
2917 | | feasible = false ; |
2918 | | } |
2919 | | } |
2920 | | OsiSolverBranch * branches = NULL; |
2921 | | anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved, |
2922 | | NULL, NULL, NULL, branches); |
2923 | | if (anyAction == -2 || newNode->objectiveValue() >= cutoff) { |
2924 | | if (anyAction != -2) { |
2925 | | // zap parent nodeInfo |
2926 | | #ifdef COIN_DEVELOP |
2927 | | printf("zapping CbcNodeInfo %x\n", reinterpret_cast<int>(newNode->nodeInfo()->parent())); |
2928 | | #endif |
2929 | | if (newNode->nodeInfo()) |
2930 | | newNode->nodeInfo()->nullParent(); |
2931 | | } |
2932 | | delete newNode ; |
2933 | | newNode = NULL ; |
2934 | | feasible = false ; |
2935 | | } |
2936 | | } |
2937 | | /* |
2938 | | At this point, the root subproblem is infeasible or fathomed by bound |
2939 | | (newNode == NULL), or we're live with an objective value that satisfies the |
2940 | | current objective cutoff. |
2941 | | */ |
2942 | | assert (!newNode || newNode->objectiveValue() <= cutoff) ; |
2943 | | // Save address of root node as we don't want to delete it |
2944 | | // initialize for print out |
2945 | | int lastDepth = 0; |
2946 | | int lastUnsatisfied = 0; |
2947 | | if (newNode) |
2948 | | lastUnsatisfied = newNode->numberUnsatisfied(); |
2949 | | /* |
2950 | | The common case is that the lp relaxation is feasible but doesn't satisfy |
2951 | | integrality (i.e., newNode->branchingObject(), indicating we've been able to |
2952 | | select a branching variable). Remove any cuts that have gone slack due to |
2953 | | forcing monotone variables. Then tack on an CbcFullNodeInfo object and full |
2954 | | basis (via createInfo()) and stash the new cuts in the nodeInfo (via |
2955 | | addCuts()). If, by some miracle, we have an integral solution at the root |
2956 | | (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds |
2957 | | a valid solution for use by setBestSolution(). |
2958 | | */ |
2959 | | CoinWarmStartBasis *lastws = NULL ; |
2960 | | if (feasible && newNode->branchingObject()) { |
2961 | | if (resolved) { |
2962 | | takeOffCuts(cuts, false, NULL) ; |
2963 | | # ifdef CHECK_CUT_COUNTS |
2964 | | { printf("Number of rows after chooseBranch fix (root)" |
2965 | | "(active only) %d\n", |
2966 | | numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ; |
2967 | | const CoinWarmStartBasis* debugws = |
2968 | | dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ; |
2969 | | debugws->print() ; |
2970 | | delete debugws ; |
2971 | | } |
2972 | | # endif |
2973 | | } |
2974 | | //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ; |
2975 | | //newNode->nodeInfo()->addCuts(cuts, |
2976 | | // newNode->numberBranches(),whichGenerator_) ; |
2977 | | if (lastws) delete lastws ; |
2978 | | lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ; |
2979 | | } |
2980 | | /* |
2981 | | Continuous data to be used later |
2982 | | */ |
2983 | | continuousObjective_ = solver_->getObjValue() * solver_->getObjSense(); |
2984 | | continuousInfeasibilities_ = 0 ; |
2985 | | if (newNode) { |
2986 | | continuousObjective_ = newNode->objectiveValue() ; |
2987 | | delete [] continuousSolution_; |
2988 | | continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(), |
2989 | | numberColumns); |
2990 | | continuousInfeasibilities_ = newNode->numberUnsatisfied() ; |
2991 | | } |
2992 | | /* |
2993 | | Bound may have changed so reset in objects |
2994 | | */ |
2995 | | { |
2996 | | int i ; |
2997 | | for (i = 0; i < numberObjects_; i++) |
2998 | | object_[i]->resetBounds(solver_) ; |
2999 | | } |
3000 | | /* |
3001 | | Feasible? Then we should have either a live node prepped for future |
3002 | | expansion (indicated by variable() >= 0), or (miracle of miracles) an |
3003 | | integral solution at the root node. |
3004 | | |
3005 | | initializeInfo sets the reference counts in the nodeInfo object. Since |
3006 | | this node is still live, push it onto the heap that holds the live set. |
3007 | | */ |
3008 | | double bestValue = 0.0 ; |
3009 | | if (newNode) { |
3010 | | bestValue = newNode->objectiveValue(); |
3011 | | if (newNode->branchingObject()) { |
3012 | | newNode->initializeInfo() ; |
3013 | | tree_->push(newNode) ; |
3014 | | if (statistics_) { |
3015 | | if (numberNodes2_ == maximumStatistics_) { |
3016 | | maximumStatistics_ = 2 * maximumStatistics_; |
3017 | | CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_]; |
3018 | | memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *)); |
3019 | | memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *)); |
3020 | | delete [] statistics_; |
3021 | | statistics_ = temp; |
3022 | | } |
3023 | | assert (!statistics_[numberNodes2_]); |
3024 | | statistics_[numberNodes2_] = new CbcStatistics(newNode, this); |
3025 | | } |
3026 | | numberNodes2_++; |
3027 | | # ifdef CHECK_NODE |
3028 | | printf("Node %x on tree\n", newNode) ; |
3029 | | # endif |
3030 | | } else { |
3031 | | // continuous is integer |
3032 | | double objectiveValue = newNode->objectiveValue(); |
3033 | | setBestSolution(CBC_SOLUTION, objectiveValue, |
3034 | | solver_->getColSolution()) ; |
3035 | | delete newNode ; |
3036 | | newNode = NULL ; |
3037 | | } |
3038 | | } |
3039 | | |
3040 | | if (printFrequency_ <= 0) { |
3041 | | printFrequency_ = 1000 ; |
3042 | | if (getNumCols() > 2000) |
3043 | | printFrequency_ = 100 ; |
3044 | | } |
3045 | | /* |
3046 | | It is possible that strong branching fixes one variable and then the code goes round |
3047 | | again and again. This can take too long. So we need to warn user - just once. |
3048 | | */ |
3049 | | numberLongStrong_ = 0; |
3050 | | CbcNode * createdNode = NULL; |
3051 | | #ifdef CBC_THREAD |
3052 | | CbcModel ** threadModel = NULL; |
3053 | | Coin_pthread_t * threadId = NULL; |
3054 | | int * threadCount = NULL; |
3055 | | pthread_mutex_t mutex; |
3056 | | pthread_cond_t condition_main; |
3057 | | pthread_mutex_t condition_mutex; |
3058 | | pthread_mutex_t * mutex2 = NULL; |
3059 | | pthread_cond_t * condition2 = NULL; |
3060 | | threadStruct * threadInfo = NULL; |
3061 | | bool locked = false; |
3062 | | int threadStats[6]; |
3063 | | int defaultParallelIterations = 400; |
3064 | | int defaultParallelNodes = 2; |
3065 | | memset(threadStats, 0, sizeof(threadStats)); |
3066 | | double timeWaiting = 0.0; |
3067 | | // For now just one model |
3068 | | if (numberThreads_) { |
3069 | | nodeCompare_->sayThreaded(); // need to use addresses |
3070 | | threadId = new Coin_pthread_t [numberThreads_]; |
3071 | | threadCount = new int [numberThreads_]; |
3072 | | CoinZeroN(threadCount, numberThreads_); |
3073 | | pthread_mutex_init(&mutex, NULL); |
3074 | | pthread_cond_init(&condition_main, NULL); |
3075 | | pthread_mutex_init(&condition_mutex, NULL); |
3076 | | threadModel = new CbcModel * [numberThreads_+1]; |
3077 | | threadInfo = new threadStruct [numberThreads_+1]; |
3078 | | mutex2 = new pthread_mutex_t [numberThreads_]; |
3079 | | condition2 = new pthread_cond_t [numberThreads_]; |
3080 | | if (parallelMode() < -1) { |
3081 | | // May need for deterministic |
3082 | | saveObjects = new OsiObject * [numberObjects_]; |
3083 | | for (int i = 0; i < numberObjects_; i++) { |
3084 | | saveObjects[i] = object_[i]->clone(); |
3085 | | } |
3086 | | } |
3087 | | // we don't want a strategy object |
3088 | | CbcStrategy * saveStrategy = strategy_; |
3089 | | strategy_ = NULL; |
3090 | | for (int i = 0; i < numberThreads_; i++) { |
3091 | | pthread_mutex_init(mutex2 + i, NULL); |
3092 | | pthread_cond_init(condition2 + i, NULL); |
3093 | | threadId[i].status = 0; |
3094 | | threadInfo[i].baseModel = this; |
3095 | | threadModel[i] = new CbcModel(*this, true); |
3096 | | threadModel[i]->synchronizeHandlers(1); |
3097 | | #ifdef COIN_HAS_CLP |
3098 | | // Solver may need to know about model |
3099 | | CbcModel * thisModel = threadModel[i]; |
3100 | | CbcOsiSolver * solver = |
3101 | | dynamic_cast<CbcOsiSolver *>(thisModel->solver()) ; |
3102 | | if (solver) |
3103 | | solver->setCbcModel(thisModel); |
3104 | | #endif |
3105 | | mutex_ = reinterpret_cast<void *> (threadInfo + i); |
3106 | | threadModel[i]->moveToModel(this, -1); |
3107 | | threadInfo[i].thisModel = threadModel[i]; |
3108 | | threadInfo[i].node = NULL; |
3109 | | threadInfo[i].createdNode = NULL; |
3110 | | threadInfo[i].threadIdOfBase.thr = pthread_self(); |
3111 | | threadInfo[i].mutex = &mutex; |
3112 | | threadInfo[i].mutex2 = mutex2 + i; |
3113 | | threadInfo[i].condition2 = condition2 + i; |
3114 | | threadInfo[i].returnCode = -1; |
3115 | | threadInfo[i].timeLocked = 0.0; |
3116 | | threadInfo[i].timeWaitingToLock = 0.0; |
3117 | | threadInfo[i].timeWaitingToStart = 0.0; |
3118 | | threadInfo[i].timeInThread = 0.0; |
3119 | | threadInfo[i].numberTimesLocked = 0; |
3120 | | threadInfo[i].numberTimesUnlocked = 0; |
3121 | | threadInfo[i].numberTimesWaitingToStart = 0; |
3122 | | threadInfo[i].dantzigState = 0; // 0 unset, -1 waiting to be set, 1 set |
3123 | | threadInfo[i].locked = false; |
3124 | | threadInfo[i].delNode = NULL; |
3125 | | threadInfo[i].maxDeleteNode = 0; |
3126 | | threadInfo[i].nDeleteNode = 0; |
3127 | | threadInfo[i].nodesThisTime = 0; |
3128 | | threadInfo[i].iterationsThisTime = 0; |
3129 | | pthread_create(&(threadId[i].thr), NULL, doNodesThread, threadInfo + i); |
3130 | | threadId[i].status = 1; |
3131 | | } |
3132 | | strategy_ = saveStrategy; |
3133 | | // Do a partial one for base model |
3134 | | threadInfo[numberThreads_].baseModel = this; |
3135 | | threadModel[numberThreads_] = this; |
3136 | | mutex_ = reinterpret_cast<void *> (threadInfo + numberThreads_); |
3137 | | threadInfo[numberThreads_].node = NULL; |
3138 | | threadInfo[numberThreads_].mutex = &mutex; |
3139 | | threadInfo[numberThreads_].condition2 = &condition_main; |
3140 | | threadInfo[numberThreads_].mutex2 = &condition_mutex; |
3141 | | threadInfo[numberThreads_].timeLocked = 0.0; |
3142 | | threadInfo[numberThreads_].timeWaitingToLock = 0.0; |
3143 | | threadInfo[numberThreads_].numberTimesLocked = 0; |
3144 | | threadInfo[numberThreads_].numberTimesUnlocked = 0; |
3145 | | threadInfo[numberThreads_].locked = false; |
3146 | | } |
3147 | | #endif |
3148 | | #ifdef COIN_HAS_CLP |
3149 | | { |
3150 | | OsiClpSolverInterface * clpSolver |
3151 | | = dynamic_cast<OsiClpSolverInterface *> (solver_); |
3152 | | if (clpSolver && !parentModel_) { |
3153 | | clpSolver->computeLargestAway(); |
3154 | | } |
3155 | | } |
3156 | | #endif |
3157 | | /* |
3158 | | At last, the actual branch-and-cut search loop, which will iterate until |
3159 | | the live set is empty or we hit some limit (integrality gap, time, node |
3160 | | count, etc.). The overall flow is to rebuild a subproblem, reoptimise using |
3161 | | solveWithCuts(), choose a branching pattern with chooseBranch(), and finally |
3162 | | add the node to the live set. |
3163 | | |
3164 | | The first action is to winnow the live set to remove nodes which are worse |
3165 | | than the current objective cutoff. |
3166 | | */ |
3167 | | if (solver_->getRowCutDebuggerAlways()) { |
3168 | | OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways()); |
3169 | | const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ; |
3170 | | if (!debugger) { |
3171 | | // infeasible!! |
3172 | | printf("before search\n"); |
3173 | | const double * lower = solver_->getColLower(); |
3174 | | const double * upper = solver_->getColUpper(); |
3175 | | const double * solution = debuggerX->optimalSolution(); |
3176 | | int numberColumns = solver_->getNumCols(); |
3177 | | for (int i = 0; i < numberColumns; i++) { |
3178 | | if (solver_->isInteger(i)) { |
3179 | | if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6) |
3180 | | printf("**** "); |
3181 | | printf("%d %g <= %g <= %g\n", |
3182 | | i, lower[i], solution[i], upper[i]); |
3183 | | } |
3184 | | } |
3185 | | //abort(); |
3186 | | } |
3187 | | } |
3188 | | { |
3189 | | // may be able to change cutoff now |
3190 | | double cutoff = getCutoff(); |
3191 | | double increment = getDblParam(CbcModel::CbcCutoffIncrement) ; |
3192 | | if (cutoff > bestObjective_ - increment) { |
3193 | | cutoff = bestObjective_ - increment ; |
3194 | | setCutoff(cutoff) ; |
3195 | | } |
3196 | | } |
3197 | | #ifdef CBC_THREAD |
3198 | | bool goneParallel = false; |
3199 | | #endif |
3200 | | #define MAX_DEL_NODE 1 |
3201 | | CbcNode * delNode[MAX_DEL_NODE+1]; |
3202 | | int nDeleteNode = 0; |
3203 | | // For Printing etc when parallel |
3204 | | int lastEvery1000 = 0; |
3205 | | int lastPrintEvery = 0; |
3206 | | while (true) { |
3207 | | #ifdef CBC_THREAD |
3208 | | if (parallelMode() > 0 && !locked) { |
3209 | | lockThread(); |
3210 | | locked = true; |
3211 | | } |
3212 | | #endif |
3213 | | #ifdef COIN_HAS_CLP |
3214 | | // Possible change of pivot method |
3215 | | if (!savePivotMethod && !parentModel_) { |
3216 | | OsiClpSolverInterface * clpSolver |
3217 | | = dynamic_cast<OsiClpSolverInterface *> (solver_); |
3218 | | if (clpSolver && numberNodes_ >= 100 && numberNodes_ < 200) { |
3219 | | if (numberIterations_ < (numberSolves_ + numberNodes_)*10) { |
3220 | | //if (numberIterations_<numberNodes_*20) { |
3221 | | ClpSimplex * simplex = clpSolver->getModelPtr(); |
3222 | | ClpDualRowPivot * pivotMethod = simplex->dualRowPivot(); |
3223 | | ClpDualRowDantzig * pivot = |
3224 | | dynamic_cast< ClpDualRowDantzig*>(pivotMethod); |
3225 | | if (!pivot) { |
3226 | | savePivotMethod = pivotMethod->clone(true); |
3227 | | ClpDualRowDantzig dantzig; |
3228 | | simplex->setDualRowPivotAlgorithm(dantzig); |
3229 | | #ifdef COIN_DEVELOP |
3230 | | printf("%d node, %d iterations ->Dantzig\n", numberNodes_, |
3231 | | numberIterations_); |
3232 | | #endif |
3233 | | #ifdef CBC_THREAD |
3234 | | for (int i = 0; i < numberThreads_; i++) { |
3235 | | threadInfo[i].dantzigState = -1; |
3236 | | } |
3237 | | #endif |
3238 | | } |
3239 | | } |
3240 | | } |
3241 | | } |
3242 | | #endif |
3243 | | if (tree_->empty()) { |
3244 | | #ifdef CBC_THREAD |
3245 | | if (parallelMode() > 0) { |
3246 | | #ifdef COIN_DEVELOP |
3247 | | printf("empty\n"); |
3248 | | #endif |
3249 | | // may still be outstanding nodes |
3250 | | int iThread; |
3251 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3252 | | if (threadId[iThread].status) { |
3253 | | if (threadInfo[iThread].returnCode == 0) |
3254 | | break; |
3255 | | } |
3256 | | } |
3257 | | if (iThread < numberThreads_) { |
3258 | | #ifdef COIN_DEVELOP |
3259 | | printf("waiting for thread %d code 0\n", iThread); |
3260 | | #endif |
3261 | | if (parallelMode() > 0) { |
3262 | | unlockThread(); |
3263 | | locked = false; |
3264 | | } |
3265 | | pthread_cond_signal(threadInfo[iThread].condition2); // unlock in case |
3266 | | while (true) { |
3267 | | pthread_mutex_lock(&condition_mutex); |
3268 | | struct timespec absTime; |
3269 | | my_gettime(&absTime); |
3270 | | double time = absTime.tv_sec + 1.0e-9 * absTime.tv_nsec; |
3271 | | absTime.tv_nsec += 1000000; // millisecond |
3272 | | if (absTime.tv_nsec >= 1000000000) { |
3273 | | absTime.tv_nsec -= 1000000000; |
3274 | | absTime.tv_sec++; |
3275 | | } |
3276 | | pthread_cond_timedwait(&condition_main, &condition_mutex, &absTime); |
3277 | | my_gettime(&absTime); |
3278 | | double time2 = absTime.tv_sec + 1.0e-9 * absTime.tv_nsec; |
3279 | | timeWaiting += time2 - time; |
3280 | | pthread_mutex_unlock(&condition_mutex); |
3281 | | if (threadInfo[iThread].returnCode != 0) |
3282 | | break; |
3283 | | pthread_cond_signal(threadInfo[iThread].condition2); // unlock |
3284 | | } |
3285 | | threadModel[iThread]->moveToModel(this, 1); |
3286 | | assert (threadInfo[iThread].returnCode == 1); |
3287 | | if (threadInfo[iThread].dantzigState == -1) { |
3288 | | // 0 unset, -1 waiting to be set, 1 set |
3289 | | threadInfo[iThread].dantzigState = 1; |
3290 | | CbcModel * model = threadInfo[iThread].thisModel; |
3291 | | OsiClpSolverInterface * clpSolver2 |
3292 | | = dynamic_cast<OsiClpSolverInterface *> (model->solver()); |
3293 | | assert (clpSolver2); |
3294 | | ClpSimplex * simplex2 = clpSolver2->getModelPtr(); |
3295 | | ClpDualRowDantzig dantzig; |
3296 | | simplex2->setDualRowPivotAlgorithm(dantzig); |
3297 | | } |
3298 | | // say available |
3299 | | threadInfo[iThread].returnCode = -1; |
3300 | | threadStats[4]++; |
3301 | | #ifdef COIN_DEVELOP |
3302 | | printf("thread %d code now -1\n", iThread); |
3303 | | #endif |
3304 | | continue; |
3305 | | } else { |
3306 | | #ifdef COIN_DEVELOP |
3307 | | printf("no threads at code 0 \n"); |
3308 | | #endif |
3309 | | // now check if any have just finished |
3310 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3311 | | if (threadId[iThread].status) { |
3312 | | if (threadInfo[iThread].returnCode == 1) |
3313 | | break; |
3314 | | } |
3315 | | } |
3316 | | if (iThread < numberThreads_) { |
3317 | | if (parallelMode() > 0) { |
3318 | | unlockThread(); |
3319 | | locked = false; |
3320 | | } |
3321 | | threadModel[iThread]->moveToModel(this, 1); |
3322 | | assert (threadInfo[iThread].returnCode == 1); |
3323 | | // say available |
3324 | | threadInfo[iThread].returnCode = -1; |
3325 | | threadStats[4]++; |
3326 | | #ifdef COIN_DEVELOP |
3327 | | printf("thread %d code now -1\n", iThread); |
3328 | | #endif |
3329 | | continue; |
3330 | | } |
3331 | | } |
3332 | | if (!tree_->empty()) { |
3333 | | #ifdef COIN_DEVELOP |
3334 | | printf("tree not empty!!!!!!\n"); |
3335 | | #endif |
3336 | | continue; |
3337 | | } |
3338 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3339 | | if (threadId[iThread].status) { |
3340 | | if (threadInfo[iThread].returnCode != -1) { |
3341 | | printf("bad end of tree\n"); |
3342 | | abort(); |
3343 | | } |
3344 | | } |
3345 | | } |
3346 | | #ifdef COIN_DEVELOP |
3347 | | printf("finished ************\n"); |
3348 | | #endif |
3349 | | unlockThread(); |
3350 | | locked = false; // not needed as break |
3351 | | } |
3352 | | #endif |
3353 | | break; |
3354 | | } |
3355 | | #ifdef CBC_THREAD |
3356 | | if (parallelMode() > 0) { |
3357 | | unlockThread(); |
3358 | | locked = false; |
3359 | | } |
3360 | | #endif |
3361 | | // If done 100 nodes see if worth trying reduction |
3362 | | if (numberNodes_ == 50 || numberNodes_ == 100) { |
3363 | | #ifdef COIN_HAS_CLP |
3364 | | OsiClpSolverInterface * clpSolver |
3365 | | = dynamic_cast<OsiClpSolverInterface *> (solver_); |
3366 | | if (clpSolver && ((specialOptions_&131072) == 0) && true) { |
3367 | | ClpSimplex * simplex = clpSolver->getModelPtr(); |
3368 | | int perturbation = simplex->perturbation(); |
3369 | | #ifdef CLP_INVESTIGATE |
3370 | | printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n", |
3371 | | numberIterations_, saveNumberIterations, |
3372 | | numberSolves_, saveNumberSolves, perturbation); |
3373 | | #endif |
3374 | | if (perturbation == 50 && (numberIterations_ - saveNumberIterations) < |
3375 | | 8*(numberSolves_ - saveNumberSolves)) { |
3376 | | // switch off perturbation |
3377 | | simplex->setPerturbation(100); |
3378 | | #ifdef PERTURB_IN_FATHOM |
3379 | | // but allow in fathom |
3380 | | specialOptions_ |= 131072; |
3381 | | #endif |
3382 | | #ifdef CLP_INVESTIGATE |
3383 | | printf("Perturbation switched off\n"); |
3384 | | #endif |
3385 | | } |
3386 | | } |
3387 | | #endif |
3388 | | if (saveSolver) { |
3389 | | bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() && |
3390 | | (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart); |
3391 | | int numberColumns = getNumCols(); |
3392 | | if (tryNewSearch) { |
3393 | | checkCutoffForRestart = getCutoff() ; |
3394 | | #ifdef CLP_INVESTIGATE |
3395 | | printf("after %d nodes, cutoff %g - looking\n", |
3396 | | numberNodes_, getCutoff()); |
3397 | | #endif |
3398 | | saveSolver->resolve(); |
3399 | | double direction = saveSolver->getObjSense() ; |
3400 | | double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ; |
3401 | | double tolerance; |
3402 | | saveSolver->getDblParam(OsiDualTolerance, tolerance) ; |
3403 | | if (gap <= 0.0) |
3404 | | gap = tolerance; |
3405 | | gap += 100.0 * tolerance; |
3406 | | double integerTolerance = getDblParam(CbcIntegerTolerance) ; |
3407 | | |
3408 | | const double *lower = saveSolver->getColLower() ; |
3409 | | const double *upper = saveSolver->getColUpper() ; |
3410 | | const double *solution = saveSolver->getColSolution() ; |
3411 | | const double *reducedCost = saveSolver->getReducedCost() ; |
3412 | | |
3413 | | int numberFixed = 0 ; |
3414 | | int numberFixed2 = 0; |
3415 | | #ifdef COIN_DEVELOP |
3416 | | printf("gap %g\n", gap); |
3417 | | #endif |
3418 | | for (int i = 0 ; i < numberIntegers_ ; i++) { |
3419 | | int iColumn = integerVariable_[i] ; |
3420 | | double djValue = direction * reducedCost[iColumn] ; |
3421 | | if (upper[iColumn] - lower[iColumn] > integerTolerance) { |
3422 | | if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) { |
3423 | | //printf("%d to lb on dj of %g - bounds %g %g\n", |
3424 | | // iColumn,djValue,lower[iColumn],upper[iColumn]); |
3425 | | saveSolver->setColUpper(iColumn, lower[iColumn]) ; |
3426 | | numberFixed++ ; |
3427 | | } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) { |
3428 | | //printf("%d to ub on dj of %g - bounds %g %g\n", |
3429 | | // iColumn,djValue,lower[iColumn],upper[iColumn]); |
3430 | | saveSolver->setColLower(iColumn, upper[iColumn]) ; |
3431 | | numberFixed++ ; |
3432 | | } |
3433 | | } else { |
3434 | | //printf("%d has dj of %g - already fixed to %g\n", |
3435 | | // iColumn,djValue,lower[iColumn]); |
3436 | | numberFixed2++; |
3437 | | } |
3438 | | } |
3439 | | #ifdef COIN_DEVELOP |
3440 | | if ((specialOptions_&1) != 0) { |
3441 | | const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ; |
3442 | | if (debugger) { |
3443 | | printf("Contains optimal\n") ; |
3444 | | OsiSolverInterface * temp = saveSolver->clone(); |
3445 | | const double * solution = debugger->optimalSolution(); |
3446 | | const double *lower = temp->getColLower() ; |
3447 | | const double *upper = temp->getColUpper() ; |
3448 | | int n = temp->getNumCols(); |
3449 | | for (int i = 0; i < n; i++) { |
3450 | | if (temp->isInteger(i)) { |
3451 | | double value = floor(solution[i] + 0.5); |
3452 | | assert (value >= lower[i] && value <= upper[i]); |
3453 | | temp->setColLower(i, value); |
3454 | | temp->setColUpper(i, value); |
3455 | | } |
3456 | | } |
3457 | | temp->writeMps("reduced_fix"); |
3458 | | delete temp; |
3459 | | saveSolver->writeMps("reduced"); |
3460 | | } else { |
3461 | | abort(); |
3462 | | } |
3463 | | } |
3464 | | printf("Restart could fix %d integers (%d already fixed)\n", |
3465 | | numberFixed + numberFixed2, numberFixed2); |
3466 | | #endif |
3467 | | numberFixed += numberFixed2; |
3468 | | if (numberFixed*10 < numberColumns) |
3469 | | tryNewSearch = false; |
3470 | | } |
3471 | | if (tryNewSearch) { |
3472 | | // back to solver without cuts? |
3473 | | OsiSolverInterface * solver2 = saveSolver->clone(); |
3474 | | const double *lower = saveSolver->getColLower() ; |
3475 | | const double *upper = saveSolver->getColUpper() ; |
3476 | | for (int i = 0 ; i < numberIntegers_ ; i++) { |
3477 | | int iColumn = integerVariable_[i] ; |
3478 | | solver2->setColLower(iColumn, lower[iColumn]); |
3479 | | solver2->setColUpper(iColumn, upper[iColumn]); |
3480 | | } |
3481 | | // swap |
3482 | | delete saveSolver; |
3483 | | saveSolver = solver2; |
3484 | | double * newSolution = new double[numberColumns]; |
3485 | | double objectiveValue = checkCutoffForRestart; |
3486 | | CbcSerendipity heuristic(*this); |
3487 | | if (bestSolution_) |
3488 | | heuristic.setInputSolution(bestSolution_, bestObjective_); |
3489 | | heuristic.setFractionSmall(0.6); |
3490 | | heuristic.setFeasibilityPumpOptions(1008013); |
3491 | | // Use numberNodes to say how many are original rows |
3492 | | heuristic.setNumberNodes(continuousSolver_->getNumRows()); |
3493 | | #ifdef COIN_DEVELOP |
3494 | | if (continuousSolver_->getNumRows() < |
3495 | | solver_->getNumRows()) |
3496 | | printf("%d rows added ZZZZZ\n", |
3497 | | solver_->getNumRows() - continuousSolver_->getNumRows()); |
3498 | | #endif |
3499 | | int returnCode = heuristic.smallBranchAndBound(saveSolver, |
3500 | | -1, newSolution, |
3501 | | objectiveValue, |
3502 | | checkCutoffForRestart, "Reduce"); |
3503 | | if (returnCode < 0) { |
3504 | | #ifdef COIN_DEVELOP |
3505 | | printf("Restart - not small enough to do search after fixing\n"); |
3506 | | #endif |
3507 | | delete [] newSolution; |
3508 | | } else { |
3509 | | if ((returnCode&1) != 0) { |
3510 | | // increment number of solutions so other heuristics can test |
3511 | | numberSolutions_++; |
3512 | | numberHeuristicSolutions_++; |
3513 | | lastHeuristic_ = NULL; |
3514 | | setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ; |
3515 | | } |
3516 | | delete [] newSolution; |
3517 | | if (tree_->size()) { |
3518 | | double dummyBest; |
3519 | | tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ; |
3520 | | } |
3521 | | break; |
3522 | | } |
3523 | | } |
3524 | | delete saveSolver; |
3525 | | saveSolver = NULL; |
3526 | | } |
3527 | | } |
3528 | | /* |
3529 | | Check for abort on limits: node count, solution count, time, integrality gap. |
3530 | | */ |
3531 | | if (!(numberNodes_ < intParam_[CbcMaxNumNode] && |
3532 | | numberSolutions_ < intParam_[CbcMaxNumSol] && |
3533 | | !maximumSecondsReached() && |
3534 | | !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 || |
3535 | | numberIterations_ < maximumNumberIterations_))) { |
3536 | | // out of loop |
3537 | | break; |
3538 | | } |
3539 | | #ifdef BONMIN |
3540 | | assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible()); |
3541 | | #endif |
3542 | | if (cutoff > getCutoff()) { |
3543 | | double newCutoff = getCutoff(); |
3544 | | if (analyzeResults_) { |
3545 | | // see if we could fix any (more) |
3546 | | int n = 0; |
3547 | | double * newLower = analyzeResults_; |
3548 | | double * objLower = newLower + numberIntegers_; |
3549 | | double * newUpper = objLower + numberIntegers_; |
3550 | | double * objUpper = newUpper + numberIntegers_; |
3551 | | for (int i = 0; i < numberIntegers_; i++) { |
3552 | | if (objLower[i] > newCutoff) { |
3553 | | n++; |
3554 | | if (objUpper[i] > newCutoff) { |
3555 | | newCutoff = -COIN_DBL_MAX; |
3556 | | break; |
3557 | | } |
3558 | | } else if (objUpper[i] > newCutoff) { |
3559 | | n++; |
3560 | | } |
3561 | | } |
3562 | | if (newCutoff == -COIN_DBL_MAX) { |
3563 | | printf("Root analysis says finished\n"); |
3564 | | } else if (n > numberFixedNow_) { |
3565 | | printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n); |
3566 | | numberFixedNow_ = n; |
3567 | | } |
3568 | | } |
3569 | | if (eventHandler) { |
3570 | | if (!eventHandler->event(CbcEventHandler::solution)) { |
3571 | | eventHappened_ = true; // exit |
3572 | | } |
3573 | | newCutoff = getCutoff(); |
3574 | | } |
3575 | | if (parallelMode() > 0) |
3576 | | lockThread(); |
3577 | | // Do from deepest |
3578 | | tree_->cleanTree(this, newCutoff, bestPossibleObjective_) ; |
3579 | | nodeCompare_->newSolution(this) ; |
3580 | | nodeCompare_->newSolution(this, continuousObjective_, |
3581 | | continuousInfeasibilities_) ; |
3582 | | tree_->setComparison(*nodeCompare_) ; |
3583 | | if (tree_->empty()) { |
3584 | | if (parallelMode() > 0) |
3585 | | unlockThread(); |
3586 | | // For threads we need to check further |
3587 | | //break; // finished |
3588 | | continue; |
3589 | | } |
3590 | | if (parallelMode() > 0) |
3591 | | unlockThread(); |
3592 | | } |
3593 | | cutoff = getCutoff() ; |
3594 | | /* |
3595 | | Periodic activities: Opportunities to |
3596 | | + tweak the nodeCompare criteria, |
3597 | | + check if we've closed the integrality gap enough to quit, |
3598 | | + print a summary line to let the user know we're working |
3599 | | */ |
3600 | | if (numberNodes_ >= lastEvery1000) { |
3601 | | if (parallelMode() > 0) |
3602 | | lockThread(); |
3603 | | #ifdef COIN_HAS_CLP |
3604 | | // Possible change of pivot method |
3605 | | if (!savePivotMethod && !parentModel_) { |
3606 | | OsiClpSolverInterface * clpSolver |
3607 | | = dynamic_cast<OsiClpSolverInterface *> (solver_); |
3608 | | if (clpSolver && numberNodes_ >= 1000 && numberNodes_ < 2000) { |
3609 | | if (numberIterations_ < (numberSolves_ + numberNodes_)*10) { |
3610 | | ClpSimplex * simplex = clpSolver->getModelPtr(); |
3611 | | ClpDualRowPivot * pivotMethod = simplex->dualRowPivot(); |
3612 | | ClpDualRowDantzig * pivot = |
3613 | | dynamic_cast< ClpDualRowDantzig*>(pivotMethod); |
3614 | | if (!pivot) { |
3615 | | savePivotMethod = pivotMethod->clone(true); |
3616 | | ClpDualRowDantzig dantzig; |
3617 | | simplex->setDualRowPivotAlgorithm(dantzig); |
3618 | | #ifdef COIN_DEVELOP |
3619 | | printf("%d node, %d iterations ->Dantzig\n", numberNodes_, |
3620 | | numberIterations_); |
3621 | | #endif |
3622 | | #ifdef CBC_THREAD |
3623 | | for (int i = 0; i < numberThreads_; i++) { |
3624 | | threadInfo[i].dantzigState = -1; |
3625 | | } |
3626 | | #endif |
3627 | | } |
3628 | | } |
3629 | | } |
3630 | | } |
3631 | | #endif |
3632 | | lastEvery1000 = numberNodes_ + 1000; |
3633 | | bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ; |
3634 | | #ifdef CHECK_CUT_SIZE |
3635 | | verifyCutSize (tree_, *this); |
3636 | | #endif |
3637 | | // redo tree if wanted |
3638 | | if (redoTree) |
3639 | | tree_->setComparison(*nodeCompare_) ; |
3640 | | if (parallelMode() > 0) |
3641 | | unlockThread(); |
3642 | | } |
3643 | | if (saveCompare && !hotstartSolution_) { |
3644 | | // hotstart switched off |
3645 | | delete nodeCompare_; // off depth first |
3646 | | nodeCompare_ = saveCompare; |
3647 | | saveCompare = NULL; |
3648 | | // redo tree |
3649 | | if (parallelMode() > 0) |
3650 | | lockThread(); |
3651 | | tree_->setComparison(*nodeCompare_) ; |
3652 | | if (parallelMode() > 0) |
3653 | | unlockThread(); |
3654 | | } |
3655 | | if (numberNodes_ >= lastPrintEvery) { |
3656 | | lastPrintEvery = numberNodes_ + printFrequency_; |
3657 | | if (parallelMode() > 0) |
3658 | | lockThread(); |
3659 | | int nNodes = tree_->size() ; |
3660 | | |
3661 | | //MODIF PIERRE |
3662 | | bestPossibleObjective_ = tree_->getBestPossibleObjective(); |
3663 | | if (parallelMode() > 0) |
3664 | | unlockThread(); |
3665 | | #ifdef CLP_INVESTIGATE |
3666 | | if (getCutoff() < 1.0e20) { |
3667 | | if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 && |
3668 | | !parentModel_) |
3669 | | printf("model cutoff in status %g, best %g, increment %g\n", |
3670 | | getCutoff(), bestObjective_, getCutoffIncrement()); |
3671 | | assert (getCutoff() < bestObjective_ - getCutoffIncrement() + |
3672 | | 1.0e-6 + 1.0e-10*fabs(bestObjective_)); |
3673 | | } |
3674 | | #endif |
3675 | | if (!intParam_[CbcPrinting]) { |
3676 | | messageHandler()->message(CBC_STATUS, messages()) |
3677 | | << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_ |
3678 | | << getCurrentSeconds() |
3679 | | << CoinMessageEol ; |
3680 | | } else if (intParam_[CbcPrinting] == 1) { |
3681 | | messageHandler()->message(CBC_STATUS2, messages()) |
3682 | | << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_ |
3683 | | << lastDepth << lastUnsatisfied << numberIterations_ |
3684 | | << getCurrentSeconds() |
3685 | | << CoinMessageEol ; |
3686 | | } else if (!numberExtraIterations_) { |
3687 | | messageHandler()->message(CBC_STATUS2, messages()) |
3688 | | << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_ |
3689 | | << lastDepth << lastUnsatisfied << numberIterations_ |
3690 | | << getCurrentSeconds() |
3691 | | << CoinMessageEol ; |
3692 | | } else { |
3693 | | messageHandler()->message(CBC_STATUS3, messages()) |
3694 | | << numberNodes_ << numberExtraNodes_ << nNodes << bestObjective_ << bestPossibleObjective_ |
3695 | | << lastDepth << lastUnsatisfied << numberIterations_ << numberExtraIterations_ |
3696 | | << getCurrentSeconds() |
3697 | | << CoinMessageEol ; |
3698 | | } |
3699 | | if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) { |
3700 | | eventHappened_ = true; // exit |
3701 | | } |
3702 | | } |
3703 | | // See if can stop on gap |
3704 | | double testGap = CoinMax(dblParam_[CbcAllowableGap], |
3705 | | CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_)) |
3706 | | * dblParam_[CbcAllowableFractionGap]); |
3707 | | if (bestObjective_ - bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) { |
3708 | | stoppedOnGap_ = true ; |
3709 | | } |
3710 | | |
3711 | | #ifdef CHECK_NODE_FULL |
3712 | | verifyTreeNodes(tree_, *this) ; |
3713 | | # endif |
3714 | | # ifdef CHECK_CUT_COUNTS |
3715 | | verifyCutCounts(tree_, *this) ; |
3716 | | # endif |
3717 | | /* |
3718 | | Now we come to the meat of the loop. To create the active subproblem, we'll |
3719 | | pop the most promising node in the live set, rebuild the subproblem it |
3720 | | represents, and then execute the current arm of the branch to create the |
3721 | | active subproblem. |
3722 | | */ |
3723 | | CbcNode * node = NULL; |
3724 | | #ifdef CBC_THREAD |
3725 | | if (!parallelMode() || parallelMode() == -1) { |
3726 | | #endif |
3727 | | node = tree_->bestNode(cutoff) ; |
3728 | | // Possible one on tree worse than cutoff |
3729 | | if (!node || node->objectiveValue() > cutoff) |
3730 | | continue; |
3731 | | // Do main work of solving node here |
3732 | | doOneNode(this, node, createdNode); |
3733 | | #ifdef CBC_THREAD |
3734 | | } else if (parallelMode() > 0) { |
3735 | | node = tree_->bestNode(cutoff) ; |
3736 | | // Possible one on tree worse than cutoff |
3737 | | if (!node || node->objectiveValue() > cutoff) |
3738 | | continue; |
3739 | | threadStats[0]++; |
3740 | | //need to think |
3741 | | int iThread; |
3742 | | // Start one off if any available |
3743 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3744 | | if (threadInfo[iThread].returnCode == -1) |
3745 | | break; |
3746 | | } |
3747 | | if (iThread < numberThreads_) { |
3748 | | threadInfo[iThread].node = node; |
3749 | | assert (threadInfo[iThread].returnCode == -1); |
3750 | | // say in use |
3751 | | threadModel[iThread]->moveToModel(this, 0); |
3752 | | // This has to be AFTER moveToModel |
3753 | | threadInfo[iThread].returnCode = 0; |
3754 | | pthread_cond_signal(threadInfo[iThread].condition2); // unlock |
3755 | | threadCount[iThread]++; |
3756 | | } |
3757 | | lockThread(); |
3758 | | locked = true; |
3759 | | // see if any finished |
3760 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3761 | | if (threadInfo[iThread].returnCode > 0) |
3762 | | break; |
3763 | | } |
3764 | | unlockThread(); |
3765 | | locked = false; |
3766 | | if (iThread < numberThreads_) { |
3767 | | threadModel[iThread]->moveToModel(this, 1); |
3768 | | assert (threadInfo[iThread].returnCode == 1); |
3769 | | // say available |
3770 | | threadInfo[iThread].returnCode = -1; |
3771 | | // carry on |
3772 | | threadStats[3]++; |
3773 | | } else { |
3774 | | // Start one off if any available |
3775 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3776 | | if (threadInfo[iThread].returnCode == -1) |
3777 | | break; |
3778 | | } |
3779 | | if (iThread < numberThreads_) { |
3780 | | lockThread(); |
3781 | | locked = true; |
3782 | | // If any on tree get |
3783 | | if (!tree_->empty()) { |
3784 | | //node = tree_->bestNode(cutoff) ; |
3785 | | //assert (node); |
3786 | | threadStats[1]++; |
3787 | | continue; // ** get another node |
3788 | | } |
3789 | | unlockThread(); |
3790 | | locked = false; |
3791 | | } |
3792 | | // wait (for debug could sleep and use test) |
3793 | | bool finished = false; |
3794 | | while (!finished) { |
3795 | | pthread_mutex_lock(&condition_mutex); |
3796 | | struct timespec absTime; |
3797 | | my_gettime(&absTime); |
3798 | | double time = absTime.tv_sec + 1.0e-9 * absTime.tv_nsec; |
3799 | | absTime.tv_nsec += 1000000; // millisecond |
3800 | | if (absTime.tv_nsec >= 1000000000) { |
3801 | | absTime.tv_nsec -= 1000000000; |
3802 | | absTime.tv_sec++; |
3803 | | } |
3804 | | pthread_cond_timedwait(&condition_main, &condition_mutex, &absTime); |
3805 | | my_gettime(&absTime); |
3806 | | double time2 = absTime.tv_sec + 1.0e-9 * absTime.tv_nsec; |
3807 | | timeWaiting += time2 - time; |
3808 | | pthread_mutex_unlock(&condition_mutex); |
3809 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3810 | | if (threadInfo[iThread].returnCode > 0) { |
3811 | | finished = true; |
3812 | | break; |
3813 | | } else if (threadInfo[iThread].returnCode == 0) { |
3814 | | pthread_cond_signal(threadInfo[iThread].condition2); // unlock |
3815 | | } |
3816 | | } |
3817 | | } |
3818 | | assert (iThread < numberThreads_); |
3819 | | // move information to model |
3820 | | threadModel[iThread]->moveToModel(this, 1); |
3821 | | node = threadInfo[iThread].node; |
3822 | | threadInfo[iThread].node = NULL; |
3823 | | assert (threadInfo[iThread].returnCode == 1); |
3824 | | // say available |
3825 | | threadInfo[iThread].returnCode = -1; |
3826 | | // carry on |
3827 | | threadStats[2]++; |
3828 | | } |
3829 | | } else { |
3830 | | // Deterministic parallel |
3831 | | if (tree_->size() < CoinMin(numberThreads_, 8) && !goneParallel) { |
3832 | | node = tree_->bestNode(cutoff) ; |
3833 | | // Possible one on tree worse than cutoff |
3834 | | if (!node || node->objectiveValue() > cutoff) |
3835 | | continue; |
3836 | | // Do main work of solving node here |
3837 | | doOneNode(this, node, createdNode); |
3838 | | assert (createdNode); |
3839 | | if (!createdNode->active()) { |
3840 | | delete createdNode; |
3841 | | createdNode = NULL; |
3842 | | } else { |
3843 | | // Say one more pointing to this |
3844 | | node->nodeInfo()->increment() ; |
3845 | | tree_->push(createdNode) ; |
3846 | | } |
3847 | | if (node->active()) { |
3848 | | assert (node->nodeInfo()); |
3849 | | if (node->nodeInfo()->numberBranchesLeft()) { |
3850 | | tree_->push(node) ; |
3851 | | } else { |
3852 | | node->setActive(false); |
3853 | | } |
3854 | | } else { |
3855 | | if (node->nodeInfo()) { |
3856 | | if (!node->nodeInfo()->numberBranchesLeft()) |
3857 | | node->nodeInfo()->allBranchesGone(); // can clean up |
3858 | | // So will delete underlying stuff |
3859 | | node->setActive(true); |
3860 | | } |
3861 | | delNode[nDeleteNode++] = node; |
3862 | | node = NULL; |
3863 | | } |
3864 | | if (nDeleteNode >= MAX_DEL_NODE) { |
3865 | | for (int i = 0; i < nDeleteNode; i++) { |
3866 | | //printf("trying to del %d %x\n",i,delNode[i]); |
3867 | | delete delNode[i]; |
3868 | | //printf("done to del %d %x\n",i,delNode[i]); |
3869 | | } |
3870 | | nDeleteNode = 0; |
3871 | | } |
3872 | | } else { |
3873 | | // Split |
3874 | | int saveTreeSize = tree_->size(); |
3875 | | goneParallel = true; |
3876 | | int nAffected = splitModel(numberThreads_, threadModel, defaultParallelNodes); |
3877 | | int iThread; |
3878 | | // do all until finished |
3879 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3880 | | // obviously tune |
3881 | | threadInfo[iThread].nDeleteNode = defaultParallelIterations; |
3882 | | } |
3883 | | // Save current state |
3884 | | int iObject; |
3885 | | for (iObject = 0; iObject < numberObjects_; iObject++) { |
3886 | | saveObjects[iObject]->updateBefore(object_[iObject]); |
3887 | | } |
3888 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3889 | | threadInfo[iThread].returnCode = 0; |
3890 | | pthread_cond_signal(threadInfo[iThread].condition2); // unlock |
3891 | | } |
3892 | | // wait |
3893 | | bool finished = false; |
3894 | | while (!finished) { |
3895 | | pthread_mutex_lock(&condition_mutex); |
3896 | | struct timespec absTime; |
3897 | | my_gettime(&absTime); |
3898 | | double time = absTime.tv_sec + 1.0e-9 * absTime.tv_nsec; |
3899 | | absTime.tv_nsec += 1000000; // millisecond |
3900 | | if (absTime.tv_nsec >= 1000000000) { |
3901 | | absTime.tv_nsec -= 1000000000; |
3902 | | absTime.tv_sec++; |
3903 | | } |
3904 | | pthread_cond_timedwait(&condition_main, &condition_mutex, &absTime); |
3905 | | my_gettime(&absTime); |
3906 | | double time2 = absTime.tv_sec + 1.0e-9 * absTime.tv_nsec; |
3907 | | timeWaiting += time2 - time; |
3908 | | pthread_mutex_unlock(&condition_mutex); |
3909 | | finished = true; |
3910 | | for (iThread = 0; iThread < numberThreads_; iThread++) { |
3911 | | if (threadInfo[iThread].returnCode <= 0) { |
3912 | | finished = false; |
3913 | | } |
3914 | | } |
3915 | | } |
3916 | | // Unmark marked |
3917 | | for (int i = 0; i < nAffected; i++) { |
3918 | | walkback_[i]->unmark(); |
3919 | | } |
3920 | | int iModel; |
3921 | | double scaleFactor = 1.0; |
3922 | | for (iModel = 0; iModel < numberThreads_; iModel++) { |
3923 | | //printf("model %d tree size %d\n",iModel,threadModel[iModel]->tree_->size()); |
3924 | | if (saveTreeSize > 4*numberThreads_*defaultParallelNodes) { |
3925 | | if (!threadModel[iModel]->tree_->size()) { |
3926 | | scaleFactor *= 1.05; |
3927 | | } |
3928 | | } |
3929 | | threadModel[iModel]->moveToModel(this, 11); |
3930 | | // Update base model |
3931 | | OsiObject ** threadObject = threadModel[iModel]->object_; |
3932 | | for (iObject = 0; iObject < numberObjects_; iObject++) { |
3933 | | object_[iObject]->updateAfter(threadObject[iObject], saveObjects[iObject]); |
3934 | | } |
3935 | | } |
3936 | | if (scaleFactor != 1.0) { |
3937 | | int newNumber = static_cast<int> (defaultParallelNodes * scaleFactor + 0.5001); |
3938 | | if (newNumber*2 < defaultParallelIterations) { |
3939 | | if (defaultParallelNodes == 1) |
3940 | | newNumber = 2; |
3941 | | if (newNumber != defaultParallelNodes) { |
3942 | | char general[200]; |
3943 | | sprintf(general, "Changing tree size from %d to %d", |
3944 | | defaultParallelNodes, newNumber); |
3945 | | messageHandler()->message(CBC_GENERAL, |
3946 | | messages()) |
3947 | | << general << CoinMessageEol ; |
3948 | | defaultParallelNodes = newNumber; |
3949 | | } |
3950 | | } |
3951 | | } |
3952 | | //printf("Tree sizes %d %d %d - affected %d\n",saveTreeSize,saveTreeSize2,tree_->size(),nAffected); |
3953 | | } |
3954 | | } |
3955 | | #endif |
3956 | | } |
3957 | | if (nDeleteNode) { |
3958 | | for (int i = 0; i < nDeleteNode; i++) { |
3959 | | delete delNode[i]; |
3960 | | } |
3961 | | nDeleteNode = 0; |
3962 | | } |
3963 | | #ifdef CBC_THREAD |
3964 | | if (numberThreads_) { |
3965 | | int i; |
3966 | | // Seems to be bug in CoinCpu on Linux - does threads as well despite documentation |
3967 | | double time = 0.0; |
3968 | | for (i = 0; i < numberThreads_; i++) |
3969 | | time += threadInfo[i].timeInThread; |
3970 | | bool goodTimer = time < (getCurrentSeconds()); |
3971 | | for (i = 0; i < numberThreads_; i++) { |
3972 | | while (threadInfo[i].returnCode == 0) { |
3973 | | pthread_cond_signal(threadInfo[i].condition2); // unlock |
3974 | | pthread_mutex_lock(&condition_mutex); |
3975 | | struct timespec absTime; |
3976 | | my_gettime(&absTime); |
3977 | | absTime.tv_nsec += 1000000; // millisecond |
3978 | | if (absTime.tv_nsec >= 1000000000) { |
3979 | | absTime.tv_nsec -= 1000000000; |
3980 | | absTime.tv_sec++; |
3981 | | } |
3982 | | pthread_cond_timedwait(&condition_main, &condition_mutex, &absTime); |
3983 | | my_gettime(&absTime); |
3984 | | pthread_mutex_unlock(&condition_mutex); |
3985 | | } |
3986 | | pthread_cond_signal(threadInfo[i].condition2); // unlock |
3987 | | pthread_mutex_lock(&condition_mutex); // not sure necessary but have had one hang on interrupt |
3988 | | threadModel[i]->numberThreads_ = 0; // say exit |
3989 | | if (parallelMode() < 0) |
3990 | | delete [] threadInfo[i].delNode; |
3991 | | threadInfo[i].returnCode = 0; |
3992 | | pthread_mutex_unlock(&condition_mutex); |
3993 | | pthread_cond_signal(threadInfo[i].condition2); // unlock |
3994 | | //if (!stopped) |
3995 | | //pthread_join(threadId[i],NULL); |
3996 | | int returnCode; |
3997 | | returnCode = pthread_join(threadId[i].thr, NULL); |
3998 | | threadId[i].status = 0; |
3999 | | assert (!returnCode); |
4000 | | //else |
4001 | | //pthread_kill(threadId[i]); // kill rather than try and synchronize |
4002 | | threadModel[i]->moveToModel(this, 2); |
4003 | | pthread_mutex_destroy (threadInfo[i].mutex2); |
4004 | | pthread_cond_destroy (threadInfo[i].condition2); |
4005 | | assert (threadInfo[i].numberTimesLocked == threadInfo[i].numberTimesUnlocked); |
4006 | | handler_->message(CBC_THREAD_STATS, messages_) |
4007 | | << "Thread"; |
4008 | | handler_->printing(true) |
4009 | | << i << threadCount[i] << threadInfo[i].timeWaitingToStart; |
4010 | | handler_->printing(goodTimer) << threadInfo[i].timeInThread; |
4011 | | handler_->printing(false) << 0.0; |
4012 | | handler_->printing(true) << threadInfo[i].numberTimesLocked |
4013 | | << threadInfo[i].timeLocked << threadInfo[i].timeWaitingToLock |
4014 | | << CoinMessageEol; |
4015 | | } |
4016 | | assert (threadInfo[numberThreads_].numberTimesLocked == threadInfo[numberThreads_].numberTimesUnlocked); |
4017 | | handler_->message(CBC_THREAD_STATS, messages_) |
4018 | | << "Main thread"; |
4019 | | handler_->printing(false) << 0 << 0 << 0.0; |
4020 | | handler_->printing(false) << 0.0; |
4021 | | handler_->printing(true) << timeWaiting; |
4022 | | handler_->printing(true) << threadInfo[numberThreads_].numberTimesLocked |
4023 | | << threadInfo[numberThreads_].timeLocked << threadInfo[numberThreads_].timeWaitingToLock |
4024 | | << CoinMessageEol; |
4025 | | pthread_mutex_destroy (&mutex); |
4026 | | pthread_cond_destroy (&condition_main); |
4027 | | pthread_mutex_destroy (&condition_mutex); |
4028 | | // delete models (here in case some point to others) |
4029 | | for (i = 0; i < numberThreads_; i++) { |
4030 | | // make sure handler will be deleted |
4031 | | threadModel[i]->defaultHandler_ = true; |
4032 | | delete threadModel[i]; |
4033 | | } |
4034 | | delete [] mutex2; |
4035 | | delete [] condition2; |
4036 | | delete [] threadId; |
4037 | | delete [] threadInfo; |
4038 | | delete [] threadModel; |
4039 | | delete [] threadCount; |
4040 | | mutex_ = NULL; |
4041 | | // adjust time to allow for children on some systems |
4042 | | dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren(); |
4043 | | } |
4044 | | #endif |
4045 | | /* |
4046 | | End of the non-abort actions. The next block of code is executed if we've |
4047 | | aborted because we hit one of the limits. Clean up by deleting the live set |
4048 | | and break out of the node processing loop. Note that on an abort, node may |
4049 | | have been pushed back onto the tree for further processing, in which case |
4050 | | it'll be deleted in cleanTree. We need to check. |
4051 | | */ |
4052 | | if (!(numberNodes_ < intParam_[CbcMaxNumNode] && |
4053 | | numberSolutions_ < intParam_[CbcMaxNumSol] && |
4054 | | !maximumSecondsReached() && |
4055 | | !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 || |
4056 | | numberIterations_ < maximumNumberIterations_))) { |
4057 | | if (tree_->size()) { |
4058 | | double dummyBest; |
4059 | | tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ; |
4060 | | } |
4061 | | delete nextRowCut_; |
4062 | | if (stoppedOnGap_) { |
4063 | | messageHandler()->message(CBC_GAP, messages()) |
4064 | | << bestObjective_ - bestPossibleObjective_ |
4065 | | << dblParam_[CbcAllowableGap] |
4066 | | << dblParam_[CbcAllowableFractionGap]*100.0 |
4067 | | << CoinMessageEol ; |
4068 | | secondaryStatus_ = 2; |
4069 | | status_ = 0 ; |
4070 | | } else if (isNodeLimitReached()) { |
4071 | | handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ; |
4072 | | secondaryStatus_ = 3; |
4073 | | status_ = 1 ; |
4074 | | } else if (maximumSecondsReached()) { |
4075 | | handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ; |
4076 | | secondaryStatus_ = 4; |
4077 | | status_ = 1 ; |
4078 | | } else if (eventHappened_) { |
4079 | | handler_->message(CBC_EVENT, messages_) << CoinMessageEol ; |
4080 | | secondaryStatus_ = 5; |
4081 | | status_ = 5 ; |
4082 | | } else { |
4083 | | handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ; |
4084 | | secondaryStatus_ = 6; |
4085 | | status_ = 1 ; |
4086 | | } |
4087 | | } |
4088 | | /* |
4089 | | That's it, we've exhausted the search tree, or broken out of the loop because |
4090 | | we hit some limit on evaluation. |
4091 | | |
4092 | | We may have got an intelligent tree so give it one more chance |
4093 | | */ |
4094 | | // Tell solver we are not in Branch and Cut |
4095 | | solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ; |
4096 | | tree_->endSearch(); |
4097 | | // If we did any sub trees - did we give up on any? |
4098 | | if ( numberStoppedSubTrees_) |
4099 | | status_ = 1; |
4100 | | numberNodes_ += numberExtraNodes_; |
4101 | | numberIterations_ += numberExtraIterations_; |
4102 | | if (eventHandler) { |
4103 | | eventHandler->event(CbcEventHandler::endSearch); |
4104 | | } |
4105 | | if (!status_) { |
4106 | | // Set best possible unless stopped on gap |
4107 | | if (secondaryStatus_ != 2) |
4108 | | bestPossibleObjective_ = bestObjective_; |
4109 | | handler_->message(CBC_END_GOOD, messages_) |
4110 | | << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds() |
4111 | | << CoinMessageEol ; |
4112 | | } else { |
4113 | | handler_->message(CBC_END, messages_) |
4114 | | << bestObjective_ << bestPossibleObjective_ |
4115 | | << numberIterations_ << numberNodes_ << getCurrentSeconds() |
4116 | | << CoinMessageEol ; |
4117 | | } |
4118 | | if (numberStrongIterations_) |
4119 | | handler_->message(CBC_STRONG_STATS, messages_) |
4120 | | << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2] |
4121 | | << strongInfo_[1] << CoinMessageEol ; |
4122 | | if (!numberExtraNodes_) |
4123 | | handler_->message(CBC_OTHER_STATS, messages_) |
4124 | | << maximumDepthActual_ |
4125 | | << numberDJFixed_ << CoinMessageEol ; |
4126 | | else |
4127 | | handler_->message(CBC_OTHER_STATS2, messages_) |
4128 | | << maximumDepthActual_ |
4129 | | << numberDJFixed_ << numberExtraNodes_ << numberExtraIterations_ |
4130 | | << CoinMessageEol ; |
4131 | | if (doStatistics == 100) { |
4132 | | for (int i = 0; i < numberObjects_; i++) { |
4133 | | CbcSimpleIntegerDynamicPseudoCost * obj = |
4134 | | dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ; |
4135 | | if (obj) |
4136 | | obj->print(); |
4137 | | } |
4138 | | } |
4139 | | if (statistics_) { |
4140 | | // report in some way |
4141 | | int * lookup = new int[numberObjects_]; |
4142 | | int i; |
4143 | | for (i = 0; i < numberObjects_; i++) |
4144 | | lookup[i] = -1; |
4145 | | bool goodIds = false; //true; |
4146 | | for (i = 0; i < numberObjects_; i++) { |
4147 | | int iColumn = object_[i]->columnNumber(); |
4148 | | if (iColumn >= 0 && iColumn < numberColumns) { |
4149 | | if (lookup[i] == -1) { |
4150 | | lookup[i] = iColumn; |
4151 | | } else { |
4152 | | goodIds = false; |
4153 | | break; |
4154 | | } |
4155 | | } else { |
4156 | | goodIds = false; |
4157 | | break; |
4158 | | } |
4159 | | } |
4160 | | if (!goodIds) { |
4161 | | delete [] lookup; |
4162 | | lookup = NULL; |
4163 | | } |
4164 | | if (doStatistics >= 3) { |
4165 | | printf(" node parent depth column value obj inf\n"); |
4166 | | for ( i = 0; i < numberNodes2_; i++) { |
4167 | | statistics_[i]->print(lookup); |
4168 | | } |
4169 | | } |
4170 | | if (doStatistics > 1) { |
4171 | | // Find last solution |
4172 | | int k; |
4173 | | for (k = numberNodes2_ - 1; k >= 0; k--) { |
4174 | | if (statistics_[k]->endingObjective() != COIN_DBL_MAX && |
4175 | | !statistics_[k]->endingInfeasibility()) |
4176 | | break; |
4177 | | } |
4178 | | if (k >= 0) { |
4179 | | int depth = statistics_[k]->depth(); |
4180 | | int * which = new int[depth+1]; |
4181 | | for (i = depth; i >= 0; i--) { |
4182 | | which[i] = k; |
4183 | | k = statistics_[k]->parentNode(); |
4184 | | } |
4185 | | printf(" node parent depth column value obj inf\n"); |
4186 | | for (i = 0; i <= depth; i++) { |
4187 | | statistics_[which[i]]->print(lookup); |
4188 | | } |
4189 | | delete [] which; |
4190 | | } |
4191 | | } |
4192 | | // now summary |
4193 | | int maxDepth = 0; |
4194 | | double averageSolutionDepth = 0.0; |
4195 | | int numberSolutions = 0; |
4196 | | double averageCutoffDepth = 0.0; |
4197 | | double averageSolvedDepth = 0.0; |
4198 | | int numberCutoff = 0; |
4199 | | int numberDown = 0; |
4200 | | int numberFirstDown = 0; |
4201 | | double averageInfDown = 0.0; |
4202 | | double averageObjDown = 0.0; |
4203 | | int numberCutoffDown = 0; |
4204 | | int numberUp = 0; |
4205 | | int numberFirstUp = 0; |
4206 | | double averageInfUp = 0.0; |
4207 | | double averageObjUp = 0.0; |
4208 | | int numberCutoffUp = 0; |
4209 | | double averageNumberIterations1 = 0.0; |
4210 | | double averageValue = 0.0; |
4211 | | for ( i = 0; i < numberNodes2_; i++) { |
4212 | | int depth = statistics_[i]->depth(); |
4213 | | int way = statistics_[i]->way(); |
4214 | | double value = statistics_[i]->value(); |
4215 | | double startingObjective = statistics_[i]->startingObjective(); |
4216 | | int startingInfeasibility = statistics_[i]->startingInfeasibility(); |
4217 | | double endingObjective = statistics_[i]->endingObjective(); |
4218 | | int endingInfeasibility = statistics_[i]->endingInfeasibility(); |
4219 | | maxDepth = CoinMax(depth, maxDepth); |
4220 | | // Only for completed |
4221 | | averageNumberIterations1 += statistics_[i]->numberIterations(); |
4222 | | averageValue += value; |
4223 | | if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) { |
4224 | | numberSolutions++; |
4225 | | averageSolutionDepth += depth; |
4226 | | } |
4227 | | if (endingObjective == COIN_DBL_MAX) { |
4228 | | numberCutoff++; |
4229 | | averageCutoffDepth += depth; |
4230 | | if (way < 0) { |
4231 | | numberDown++; |
4232 | | numberCutoffDown++; |
4233 | | if (way == -1) |
4234 | | numberFirstDown++; |
4235 | | } else { |
4236 | | numberUp++; |
4237 | | numberCutoffUp++; |
4238 | | if (way == 1) |
4239 | | numberFirstUp++; |
4240 | | } |
4241 | | } else { |
4242 | | averageSolvedDepth += depth; |
4243 | | if (way < 0) { |
4244 | | numberDown++; |
4245 | | averageInfDown += startingInfeasibility - endingInfeasibility; |
4246 | | averageObjDown += endingObjective - startingObjective; |
4247 | | if (way == -1) |
4248 | | numberFirstDown++; |
4249 | | } else { |
4250 | | numberUp++; |
4251 | | averageInfUp += startingInfeasibility - endingInfeasibility; |
4252 | | averageObjUp += endingObjective - startingObjective; |
4253 | | if (way == 1) |
4254 | | numberFirstUp++; |
4255 | | } |
4256 | | } |
4257 | | } |
4258 | | // Now print |
4259 | | if (numberSolutions) |
4260 | | averageSolutionDepth /= static_cast<double> (numberSolutions); |
4261 | | int numberSolved = numberNodes2_ - numberCutoff; |
4262 | | double averageNumberIterations2 = numberIterations_ - averageNumberIterations1 |
4263 | | - numberIterationsAtContinuous; |
4264 | | if (numberCutoff) { |
4265 | | averageCutoffDepth /= static_cast<double> (numberCutoff); |
4266 | | averageNumberIterations2 /= static_cast<double> (numberCutoff); |
4267 | | } |
4268 | | if (numberNodes2_) |
4269 | | averageValue /= static_cast<double> (numberNodes2_); |
4270 | | if (numberSolved) { |
4271 | | averageNumberIterations1 /= static_cast<double> (numberSolved); |
4272 | | averageSolvedDepth /= static_cast<double> (numberSolved); |
4273 | | } |
4274 | | printf("%d solution(s) were found (by branching) at an average depth of %g\n", |
4275 | | numberSolutions, averageSolutionDepth); |
4276 | | printf("average value of variable being branched on was %g\n", |
4277 | | averageValue); |
4278 | | printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n", |
4279 | | numberCutoff, averageCutoffDepth, averageNumberIterations2); |
4280 | | printf("%d nodes were solved at an average depth of %g with iteration count of %g\n", |
4281 | | numberSolved, averageSolvedDepth, averageNumberIterations1); |
4282 | | if (numberDown) { |
4283 | | averageInfDown /= static_cast<double> (numberDown); |
4284 | | averageObjDown /= static_cast<double> (numberDown); |
4285 | | } |
4286 | | printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n", |
4287 | | numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown, |
4288 | | averageInfDown, averageObjDown); |
4289 | | if (numberUp) { |
4290 | | averageInfUp /= static_cast<double> (numberUp); |
4291 | | averageObjUp /= static_cast<double> (numberUp); |
4292 | | } |
4293 | | printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n", |
4294 | | numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp, |
4295 | | averageInfUp, averageObjUp); |
4296 | | for ( i = 0; i < numberNodes2_; i++) |
4297 | | delete statistics_[i]; |
4298 | | delete [] statistics_; |
4299 | | statistics_ = NULL; |
4300 | | maximumStatistics_ = 0; |
4301 | | delete [] lookup; |
4302 | | } |
4303 | | /* |
4304 | | If we think we have a solution, restore and confirm it with a call to |
4305 | | setBestSolution(). We need to reset the cutoff value so as not to fathom |
4306 | | the solution on bounds. Note that calling setBestSolution( ..., true) |
4307 | | leaves the continuousSolver_ bounds vectors fixed at the solution value. |
4308 | | |
4309 | | Running resolve() here is a failsafe --- setBestSolution has already |
4310 | | reoptimised using the continuousSolver_. If for some reason we fail to |
4311 | | prove optimality, run the problem again after instructing the solver to |
4312 | | tell us more. |
4313 | | |
4314 | | If all looks good, replace solver_ with continuousSolver_, so that the |
4315 | | outside world will be able to obtain information about the solution using |
4316 | | public methods. |
4317 | | */ |
4318 | | if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4)) { |
4319 | | setCutoff(1.0e50) ; // As best solution should be worse than cutoff |
4320 | | phase_ = 5; |
4321 | | double increment = getDblParam(CbcModel::CbcCutoffIncrement) ; |
4322 | | if ((specialOptions_&4) == 0) |
4323 | | bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve |
4324 | | setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ; |
4325 | | continuousSolver_->resolve() ; |
4326 | | if (!continuousSolver_->isProvenOptimal()) { |
4327 | | continuousSolver_->messageHandler()->setLogLevel(2) ; |
4328 | | continuousSolver_->initialSolve() ; |
4329 | | } |
4330 | | delete solver_ ; |
4331 | | // above deletes solverCharacteristics_ |
4332 | | solverCharacteristics_ = NULL; |
4333 | | solver_ = continuousSolver_ ; |
4334 | | setPointers(solver_); |
4335 | | continuousSolver_ = NULL ; |
4336 | | } |
4337 | | /* |
4338 | | Clean up dangling objects. continuousSolver_ may already be toast. |
4339 | | */ |
4340 | | delete lastws ; |
4341 | | if (saveObjects) { |
4342 | | for (int i = 0; i < numberObjects_; i++) |
4343 | | delete saveObjects[i]; |
4344 | | delete [] saveObjects; |
4345 | | } |
4346 | | numberStrong_ = saveNumberStrong; |
4347 | | numberBeforeTrust_ = saveNumberBeforeTrust; |
4348 | | delete [] whichGenerator_ ; |
4349 | | whichGenerator_ = NULL; |
4350 | | delete [] lowerBefore ; |
4351 | | delete [] upperBefore ; |
4352 | | delete [] walkback_ ; |
4353 | | walkback_ = NULL ; |
4354 | | delete [] lastNodeInfo_ ; |
4355 | | lastNodeInfo_ = NULL; |
4356 | | delete [] lastNumberCuts_ ; |
4357 | | lastNumberCuts_ = NULL; |
4358 | | delete [] lastCut_; |
4359 | | lastCut_ = NULL; |
4360 | | delete [] addedCuts_ ; |
4361 | | addedCuts_ = NULL ; |
4362 | | //delete persistentInfo; |
4363 | | // Get rid of characteristics |
4364 | | solverCharacteristics_ = NULL; |
4365 | | if (continuousSolver_) { |
4366 | | delete continuousSolver_ ; |
4367 | | continuousSolver_ = NULL ; |
4368 | | } |
4369 | | /* |
4370 | | Destroy global cuts by replacing with an empty OsiCuts object. |
4371 | | */ |
4372 | | globalCuts_ = OsiCuts() ; |
4373 | | if (!bestSolution_) { |
4374 | | // make sure lp solver is infeasible |
4375 | | int numberColumns = solver_->getNumCols(); |
4376 | | const double * columnLower = solver_->getColLower(); |
4377 | | int iColumn; |
4378 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
4379 | | if (solver_->isInteger(iColumn)) |
4380 | | solver_->setColUpper(iColumn, columnLower[iColumn]); |
4381 | | } |
4382 | | solver_->initialSolve(); |
4383 | | } |
4384 | | #ifdef COIN_HAS_CLP |
4385 | | { |
4386 | | OsiClpSolverInterface * clpSolver |
4387 | | = dynamic_cast<OsiClpSolverInterface *> (solver_); |
4388 | | if (clpSolver) { |
4389 | | // Possible restore of pivot method |
4390 | | if (savePivotMethod) { |
4391 | | // model may have changed |
4392 | | savePivotMethod->setModel(NULL); |
4393 | | clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod); |
4394 | | delete savePivotMethod; |
4395 | | } |
4396 | | clpSolver->setLargestAway(-1.0); |
4397 | | } |
4398 | | } |
4399 | | #endif |
4400 | | if (fastNodeDepth_ >= 1000000 && !parentModel_) { |
4401 | | // delete object off end |
4402 | | delete object_[numberObjects_]; |
4403 | | fastNodeDepth_ -= 1000000; |
4404 | | } |
4405 | | delete saveSolver; |
4406 | | if (strategy_ && strategy_->preProcessState() > 0) { |
4407 | | // undo preprocessing |
| 1234 | int nOrig = originalSolver->getNumCols(); |