Changeset 520 for branches/devel/Cbc


Ignore:
Timestamp:
Jan 18, 2007 2:37:10 PM (13 years ago)
Author:
forrest
Message:

fix bug in ampl interface and try slp

Location:
branches/devel/Cbc/src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Cbc/src/CbcLinked.cpp

    r503 r520  
    561561              for ( i=rowStart[objectiveRow_];i<rowStart[objectiveRow_+1];i++)
    562562                gradient[column2[i]] = element[i];
    563               const double * columnLower = modelPtr_->columnLower();
    564               const double * columnUpper = modelPtr_->columnUpper();
     563              //const double * columnLower = modelPtr_->columnLower();
     564              //const double * columnUpper = modelPtr_->columnUpper();
    565565              for ( i =0;i<numberObjects_;i++) {
    566566                OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
     
    625625            double * gradient = new double [numberColumns2];
    626626            int * column = new int[numberColumns2];
    627             const double * columnLower = modelPtr_->columnLower();
    628             const double * columnUpper = modelPtr_->columnUpper();
     627            //const double * columnLower = modelPtr_->columnLower();
     628            //const double * columnUpper = modelPtr_->columnUpper();
    629629            for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
    630630              int iRow = rowNonLinear_[iNon];
     
    14221422  }
    14231423  delete [] which;
     1424}
     1425/* Solves nonlinear problem from CoinModel using SLP - may be used as crash
     1426   for other algorithms when number of iterations small.
     1427   Also exits if all problematical variables are changing
     1428   less than deltaTolerance
     1429   Returns solution array
     1430*/
     1431double *
     1432OsiSolverLink::nonlinearSLP(int numberPasses,double deltaTolerance)
     1433{
     1434  if (!coinModel_.numberRows()) {
     1435    printf("Model not set up or nonlinear arrays not created!\n");
     1436    return NULL;
     1437  }
     1438  // first check and set up arrays
     1439  int numberColumns = coinModel_.numberColumns();
     1440  int numberRows = coinModel_.numberRows();
     1441  char * markNonlinear = new char [numberColumns+numberRows];
     1442  CoinZeroN(markNonlinear,numberColumns+numberRows);
     1443  // List of nonlinear entries
     1444  int * listNonLinearColumn = new int[numberColumns];
     1445  // List of nonlinear constraints
     1446  int * whichRow = new int [numberRows];
     1447  CoinZeroN(whichRow,numberRows);
     1448  int numberNonLinearColumns=0;
     1449  int iColumn;
     1450  CoinModel coinModel = coinModel_;
     1451  //const CoinModelHash * stringArray = coinModel.stringArray();
     1452  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1453    CoinModelLink triple=coinModel.firstInColumn(iColumn);
     1454    bool linear=true;
     1455    int n=0;
     1456    // See if nonlinear objective
     1457    const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
     1458    if (strcmp(expr,"Numeric")) {
     1459      linear=false;
     1460      // try and see which columns
     1461      assert (strlen(expr)<20000);
     1462      char temp[20000];
     1463      strcpy(temp,expr);
     1464      char * pos = temp;
     1465      bool ifFirst=true;
     1466      double linearTerm=0.0;
     1467      while (*pos) {
     1468        double value;
     1469        int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
     1470        // must be column unless first when may be linear term
     1471        if (jColumn>=0) {
     1472          markNonlinear[jColumn]=1;
     1473        } else if (jColumn==-2) {
     1474          linearTerm = value;
     1475        } else {
     1476          printf("bad nonlinear term %s\n",temp);
     1477          abort();
     1478        }
     1479        ifFirst=false;
     1480      }
     1481    }
     1482    while (triple.row()>=0) {
     1483      int iRow = triple.row();
     1484      const char * expr = coinModel.getElementAsString(iRow,iColumn);
     1485      if (strcmp(expr,"Numeric")) {
     1486        linear=false;
     1487        whichRow[iRow]++;
     1488        // try and see which columns
     1489        assert (strlen(expr)<20000);
     1490        char temp[20000];
     1491        strcpy(temp,expr);
     1492        char * pos = temp;
     1493        bool ifFirst=true;
     1494        double linearTerm=0.0;
     1495        while (*pos) {
     1496          double value;
     1497          int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
     1498          // must be column unless first when may be linear term
     1499          if (jColumn>=0) {
     1500            markNonlinear[jColumn]=1;
     1501          } else if (jColumn==-2) {
     1502            linearTerm = value;
     1503          } else {
     1504            printf("bad nonlinear term %s\n",temp);
     1505            abort();
     1506          }
     1507          ifFirst=false;
     1508        }
     1509      }
     1510      triple=coinModel.next(triple);
     1511      n++;
     1512    }
     1513    if (!linear) {
     1514      markNonlinear[iColumn]=1;
     1515    }
     1516  }
     1517  //int xxxx[]={3,2,0,4,3,0};
     1518  //double initialSolution[6];
     1519  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1520    if (markNonlinear[iColumn]) {
     1521      // put in something
     1522      double lower = coinModel.columnLower(iColumn);
     1523      double upper = CoinMin(coinModel.columnUpper(iColumn),lower+1000.0);
     1524      coinModel.associateElement(coinModel.columnName(iColumn),0.5*(lower+upper));
     1525      //coinModel.associateElement(coinModel.columnName(iColumn),xxxx[iColumn]);
     1526      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
     1527      //initialSolution[iColumn]=xxxx[iColumn];
     1528    }
     1529  }
     1530  // if nothing just solve
     1531  if (!numberNonLinearColumns) {
     1532    delete [] listNonLinearColumn;
     1533    delete [] whichRow;
     1534    delete [] markNonlinear;
     1535    ClpSimplex tempModel;
     1536    tempModel.loadProblem(coinModel,true);
     1537    tempModel.initialSolve();
     1538    double * solution = CoinCopyOfArray(tempModel.getColSolution(),numberColumns);
     1539    return solution;
     1540  }
     1541  // Create artificials
     1542  ClpSimplex tempModel;
     1543  tempModel.loadProblem(coinModel,true);
     1544  const double * rowLower = tempModel.rowLower();
     1545  const double * rowUpper = tempModel.rowUpper();
     1546  bool takeAll=false;
     1547  int iRow;
     1548  int numberArtificials=0;
     1549  for (iRow=0;iRow<numberRows;iRow++) {
     1550    if (whichRow[iRow]||takeAll) {
     1551      if (rowLower[iRow]>-1.0e30)
     1552        numberArtificials++;
     1553      if (rowUpper[iRow]<1.0e30)
     1554        numberArtificials++;
     1555    }
     1556  }
     1557  CoinBigIndex * startArtificial = new CoinBigIndex [numberArtificials+1];
     1558  int * rowArtificial = new int [numberArtificials];
     1559  double * elementArtificial = new double [numberArtificials];
     1560  double * objectiveArtificial = new double [numberArtificials];
     1561  numberArtificials=0;
     1562  startArtificial[0]=0;
     1563  double artificialCost =1.0e9;
     1564  for (iRow=0;iRow<numberRows;iRow++) {
     1565    if (whichRow[iRow]||takeAll) {
     1566      if (rowLower[iRow]>-1.0e30) {
     1567        rowArtificial[numberArtificials]=iRow;
     1568        elementArtificial[numberArtificials]=1.0;
     1569        objectiveArtificial[numberArtificials]=artificialCost;
     1570        numberArtificials++;
     1571        startArtificial[numberArtificials]=numberArtificials;
     1572      }
     1573      if (rowUpper[iRow]<1.0e30) {
     1574        rowArtificial[numberArtificials]=iRow;
     1575        elementArtificial[numberArtificials]=-1.0;
     1576        objectiveArtificial[numberArtificials]=artificialCost;
     1577        numberArtificials++;
     1578        startArtificial[numberArtificials]=numberArtificials;
     1579      }
     1580    }
     1581  }
     1582  // Get first solution
     1583  int numberColumnsSmall=numberColumns;
     1584  ClpSimplex model;
     1585  model.loadProblem(coinModel,true);
     1586  model.addColumns(numberArtificials,NULL,NULL,objectiveArtificial,
     1587                       startArtificial,rowArtificial,elementArtificial);
     1588  double * columnLower = model.columnLower();
     1589  double * columnUpper = model.columnUpper();
     1590  double * trueLower = new double[numberNonLinearColumns];
     1591  double * trueUpper = new double[numberNonLinearColumns];
     1592  int jNon;
     1593  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1594    iColumn=listNonLinearColumn[jNon];
     1595    trueLower[jNon]=columnLower[iColumn];
     1596    trueUpper[jNon]=columnUpper[iColumn];
     1597    //columnLower[iColumn]=initialSolution[iColumn];
     1598    //columnUpper[iColumn]=initialSolution[iColumn];
     1599  }
     1600  model.initialSolve();
     1601  model.writeMps("bad.mps");
     1602  // redo number of columns
     1603  numberColumns = model.numberColumns();
     1604  int * last[3];
     1605  double * solution = model.primalColumnSolution();
     1606 
     1607  double * trust = new double[numberNonLinearColumns];
     1608  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1609    iColumn=listNonLinearColumn[jNon];
     1610    trust[jNon]=0.5;
     1611    if (solution[iColumn]<trueLower[jNon])
     1612      solution[iColumn]=trueLower[jNon];
     1613    else if (solution[iColumn]>trueUpper[jNon])
     1614      solution[iColumn]=trueUpper[jNon];
     1615  }
     1616  int iPass;
     1617  double lastObjective=1.0e31;
     1618  double * saveSolution = new double [numberColumns];
     1619  double * saveRowSolution = new double [numberRows];
     1620  memset(saveRowSolution,0,numberRows*sizeof(double));
     1621  double * savePi = new double [numberRows];
     1622  double * safeSolution = new double [numberColumns];
     1623  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
     1624  double targetDrop=1.0e31;
     1625  //double objectiveOffset;
     1626  //model.getDblParam(ClpObjOffset,objectiveOffset);
     1627  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
     1628  for (iPass=0;iPass<3;iPass++) {
     1629    last[iPass]=new int[numberNonLinearColumns];
     1630    for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     1631      last[iPass][jNon]=0;
     1632  }
     1633  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
     1634  int goodMove=-2;
     1635  char * statusCheck = new char[numberColumns];
     1636  double * changeRegion = new double [numberColumns];
     1637  int logLevel=63;
     1638  double dualTolerance = model.dualTolerance();
     1639  double primalTolerance = model.primalTolerance();
     1640  int lastGoodMove=1;
     1641  for (iPass=0;iPass<numberPasses;iPass++) {
     1642    lastGoodMove=goodMove;
     1643    columnLower = model.columnLower();
     1644    columnUpper = model.columnUpper();
     1645    solution = model.primalColumnSolution();
     1646    double * rowActivity = model.primalRowSolution();
     1647    // redo objective
     1648    ClpSimplex tempModel;
     1649    // load new values
     1650    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1651      iColumn=listNonLinearColumn[jNon];
     1652      coinModel.associateElement(coinModel.columnName(iColumn),solution[iColumn]);
     1653    }
     1654    tempModel.loadProblem(coinModel);
     1655    double objectiveOffset;
     1656    tempModel.getDblParam(ClpObjOffset,objectiveOffset);
     1657    double objValue=-objectiveOffset;
     1658    const double * objective = tempModel.objective();
     1659    for (iColumn=0;iColumn<numberColumnsSmall;iColumn++)
     1660      objValue += solution[iColumn]*objective[iColumn];
     1661    double * rowActivity2 = tempModel.primalRowSolution();
     1662    const double * rowLower2 = tempModel.rowLower();
     1663    const double * rowUpper2 = tempModel.rowUpper();
     1664    memset(rowActivity2,0,numberRows*sizeof(double));
     1665    tempModel.times(1.0,solution,rowActivity2);
     1666    for (iRow=0;iRow<numberRows;iRow++) {
     1667      if (rowActivity2[iRow]<rowLower2[iRow]-primalTolerance)
     1668        objValue += (rowLower2[iRow]-rowActivity2[iRow]-primalTolerance)*artificialCost;
     1669      else if (rowActivity2[iRow]>rowUpper2[iRow]+primalTolerance)
     1670        objValue -= (rowUpper2[iRow]-rowActivity2[iRow]+primalTolerance)*artificialCost;
     1671    }
     1672    double theta=-1.0;
     1673    double maxTheta=COIN_DBL_MAX;
     1674    if (objValue<=lastObjective+1.0e-15*fabs(lastObjective)||!iPass)
     1675      goodMove=1;
     1676    else
     1677      goodMove=-1;
     1678    //maxTheta=1.0;
     1679    if (iPass) {
     1680      int jNon=0;
     1681      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1682        changeRegion[iColumn]=solution[iColumn]-saveSolution[iColumn];
     1683        double alpha = changeRegion[iColumn];
     1684        double oldValue = saveSolution[iColumn];
     1685        if (markNonlinear[iColumn]==0) {
     1686          // linear
     1687          if (alpha<-1.0e-15) {
     1688            // variable going towards lower bound
     1689            double bound = columnLower[iColumn];
     1690            oldValue -= bound;
     1691            if (oldValue+maxTheta*alpha<0.0) {
     1692              maxTheta = CoinMax(0.0,oldValue/(-alpha));
     1693            }
     1694          } else if (alpha>1.0e-15) {
     1695            // variable going towards upper bound
     1696            double bound = columnUpper[iColumn];
     1697            oldValue = bound-oldValue;
     1698            if (oldValue-maxTheta*alpha<0.0) {
     1699              maxTheta = CoinMax(0.0,oldValue/alpha);
     1700            }
     1701          }
     1702        } else {
     1703          // nonlinear
     1704          if (alpha<-1.0e-15) {
     1705            // variable going towards lower bound
     1706            double bound = trueLower[jNon];
     1707            oldValue -= bound;
     1708            if (oldValue+maxTheta*alpha<0.0) {
     1709              maxTheta = CoinMax(0.0,oldValue/(-alpha));
     1710            }
     1711          } else if (alpha>1.0e-15) {
     1712            // variable going towards upper bound
     1713            double bound = trueUpper[jNon];
     1714            oldValue = bound-oldValue;
     1715            if (oldValue-maxTheta*alpha<0.0) {
     1716              maxTheta = CoinMax(0.0,oldValue/alpha);
     1717            }
     1718          }
     1719          jNon++;
     1720        }
     1721      }
     1722      // make sure both accurate
     1723      memset(rowActivity,0,numberRows*sizeof(double));
     1724      model.times(1.0,solution,rowActivity);
     1725      memset(saveRowSolution,0,numberRows*sizeof(double));
     1726      model.times(1.0,saveSolution,saveRowSolution);
     1727      for (int iRow=0;iRow<numberRows;iRow++) {
     1728        double alpha =rowActivity[iRow]-saveRowSolution[iRow];
     1729        double oldValue = saveRowSolution[iRow];
     1730        if (alpha<-1.0e-15) {
     1731          // variable going towards lower bound
     1732          double bound = rowLower[iRow];
     1733          oldValue -= bound;
     1734          if (oldValue+maxTheta*alpha<0.0) {
     1735            maxTheta = CoinMax(0.0,oldValue/(-alpha));
     1736          }
     1737        } else if (alpha>1.0e-15) {
     1738          // variable going towards upper bound
     1739          double bound = rowUpper[iRow];
     1740          oldValue = bound-oldValue;
     1741          if (oldValue-maxTheta*alpha<0.0) {
     1742            maxTheta = CoinMax(0.0,oldValue/alpha);
     1743          }
     1744        }
     1745      }
     1746    } else {
     1747      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1748        changeRegion[iColumn]=0.0;
     1749        saveSolution[iColumn]=solution[iColumn];
     1750      }
     1751      memcpy(saveRowSolution,rowActivity,numberRows*sizeof(double));
     1752    }
     1753    if (goodMove>=0) {
     1754      //theta = CoinMin(theta2,maxTheta);
     1755      theta = maxTheta;
     1756      if (theta>0.0&&theta<=1.0) {
     1757        // update solution
     1758        double lambda = 1.0-theta;
     1759        for (iColumn=0;iColumn<numberColumns;iColumn++)
     1760          solution[iColumn] = lambda * saveSolution[iColumn]
     1761            + theta * solution[iColumn];
     1762        memset(rowActivity,0,numberRows*sizeof(double));
     1763        model.times(1.0,solution,rowActivity);
     1764        if (lambda>0.999) {
     1765          memcpy(model.dualRowSolution(),savePi,numberRows*sizeof(double));
     1766          memcpy(model.statusArray(),saveStatus,numberRows+numberColumns);
     1767        }
     1768        // redo rowActivity
     1769        memset(rowActivity,0,numberRows*sizeof(double));
     1770        model.times(1.0,solution,rowActivity);
     1771      }
     1772    }
     1773    // load new values
     1774    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1775      iColumn=listNonLinearColumn[jNon];
     1776      coinModel.associateElement(coinModel.columnName(iColumn),solution[iColumn]);
     1777    }
     1778    double * sol2 = CoinCopyOfArray(model.primalColumnSolution(),numberColumns);
     1779    unsigned char * status2 = CoinCopyOfArray(model.statusArray(),numberColumns);
     1780    model.loadProblem(coinModel);
     1781    model.addColumns(numberArtificials,NULL,NULL,objectiveArtificial,
     1782                     startArtificial,rowArtificial,elementArtificial);
     1783    memcpy(model.primalColumnSolution(),sol2,numberColumns*sizeof(double));
     1784    memcpy(model.statusArray(),status2,numberColumns);
     1785    delete [] sol2;
     1786    delete [] status2;
     1787    columnLower = model.columnLower();
     1788    columnUpper = model.columnUpper();
     1789    solution = model.primalColumnSolution();
     1790    rowActivity = model.primalRowSolution();
     1791    int * temp=last[2];
     1792    last[2]=last[1];
     1793    last[1]=last[0];
     1794    last[0]=temp;
     1795    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1796      iColumn=listNonLinearColumn[jNon];
     1797      double change = solution[iColumn]-saveSolution[iColumn];
     1798      if (change<-1.0e-5) {
     1799        if (fabs(change+trust[jNon])<1.0e-5)
     1800          temp[jNon]=-1;
     1801        else
     1802          temp[jNon]=-2;
     1803      } else if(change>1.0e-5) {
     1804        if (fabs(change-trust[jNon])<1.0e-5)
     1805          temp[jNon]=1;
     1806        else
     1807          temp[jNon]=2;
     1808      } else {
     1809        temp[jNon]=0;
     1810      }
     1811    }
     1812    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
     1813    double maxDelta=0.0;
     1814    if (goodMove>=0) {
     1815      if (objValue<=lastObjective+1.0e-15*fabs(lastObjective))
     1816        goodMove=1;
     1817      else
     1818        goodMove=0;
     1819    } else {
     1820      maxDelta=1.0e10;
     1821    }
     1822    double maxGap=0.0;
     1823    int numberSmaller=0;
     1824    int numberSmaller2=0;
     1825    int numberLarger=0;
     1826    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1827      iColumn=listNonLinearColumn[jNon];
     1828      maxDelta = CoinMax(maxDelta,
     1829                     fabs(solution[iColumn]-saveSolution[iColumn]));
     1830      if (goodMove>0) {
     1831        if (last[0][jNon]*last[1][jNon]<0) {
     1832          // halve
     1833          trust[jNon] *= 0.5;
     1834          numberSmaller2++;
     1835        } else {
     1836          if (last[0][jNon]==last[1][jNon]&&
     1837              last[0][jNon]==last[2][jNon])
     1838            trust[jNon] = CoinMin(1.5*trust[jNon],1.0e6);
     1839          numberLarger++;
     1840        }
     1841      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
     1842        trust[jNon] *= 0.2;
     1843        numberSmaller++;
     1844      }
     1845      maxGap = CoinMax(maxGap,trust[jNon]);
     1846    }
     1847#ifdef CLP_DEBUG
     1848    if (logLevel&32)
     1849      std::cout<<"largest gap is "<<maxGap<<" "
     1850               <<numberSmaller+numberSmaller2<<" reduced ("
     1851               <<numberSmaller<<" badMove ), "
     1852               <<numberLarger<<" increased"<<std::endl;
     1853#endif
     1854    if (iPass>10000) {
     1855      for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     1856        trust[jNon] *=0.0001;
     1857    }
     1858    printf("last good %d goodMove %d\n",lastGoodMove,goodMove);
     1859    if (goodMove>0) {
     1860      double drop = lastObjective-objValue;
     1861      printf("Pass %d, objective %g - drop %g maxDelta %g\n",iPass,objValue,drop,maxDelta);
     1862      if (iPass>20&&drop<1.0e-12*fabs(objValue)&&lastGoodMove>0)
     1863        drop=0.999e-4; // so will exit
     1864      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&theta<0.99999&&lastGoodMove>0) {
     1865        if (logLevel>1)
     1866          std::cout<<"Exiting as maxDelta < tolerance and small drop"<<std::endl;
     1867        break;
     1868      }
     1869    } else if (!numberSmaller&&iPass>1) {
     1870      if (logLevel>1)
     1871          std::cout<<"Exiting as all gaps small"<<std::endl;
     1872        break;
     1873    }
     1874    if (!iPass)
     1875      goodMove=1;
     1876    targetDrop=0.0;
     1877    double * r = model.dualColumnSolution();
     1878    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1879      iColumn=listNonLinearColumn[jNon];
     1880      columnLower[iColumn]=CoinMax(solution[iColumn]
     1881                                   -trust[jNon],
     1882                                   trueLower[jNon]);
     1883      columnUpper[iColumn]=CoinMin(solution[iColumn]
     1884                                   +trust[jNon],
     1885                                   trueUpper[jNon]);
     1886    }
     1887    if (iPass) {
     1888      // get reduced costs
     1889      model.matrix()->transposeTimes(savePi,
     1890                                     model.dualColumnSolution());
     1891      const double * objective = model.objective();
     1892      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1893        iColumn=listNonLinearColumn[jNon];
     1894        double dj = objective[iColumn]-r[iColumn];
     1895        r[iColumn]=dj;
     1896        if (dj<-dualTolerance)
     1897          targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
     1898        else if (dj>dualTolerance)
     1899          targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
     1900      }
     1901    } else {
     1902      memset(r,0,numberColumns*sizeof(double));
     1903    }
     1904#if 0
     1905    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1906      iColumn=listNonLinearColumn[jNon];
     1907      if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
     1908        columnLower[iColumn]=CoinMax(solution[iColumn],
     1909                                     trueLower[jNon]);
     1910        columnUpper[iColumn]=CoinMin(solution[iColumn]
     1911                                     +trust[jNon],
     1912                                     trueUpper[jNon]);
     1913      } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
     1914        columnLower[iColumn]=CoinMax(solution[iColumn]
     1915                                     -trust[jNon],
     1916                                     trueLower[jNon]);
     1917        columnUpper[iColumn]=CoinMin(solution[iColumn],
     1918                                     trueUpper[jNon]);
     1919      } else {
     1920        columnLower[iColumn]=CoinMax(solution[iColumn]
     1921                                     -trust[jNon],
     1922                                     trueLower[jNon]);
     1923        columnUpper[iColumn]=CoinMin(solution[iColumn]
     1924                                     +trust[jNon],
     1925                                     trueUpper[jNon]);
     1926      }
     1927    }
     1928#endif
     1929    if (goodMove>0) {
     1930      memcpy(saveSolution,solution,numberColumns*sizeof(double));
     1931      memcpy(saveRowSolution,rowActivity,numberRows*sizeof(double));
     1932      memcpy(savePi,model.dualRowSolution(),numberRows*sizeof(double));
     1933      memcpy(saveStatus,model.statusArray(),numberRows+numberColumns);
     1934     
     1935#ifdef CLP_DEBUG
     1936      if (logLevel&32)
     1937        std::cout<<"Pass - "<<iPass
     1938                 <<", target drop is "<<targetDrop
     1939                 <<std::endl;
     1940#endif
     1941      lastObjective = objValue;
     1942      if (targetDrop<CoinMax(1.0e-8,CoinMin(1.0e-6,1.0e-6*fabs(objValue)))&&lastGoodMove&&iPass>3) {
     1943        if (logLevel>1)
     1944          printf("Exiting on target drop %g\n",targetDrop);
     1945        break;
     1946      }
     1947#ifdef CLP_DEBUG
     1948      {
     1949        double * r = model.dualColumnSolution();
     1950        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1951          iColumn=listNonLinearColumn[jNon];
     1952          if (logLevel&32)
     1953            printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
     1954                   jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
     1955                   r[iColumn],statusCheck[iColumn],columnLower[iColumn],
     1956                   columnUpper[iColumn]);
     1957        }
     1958      }
     1959#endif
     1960      model.scaling(false);
     1961      model.primal(1);
     1962      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1963        iColumn=listNonLinearColumn[jNon];
     1964        printf("%d bounds etc %g %g %g\n",iColumn, columnLower[iColumn],solution[iColumn],columnUpper[iColumn]);
     1965      }
     1966      char temp[20];
     1967      sprintf(temp,"pass%d.mps",iPass);
     1968      model.writeMps(temp);
     1969#ifdef CLP_DEBUG
     1970      if (model.status()) {
     1971        model.writeMps("xx.mps");
     1972      }
     1973#endif
     1974      if (model.status()==1) {
     1975        // not feasible ! - backtrack and exit
     1976        // use safe solution
     1977        memcpy(solution,safeSolution,numberColumns*sizeof(double));
     1978        memcpy(saveSolution,solution,numberColumns*sizeof(double));
     1979        memset(rowActivity,0,numberRows*sizeof(double));
     1980        model.times(1.0,solution,rowActivity);
     1981        memcpy(saveRowSolution,rowActivity,numberRows*sizeof(double));
     1982        memcpy(model.dualRowSolution(),savePi,numberRows*sizeof(double));
     1983        memcpy(model.statusArray(),saveStatus,numberRows+numberColumns);
     1984        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     1985          iColumn=listNonLinearColumn[jNon];
     1986          columnLower[iColumn]=CoinMax(solution[iColumn]
     1987                                       -trust[jNon],
     1988                                       trueLower[jNon]);
     1989          columnUpper[iColumn]=CoinMin(solution[iColumn]
     1990                                       +trust[jNon],
     1991                                       trueUpper[jNon]);
     1992        }
     1993        break;
     1994      } else {
     1995        // save in case problems
     1996        memcpy(safeSolution,solution,numberColumns*sizeof(double));
     1997      }
     1998      goodMove=1;
     1999    } else {
     2000      // bad pass - restore solution
     2001#ifdef CLP_DEBUG
     2002      if (logLevel&32)
     2003        printf("Backtracking\n");
     2004#endif
     2005      // load old values
     2006      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     2007        iColumn=listNonLinearColumn[jNon];
     2008        coinModel.associateElement(coinModel.columnName(iColumn),saveSolution[iColumn]);
     2009      }
     2010      model.loadProblem(coinModel);
     2011      model.addColumns(numberArtificials,NULL,NULL,objectiveArtificial,
     2012                     startArtificial,rowArtificial,elementArtificial);
     2013      solution = model.primalColumnSolution();
     2014      rowActivity = model.primalRowSolution();
     2015      memcpy(solution,saveSolution,numberColumns*sizeof(double));
     2016      memcpy(rowActivity,saveRowSolution,numberRows*sizeof(double));
     2017      memcpy(model.dualRowSolution(),savePi,numberRows*sizeof(double));
     2018      memcpy(model.statusArray(),saveStatus,numberRows+numberColumns);
     2019      columnLower = model.columnLower();
     2020      columnUpper = model.columnUpper();
     2021      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     2022        iColumn=listNonLinearColumn[jNon];
     2023        columnLower[iColumn]=solution[iColumn];
     2024        columnUpper[iColumn]=solution[iColumn];
     2025      }
     2026      model.primal(1);
     2027      model.writeMps("xx.mps");
     2028      iPass--;
     2029      goodMove=-1;
     2030    }
     2031  }
     2032  // restore solution
     2033  memcpy(solution,saveSolution,numberColumns*sizeof(double));
     2034  delete [] statusCheck;
     2035  delete [] savePi;
     2036  delete [] saveStatus;
     2037  // load new values
     2038  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     2039    iColumn=listNonLinearColumn[jNon];
     2040    coinModel.associateElement(coinModel.columnName(iColumn),solution[iColumn]);
     2041  }
     2042  double * sol2 = CoinCopyOfArray(model.primalColumnSolution(),numberColumns);
     2043  unsigned char * status2 = CoinCopyOfArray(model.statusArray(),numberColumns);
     2044  model.loadProblem(coinModel);
     2045  model.addColumns(numberArtificials,NULL,NULL,objectiveArtificial,
     2046                   startArtificial,rowArtificial,elementArtificial);
     2047  memcpy(model.primalColumnSolution(),sol2,numberColumns*sizeof(double));
     2048  memcpy(model.statusArray(),status2,numberColumns);
     2049  delete [] sol2;
     2050  delete [] status2;
     2051  columnLower = model.columnLower();
     2052  columnUpper = model.columnUpper();
     2053  solution = model.primalColumnSolution();
     2054  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     2055    iColumn=listNonLinearColumn[jNon];
     2056    columnLower[iColumn]=CoinMax(solution[iColumn],
     2057                                 trueLower[jNon]);
     2058    columnUpper[iColumn]=CoinMin(solution[iColumn],
     2059                                 trueUpper[jNon]);
     2060  }
     2061  model.primal(1);
     2062  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     2063    iColumn=listNonLinearColumn[jNon];
     2064    columnLower[iColumn]= trueLower[jNon];
     2065    columnUpper[iColumn]= trueUpper[jNon];
     2066  }
     2067  delete [] saveSolution;
     2068  delete [] safeSolution;
     2069  delete [] saveRowSolution;
     2070  for (iPass=0;iPass<3;iPass++)
     2071    delete [] last[iPass];
     2072  delete [] trust;
     2073  delete [] trueUpper;
     2074  delete [] trueLower;
     2075  delete [] changeRegion;
     2076  delete [] startArtificial;
     2077  delete [] rowArtificial;
     2078  delete [] elementArtificial;
     2079  delete [] objectiveArtificial;
     2080  delete [] listNonLinearColumn;
     2081  delete [] whichRow;
     2082  delete [] markNonlinear;
     2083  return CoinCopyOfArray(solution,coinModel.numberColumns());
    14242084}
    14252085// Analyze constraints to see which are convex (quadratic)
  • branches/devel/Cbc/src/CbcLinked.hpp

    r501 r520  
    4646  */
    4747  virtual int fathom(bool allFixed) {return 0;};
     48  /** Solves nonlinear problem from CoinModel using SLP - may be used as crash
     49      for other algorithms when number of iterations small.
     50      Also exits if all problematical variables are changing
     51      less than deltaTolerance
     52      Returns solution array
     53  */
     54  double * nonlinearSLP(int numberPasses,double deltaTolerance);
    4855  //@}
    4956 
     
    152159  inline void setCbcModel(CbcModel * model)
    153160  { cbcModel_=model;};
     161  /// Return CoinModel
     162  inline const CoinModel * coinModel() const
     163  { return &coinModel_;};
    154164  //@}
    155165 
     
    162172  //@{
    163173  /// Do real work of initialize
    164   void initialize(ClpSimplex * & solver, OsiObject ** & object) const;
     174  //void initialize(ClpSimplex * & solver, OsiObject ** & object) const;
    165175  /// Do real work of delete
    166176  void gutsOfDestructor(bool justNullify=false);
  • branches/devel/Cbc/src/Cbc_ampl.cpp

    r509 r520  
    715715  int numberBinary=-1;
    716716  int numberIntegers=-1;
     717  int numberAllNonLinearBoth=0;
     718  int numberIntegerNonLinearBoth=0;
     719  int numberAllNonLinearConstraints=0;
     720  int numberIntegerNonLinearConstraints=0;
     721  int numberAllNonLinearObjective=0;
     722  int numberIntegerNonLinearObjective=0;
    717723  double * primalSolution=NULL;
    718724  double direction=1.0;
     
    923929    numberElements=nzc;;
    924930    numberBinary=nbv;
    925     numberIntegers=niv+nlvci;
     931    numberIntegers=niv;
     932    numberAllNonLinearBoth=nlvb;
     933    numberIntegerNonLinearBoth=nlvbi;
     934    numberAllNonLinearConstraints=nlvc;
     935    numberIntegerNonLinearConstraints=nlvci;
     936    numberAllNonLinearObjective=nlvo;
     937    numberIntegerNonLinearObjective=nlvoi;
    926938    /* say we want primal solution */
    927939    want_xpi0=1;
     
    951963        i<numberColumns;i++) {
    952964      setColumnIsInteger(i,true);;
     965  }
     966  // and non linear
     967  for (i=numberAllNonLinearBoth-numberIntegerNonLinearBoth;
     968       i<numberAllNonLinearBoth;i++) {
     969    setColumnIsInteger(i,true);;
     970  }
     971  for (i=numberAllNonLinearConstraints-numberIntegerNonLinearConstraints;
     972       i<numberAllNonLinearConstraints;i++) {
     973    setColumnIsInteger(i,true);;
     974  }
     975  for (i=numberAllNonLinearObjective-numberIntegerNonLinearObjective;
     976       i<numberAllNonLinearObjective;i++) {
     977    setColumnIsInteger(i,true);;
    953978  }
    954979  // do names
  • branches/devel/Cbc/src/CoinSolve.cpp

    r505 r520  
    15881588                solveOptions.setSpecialOption(4,barrierOptions);
    15891589              }
    1590               model2->initialSolve(solveOptions);
     1590              OsiSolverInterface * coinSolver = model.solver();
     1591              OsiSolverLink * linkSolver = dynamic_cast< OsiSolverLink*> (coinSolver);
     1592              if (!linkSolver) {
     1593                model2->initialSolve(solveOptions);
     1594              } else {
     1595                // special solver
     1596                double * solution = linkSolver->nonlinearSLP(slpValue,1.0e-5);
     1597                if (solution) {
     1598                  memcpy(model2->primalColumnSolution(),solution,
     1599                         CoinMin(model2->numberColumns(),linkSolver->coinModel()->numberColumns())*sizeof(double));
     1600                  delete [] solution;
     1601                } else {
     1602                  printf("No nonlinear solution\n");
     1603                }
     1604              }
    15911605              basisHasValues=1;
    15921606              if (dualize) {
Note: See TracChangeset for help on using the changeset viewer.