Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcSimpleIntegerDynamicPseudoCost.cpp

    r1902 r1943  
    2424#include "CoinError.hpp"
    2525#include "CbcSimpleIntegerDynamicPseudoCost.hpp"
     26#ifdef COIN_HAS_CLP
     27#include "OsiClpSolverInterface.hpp"
     28#endif
    2629#ifdef COIN_DEVELOP
    2730typedef struct {
     
    840843    if (!info->hotstartSolution_ && priority_ != -999) {
    841844#ifndef NDEBUG
     845#ifndef SWITCH_VARIABLES
    842846        double nearest = floor(value + 0.5);
    843847        assert (fabs(value - nearest) > info->integerTolerance_);
     848#endif
    844849#endif
    845850    } else if (info->hotstartSolution_) {
     
    10161021            if (fabs(value - nearest) > integerTolerance)
    10171022                unsatisfied++;
     1023#ifdef SWITCH_VARIABLES
     1024            const CbcSwitchingBinary * sObject = dynamic_cast<const CbcSwitchingBinary *> (this);
     1025            if (sObject) {
     1026              int state[3],nBadFixed;
     1027              unsatisfied +=
     1028                sObject->checkAssociatedBounds(solver,solution,0,
     1029                                               state,nBadFixed);
     1030            }
     1031#endif
    10181032        }
    10191033    }
     
    13981412    return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
    13991413}
    1400 
     1414#ifdef SWITCH_VARIABLES
     1415/** Default Constructor
     1416
     1417  Equivalent to an unspecified binary variable.
     1418*/
     1419CbcSwitchingBinary::CbcSwitchingBinary ()
     1420        : CbcSimpleIntegerDynamicPseudoCost(),
     1421          zeroLowerBound_(NULL),
     1422          oneLowerBound_(NULL),
     1423          zeroUpperBound_(NULL),
     1424          oneUpperBound_(NULL),
     1425          otherVariable_(NULL),
     1426          numberOther_(0),
     1427          type_(0)
     1428{
     1429}
     1430
     1431/** Useful constructor
     1432*/
     1433CbcSwitchingBinary::CbcSwitchingBinary (CbcSimpleIntegerDynamicPseudoCost * oldObject,
     1434                                        int nOdd,const int * other, const int * otherRow)
     1435  : CbcSimpleIntegerDynamicPseudoCost(*oldObject),
     1436    zeroLowerBound_(NULL),
     1437    oneLowerBound_(NULL),
     1438    zeroUpperBound_(NULL),
     1439    oneUpperBound_(NULL),
     1440    otherVariable_(NULL),
     1441    numberOther_(0),
     1442    type_(0)
     1443{
     1444  if (nOdd)
     1445    type_=2;
     1446  const CoinPackedMatrix * rowCopy = model_->solver()->getMatrixByRow();
     1447  const int * column = rowCopy->getIndices();
     1448  //const int * rowLength = rowCopy->getVectorLengths();
     1449  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     1450  //const double * rowLower = model_->solver()->getRowLower();
     1451  const double * rowUpper = model_->solver()->getRowUpper();
     1452  const double * columnLower = model_->solver()->getColLower();
     1453  const double * columnUpper = model_->solver()->getColUpper();
     1454  const double * element = rowCopy->getElements();
     1455  int last = other[0];
     1456  int nPair=0;
     1457  int nInGroup=1;
     1458  for (int i=1;i<=nOdd;i++) {
     1459    if (other[i]==last) {
     1460      nInGroup++;
     1461    } else {
     1462      if (nInGroup>2 && model_->logLevel()>2)
     1463        printf("%d in group for column %d - some redundancy\n",
     1464               nInGroup,columnNumber_);
     1465      nPair++;
     1466      last=other[i];
     1467      nInGroup=1;
     1468    }
     1469  }
     1470  zeroLowerBound_ = new double [4*nPair];
     1471  oneLowerBound_ = zeroLowerBound_+nPair;
     1472  zeroUpperBound_ = oneLowerBound_+nPair;
     1473  oneUpperBound_ = zeroUpperBound_+nPair;
     1474  otherVariable_ = new int [nPair];
     1475  numberOther_ = nPair;
     1476  if (nPair>1&&model_->logLevel()>2)
     1477    printf("%d pairs for column %d\n",
     1478               nPair,columnNumber_);
     1479  // Now fill
     1480  last = other[0];
     1481  nPair=0;
     1482  int rows[20];
     1483  rows[0]= otherRow[0];
     1484  nInGroup=1;
     1485  for (int i=1;i<=nOdd;i++) {
     1486    if (other[i]==last) {
     1487      rows[nInGroup++]=otherRow[i];
     1488    } else {
     1489      double newLowerZero=0.0;
     1490      double newUpperZero=COIN_DBL_MAX;
     1491      double newLowerOne=0.0;
     1492      double newUpperOne=COIN_DBL_MAX;
     1493      int cColumn=-1;
     1494      for (int j=0;j<nInGroup;j++) {
     1495        int iRow = rows[j];
     1496        CoinBigIndex k = rowStart[iRow];
     1497        double bValue, cValue;
     1498        if (column[k]==columnNumber_) {
     1499          bValue=element[k];
     1500          cValue=element[k+1];
     1501          cColumn=column[k+1];
     1502        } else {
     1503          bValue=element[k+1];
     1504          cValue=element[k];
     1505          cColumn=column[k];
     1506        }
     1507        if (rowUpper[iRow]) {
     1508          // G row - convert to L
     1509          bValue=-bValue;
     1510          cValue=-cValue;
     1511        }
     1512        if (bValue>0.0) {
     1513          // binary*abs(bValue) <= continuous*abs(cValue);
     1514          newLowerOne = -bValue/cValue;
     1515        } else {
     1516          // binary*abs(bValue) >= continuous*abs(cValue);
     1517          newUpperOne = -bValue/cValue;
     1518          newUpperZero = 0.0;
     1519        }
     1520      }
     1521      zeroLowerBound_[nPair]=newLowerZero;
     1522      oneLowerBound_[nPair]=newLowerOne;
     1523      zeroUpperBound_[nPair]=newUpperZero;
     1524      oneUpperBound_[nPair]=newUpperOne;
     1525      // make current bounds tight
     1526      double newLower = CoinMin(newLowerZero,newLowerOne);
     1527      if (newLower>columnLower[cColumn])
     1528        model_->solver()->setColLower(cColumn,newLower);
     1529      double newUpper = CoinMax(newUpperZero,newUpperOne);
     1530      if (newUpper<columnUpper[cColumn])
     1531        model_->solver()->setColUpper(cColumn,newUpper);
     1532      otherVariable_[nPair++]=cColumn;
     1533      last=other[i];
     1534      rows[0] = otherRow[i];
     1535      nInGroup=1;
     1536    }
     1537  }
     1538}
     1539// Copy constructor
     1540CbcSwitchingBinary::CbcSwitchingBinary ( const CbcSwitchingBinary & rhs)
     1541        : CbcSimpleIntegerDynamicPseudoCost(rhs),
     1542          numberOther_(rhs.numberOther_),
     1543          type_(rhs.type_)
     1544{
     1545  zeroLowerBound_ = CoinCopyOfArray(rhs.zeroLowerBound_,4*numberOther_);
     1546  oneLowerBound_ = zeroLowerBound_+numberOther_;
     1547  zeroUpperBound_ = oneLowerBound_+numberOther_;
     1548  oneUpperBound_ = zeroUpperBound_+numberOther_;
     1549  otherVariable_ = CoinCopyOfArray(rhs.otherVariable_,numberOther_);
     1550}
     1551
     1552// Clone
     1553CbcObject *
     1554CbcSwitchingBinary::clone() const
     1555{
     1556    return new CbcSwitchingBinary(*this);
     1557}
     1558
     1559// Assignment operator
     1560CbcSwitchingBinary &
     1561CbcSwitchingBinary::operator=( const CbcSwitchingBinary & rhs)
     1562{
     1563    if (this != &rhs) {
     1564        CbcSimpleIntegerDynamicPseudoCost::operator=(rhs);
     1565        numberOther_=rhs.numberOther_;
     1566        type_ = rhs.type_;
     1567        delete [] zeroLowerBound_;
     1568        delete [] otherVariable_;
     1569        zeroLowerBound_ = CoinCopyOfArray(rhs.zeroLowerBound_,4*numberOther_);
     1570        oneLowerBound_ = zeroLowerBound_+numberOther_;
     1571        zeroUpperBound_ = oneLowerBound_+numberOther_;
     1572        oneUpperBound_ = zeroUpperBound_+numberOther_;
     1573        otherVariable_ = CoinCopyOfArray(rhs.otherVariable_,numberOther_);
     1574    }
     1575    return *this;
     1576}
     1577
     1578// Destructor
     1579CbcSwitchingBinary::~CbcSwitchingBinary ()
     1580{
     1581  delete [] zeroLowerBound_;
     1582  delete [] otherVariable_;
     1583}
     1584// Add in zero switches
     1585void
     1586CbcSwitchingBinary::addZeroSwitches(int nAdd,const int * columns)
     1587{
     1588  type_ |= 1;
     1589  int nNew = numberOther_+nAdd;
     1590  double * bounds = new double[4*nNew];
     1591  int * other = new int [nNew];
     1592  memcpy(other,otherVariable_,numberOther_*sizeof(int));
     1593  delete [] otherVariable_;
     1594  otherVariable_=other;
     1595  memcpy(bounds,zeroLowerBound_,numberOther_*sizeof(double));
     1596  memcpy(bounds+nNew,oneLowerBound_,numberOther_*sizeof(double));
     1597  memcpy(bounds+2*nNew,zeroUpperBound_,numberOther_*sizeof(double));
     1598  memcpy(bounds+3*nNew,oneUpperBound_,numberOther_*sizeof(double));
     1599  delete [] zeroLowerBound_;
     1600  zeroLowerBound_ = bounds;
     1601  oneLowerBound_ = zeroLowerBound_+nNew;
     1602  zeroUpperBound_ = oneLowerBound_+nNew;
     1603  oneUpperBound_ = zeroUpperBound_+nNew;
     1604  for (int i=0;i<nAdd;i++) {
     1605    zeroLowerBound_[numberOther_]=0.0;
     1606    oneLowerBound_[numberOther_]=0.0;
     1607    zeroUpperBound_[numberOther_]=0.0;
     1608    oneUpperBound_[numberOther_]=COIN_DBL_MAX;
     1609    otherVariable_[numberOther_++]=columns[i];
     1610  }
     1611}
     1612// Same - returns true if contents match(ish)
     1613bool
     1614CbcSwitchingBinary::same(const CbcSwitchingBinary * otherObject) const
     1615{
     1616  bool okay = CbcSimpleIntegerDynamicPseudoCost::same(otherObject);
     1617    return okay;
     1618}
     1619double
     1620CbcSwitchingBinary::infeasibility(const OsiBranchingInformation * info,
     1621        int &preferredWay) const
     1622{
     1623    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
     1624    double * solution = const_cast<double *>(model_->testSolution());
     1625    const double * lower = model_->getCbcColLower();
     1626    const double * upper = model_->getCbcColUpper();
     1627    double saveValue = solution[columnNumber_];
     1628    if (!lower[columnNumber_]&&upper[columnNumber_]==1.0) {
     1629      double integerTolerance =
     1630        model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1631      if (saveValue<integerTolerance) {
     1632        // check others OK
     1633        bool allGood=true;
     1634        double tolerance;
     1635        model_->solver()->getDblParam(OsiPrimalTolerance, tolerance) ;
     1636        for (int i=0;i<numberOther_;i++) {
     1637          int otherColumn = otherVariable_[i];
     1638          double value = solution[otherColumn];
     1639          if (value<zeroLowerBound_[i]-tolerance||
     1640              value>zeroUpperBound_[i]+tolerance)
     1641            allGood=false;
     1642        }
     1643        if (!allGood)
     1644          solution[columnNumber_]=2.0*integerTolerance;
     1645      } else if (saveValue>1.0-integerTolerance) {
     1646        // check others OK
     1647        bool allGood=true;
     1648        double tolerance;
     1649        model_->solver()->getDblParam(OsiPrimalTolerance, tolerance) ;
     1650        for (int i=0;i<numberOther_;i++) {
     1651          int otherColumn = otherVariable_[i];
     1652          double value = solution[otherColumn];
     1653          if (value<oneLowerBound_[i]-tolerance||
     1654              value>oneUpperBound_[i]+tolerance)
     1655            allGood=false;
     1656        }
     1657        if (!allGood)
     1658          solution[columnNumber_]=1.0-2.0*integerTolerance;
     1659      }
     1660    }
     1661    double inf = CbcSimpleIntegerDynamicPseudoCost::infeasibility(info,preferredWay);
     1662    solution[columnNumber_]=saveValue;
     1663    return inf;
     1664}
     1665// Set associated bounds
     1666int
     1667CbcSwitchingBinary::setAssociatedBounds(OsiSolverInterface * solver,
     1668                                        int cleanBasis) const
     1669{
     1670  if (!solver)
     1671    solver = model_->solver();
     1672#ifdef COIN_HAS_CLP
     1673  OsiClpSolverInterface * clpSolver
     1674    = dynamic_cast<OsiClpSolverInterface *> (solver);
     1675  if (cleanBasis!=1)
     1676    clpSolver=NULL;
     1677#endif
     1678  const double * columnLower = solver->getColLower();
     1679  const double * columnUpper = solver->getColUpper();
     1680  int nChanged=0;
     1681  if (!columnUpper[columnNumber_]) {
     1682#ifdef COIN_HAS_CLP
     1683    if (clpSolver)
     1684      clpSolver->setColumnStatus(columnNumber_,ClpSimplex::isFixed);
     1685#endif
     1686    for (int i=0;i<numberOther_;i++) {
     1687      int otherColumn = otherVariable_[i];
     1688      if (zeroLowerBound_[i]>columnLower[otherColumn]) {
     1689        solver->setColLower(otherColumn,zeroLowerBound_[i]);
     1690        nChanged++;
     1691      }
     1692      if (zeroUpperBound_[i]<columnUpper[otherColumn]) {
     1693        solver->setColUpper(otherColumn,zeroUpperBound_[i]);
     1694#ifdef COIN_DEVELOP
     1695        const double * solution = solver->getColSolution();
     1696        double value = solution[otherColumn];
     1697        if (value - zeroUpperBound_[i] > 1.0e-5 && model_->logLevel()>1)
     1698          printf("value for continuous %d %g - above %g - switch %d is %.12g (ub 0)\n",
     1699                 otherColumn, value, zeroUpperBound_[i],columnNumber_,solution[columnNumber_]);
     1700#endif
     1701        nChanged++;
     1702      }
     1703    }
     1704  } else if (columnLower[columnNumber_]==1.0) {
     1705#ifdef COIN_HAS_CLP
     1706    if (clpSolver)
     1707      clpSolver->setColumnStatus(columnNumber_,ClpSimplex::isFixed);
     1708#endif
     1709    for (int i=0;i<numberOther_;i++) {
     1710      int otherColumn = otherVariable_[i];
     1711      if (oneLowerBound_[i]>columnLower[otherColumn]) {
     1712        solver->setColLower(otherColumn,oneLowerBound_[i]);
     1713        nChanged++;
     1714      }
     1715      if (oneUpperBound_[i]<columnUpper[otherColumn]) {
     1716        solver->setColUpper(otherColumn,oneUpperBound_[i]);
     1717        nChanged++;
     1718      }
     1719    }
     1720  } else if (cleanBasis>=2) {
     1721    // if all OK then can fix
     1722    int state[3];
     1723    int nBadFixed;
     1724    const double * solution = solver->getColSolution();
     1725    if (!checkAssociatedBounds(solver,solution,
     1726                               0,state,nBadFixed)) {
     1727      const double *reducedCost = solver->getReducedCost() ;
     1728      double good=true;
     1729      double integerTolerance =
     1730        model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1731      if (solution[columnNumber_]<integerTolerance) {
     1732        if (cleanBasis==2||reducedCost[columnNumber_]>1.0e-6)
     1733          solver->setColUpper(columnNumber_,0.0);
     1734        else
     1735          good=false;
     1736      } else if (solution[columnNumber_]>1.0-integerTolerance) {
     1737        if (cleanBasis==2||reducedCost[columnNumber_]<-1.0e-6)
     1738          solver->setColLower(columnNumber_,1.0);
     1739        else
     1740          good=false;
     1741      }
     1742      if (good)
     1743        nChanged=setAssociatedBounds(solver,0);
     1744    }
     1745  } else {
     1746    // see if any continuous bounds force binary
     1747    for (int i=0;i<numberOther_;i++) {
     1748      int otherColumn = otherVariable_[i];
     1749      if (columnLower[otherColumn]>zeroUpperBound_[i]) {
     1750        // can't be zero
     1751        solver->setColLower(columnNumber_,1.0);
     1752        nChanged++;
     1753      } else if (columnLower[otherColumn]>oneUpperBound_[i]) {
     1754        // can't be one
     1755        solver->setColUpper(columnNumber_,0.0);
     1756        nChanged++;
     1757      }
     1758      if (columnUpper[otherColumn]<zeroLowerBound_[i]) {
     1759        // can't be zero
     1760        solver->setColLower(columnNumber_,1.0);
     1761        nChanged++;
     1762      } else if (columnUpper[otherColumn]<oneLowerBound_[i]) {
     1763        // can't be one
     1764        solver->setColUpper(columnNumber_,0.0);
     1765        nChanged++;
     1766      }
     1767    }
     1768  }
     1769  return nChanged;
     1770}
     1771// Check associated bounds
     1772int
     1773CbcSwitchingBinary::checkAssociatedBounds(const OsiSolverInterface * solver,
     1774                                          const double * solution, int printLevel,
     1775                                          int state[3], int & nBadFixed) const
     1776{
     1777  state[0] = 0;
     1778  int nBad=0;
     1779  if (!solver)
     1780    solver = model_->solver();
     1781  double tolerance;
     1782  solver->getDblParam(OsiPrimalTolerance, tolerance) ;
     1783  const double * columnLower = solver->getColLower();
     1784  const double * columnUpper = solver->getColUpper();
     1785  double integerTolerance =
     1786    model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1787  bool printIt = printLevel>2 && model_->logLevel()>1;
     1788  if (solution[columnNumber_]<integerTolerance) {
     1789    state[0] = -1;
     1790    for (int i=0;i<numberOther_;i++) {
     1791      int otherColumn = otherVariable_[i];
     1792      if (zeroLowerBound_[i]>solution[otherColumn]+tolerance*5.0) {
     1793        nBad++;
     1794        if (columnUpper[columnNumber_]==0.0) {
     1795          nBadFixed++;
     1796          //printIt=true;
     1797        }
     1798        if (printIt)
     1799          printf("switch %d at zero, other %d at %.12g below bound of %.12g\n",
     1800                 columnNumber_,otherColumn,solution[otherColumn],zeroLowerBound_[i]);
     1801      }
     1802      if (zeroUpperBound_[i]<solution[otherColumn]-tolerance*5.0) {
     1803        nBad++;
     1804        if (columnUpper[columnNumber_]==0.0) {
     1805          nBadFixed++;
     1806          //printIt=true;
     1807        }
     1808        if (printIt)
     1809          printf("switch %d at zero, other %d at %.12g above bound of %.12g\n",
     1810                 columnNumber_,otherColumn,solution[otherColumn],zeroUpperBound_[i]);
     1811      }
     1812    }
     1813  } else if (solution[columnNumber_]>1.0-integerTolerance) {
     1814    state[0] = 1;
     1815    for (int i=0;i<numberOther_;i++) {
     1816      int otherColumn = otherVariable_[i];
     1817      if (oneLowerBound_[i]>solution[otherColumn]+tolerance*5.0) {
     1818        nBad++;
     1819        if (columnLower[columnNumber_]==1.0) {
     1820          nBadFixed++;
     1821          //printIt=true;
     1822        }
     1823        if (printIt)
     1824          printf("switch %d at one, other %d at %.12g below bound of %.12g\n",
     1825                 columnNumber_,otherColumn,solution[otherColumn],oneLowerBound_[i]);
     1826      }
     1827      if (oneUpperBound_[i]<solution[otherColumn]-tolerance*5.0) {
     1828        nBad++;
     1829        if (columnLower[columnNumber_]==1.0) {
     1830          nBadFixed++;
     1831          //printIt=true;
     1832        }
     1833        if (printIt)
     1834          printf("switch %d at one, other %d at %.12g above bound of %.12g\n",
     1835                 columnNumber_,otherColumn,solution[otherColumn],oneUpperBound_[i]);
     1836      }
     1837    }
     1838  } else {
     1839    // in between - compute tight variables
     1840    state[1]=0;
     1841    state[2]=0;
     1842    // for now just compute ones away from bounds
     1843    for (int i=0;i<numberOther_;i++) {
     1844      int otherColumn = otherVariable_[i];
     1845      double otherValue = solution[otherColumn];
     1846      if (otherValue>columnLower[otherColumn]+tolerance&&
     1847          otherValue<columnUpper[otherColumn]-tolerance)
     1848        state[1]++;
     1849    }
     1850  }
     1851  return nBad;
     1852}
     1853#endif
Note: See TracChangeset for help on using the changeset viewer.