Changeset 501 for branches/devel/Cbc/src


Ignore:
Timestamp:
Dec 12, 2006 4:22:19 PM (13 years ago)
Author:
forrest
Message:

qudaratic stuff

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

Legend:

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

    r500 r501  
    557557      }
    558558      // ???  - try
     559      // But skip if strong branching
     560      CbcModel * cbcModel = (modelPtr_->maximumIterations()<10000) ? cbcModel_ : NULL;
    559561      if ((specialOptions2_&2)!=0) {
    560562        // If model has stored then add cut (if convex)
    561         if (cbcModel_&&(specialOptions2_&4)!=0&&quadraticModel_) {
     563        if (cbcModel&&(specialOptions2_&4)!=0&&quadraticModel_) {
    562564          int numberGenerators = cbcModel_->numberCutGenerators();
    563565          int iGenerator;
     
    627629          }
    628630        }
     631      } else if (cbcModel&&(specialOptions2_&8)==8) {
     632        // convex and nonlinear in constraints
     633        int numberGenerators = cbcModel_->numberCutGenerators();
     634        int iGenerator;
     635        for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
     636          CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
     637          CglCutGenerator * gen = generator->generator();
     638          CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
     639          if (gen2) {
     640            const double * solution = getColSolution();
     641            const double * rowUpper = getRowUpper();
     642            const double * rowLower = getRowLower();
     643            const double * element = originalRowCopy_->getElements();
     644            const int * column2 = originalRowCopy_->getIndices();
     645            const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
     646            //const int * rowLength = originalRowCopy_->getVectorLengths();
     647            int numberColumns2 = CoinMax(coinModel_.numberColumns(),objectiveVariable_+1);
     648            double * gradient = new double [numberColumns2];
     649            int * column = new int[numberColumns2];
     650            for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
     651              int iRow = rowNonLinear_[iNon];
     652              bool convex = convex_[iNon]>0;
     653              if (!convex_[iNon])
     654                continue; // can't use this row
     655              // add OA cuts
     656              double offset=0.0;
     657              // gradient from bilinear
     658              int i;
     659              CoinZeroN(gradient,numberColumns2);
     660              //const double * objective = modelPtr_->objective();
     661              for ( i=rowStart[iRow];i<rowStart[iRow+1];i++)
     662                gradient[column2[i]] = element[i];
     663              for ( i =startNonLinear_[iNon];i<startNonLinear_[iNon+1];i++) {
     664                OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[whichNonLinear_[i]]);
     665                assert (obj);
     666                int xColumn = obj->xColumn();
     667                int yColumn = obj->yColumn();
     668                if (xColumn!=yColumn) {
     669                  double coefficient = 2.0*obj->coefficient();
     670                  gradient[xColumn] += coefficient*solution[yColumn];
     671                  gradient[yColumn] += coefficient*solution[xColumn];
     672                  offset += coefficient*solution[xColumn]*solution[yColumn];
     673                } else {
     674                  double coefficient = obj->coefficient();
     675                  gradient[xColumn] += 2.0*coefficient*solution[yColumn];
     676                  offset += coefficient*solution[xColumn]*solution[yColumn];
     677                }
     678              }
     679              // assume convex
     680              double rhs = 0.0;
     681              int n=0;
     682              for (int i=0;i<numberColumns2;i++) {
     683                double value = gradient[i];
     684                if (fabs(value)>1.0e-12) {
     685                  gradient[n]=value;
     686                  rhs += value*solution[i];
     687                  column[n++]=i;
     688                }
     689              }
     690              if (iRow==objectiveRow_) {
     691                gradient[n]=-1.0;
     692                assert (objectiveVariable_>=0);
     693                rhs -= solution[objectiveVariable_];
     694                column[n++]=objectiveVariable_;
     695                assert (convex);
     696              } else if (convex) {
     697                offset += rowUpper[iRow];
     698              } else if (!convex) {
     699                offset += rowLower[iRow];
     700              }
     701              if (convex&&rhs>offset+1.0e-5)
     702                gen2->addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
     703              else if (!convex&&rhs<offset-1.0e-5)
     704                gen2->addCut(offset-1.0e-7,COIN_DBL_MAX,n,column,gradient);
     705            }
     706            delete [] gradient;
     707            delete [] column;
     708            break;
     709          }
     710        }
    629711      }
    630712    }
     
    826908  int * which = new int[numberColumns];
    827909  numberVariables_=0;
    828   specialOptions2_=0;
     910  //specialOptions2_=0;
    829911  int iColumn;
    830912  int numberErrors=0;
     
    10241106      elementQuadratic = new double [numberQuadratic];
    10251107      numberQuadratic=0;
     1108    }
     1109    if (quadraticObjective||((specialOptions2_&8)!=0&&saveNBi)) {
    10261110      // add in objective as constraint
    10271111      objectiveVariable_= numberColumns;
     
    12451329                                              numberRows2,whichRows,
    12461330                                              numberColumns,whichColumns);
     1331      delete [] whichColumns;
     1332      numberNonLinearRows_=0;
     1333      CoinZeroN(whichRows,numberRows2);
     1334      for ( i =0;i<numberObjects_;i++) {
     1335        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
     1336        if (obj) {
     1337          int xyRow = obj->xyRow();
     1338          assert (xyRow>=0&&xyRow<numberRows2); // even if obj we should move
     1339          whichRows[xyRow]++;
     1340        }
     1341      }
     1342      int * pos = new int [numberRows2];
     1343      int n=0;
     1344      for (i=0;i<numberRows2;i++) {
     1345        if (whichRows[i]) {
     1346          pos[numberNonLinearRows_]=n;
     1347          n+=whichRows[i];
     1348          whichRows[i]=numberNonLinearRows_;
     1349          numberNonLinearRows_++;
     1350        } else {
     1351          whichRows[i]=-1;
     1352        }
     1353      }
     1354      startNonLinear_ = new int [numberNonLinearRows_+1];
     1355      memcpy(startNonLinear_,pos,numberNonLinearRows_*sizeof(int));
     1356      startNonLinear_[numberNonLinearRows_]=n;
     1357      rowNonLinear_ = new int [numberNonLinearRows_];
     1358      convex_ = new int [numberNonLinearRows_];
     1359      // do row numbers now
     1360      numberNonLinearRows_=0;
     1361      for (i=0;i<numberRows2;i++) {
     1362        if (whichRows[i]>=0) {
     1363          rowNonLinear_[numberNonLinearRows_++]=i;
     1364        }
     1365      }
     1366      whichNonLinear_ = new int [n];
     1367      for ( i =0;i<numberObjects_;i++) {
     1368        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
     1369        if (obj) {
     1370          int xyRow = obj->xyRow();
     1371          int k=whichRows[xyRow];
     1372          int put = pos[k];
     1373          pos[k]++;
     1374          whichNonLinear_[put]=i;
     1375        }
     1376      }
     1377      delete [] pos;
    12471378      delete [] whichRows;
    1248       delete [] whichColumns;
     1379      analyzeObjects();
    12491380    }
    12501381  }
     
    13121443  delete [] which;
    13131444}
     1445// Analyze constraints to see which are convex (quadratic)
     1446void
     1447OsiSolverLink::analyzeObjects()
     1448{
     1449  // space for starts
     1450  int numberColumns = coinModel_.numberColumns();
     1451  int * start = new int [numberColumns+1];
     1452  const double * rowLower = getRowLower();
     1453  const double * rowUpper = getRowUpper();
     1454  for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
     1455    int iRow = rowNonLinear_[iNon];
     1456    int numberElements = startNonLinear_[iNon+1]-startNonLinear_[iNon];
     1457    // triplet arrays
     1458    int * iColumn = new int [2*numberElements+1];
     1459    int * jColumn = new int [2*numberElements];
     1460    double * element = new double [2*numberElements];
     1461    int i;
     1462    int n=0;
     1463    for ( i =startNonLinear_[iNon];i<startNonLinear_[iNon+1];i++) {
     1464      OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[whichNonLinear_[i]]);
     1465      assert (obj);
     1466      int xColumn = obj->xColumn();
     1467      int yColumn = obj->yColumn();
     1468      double coefficient = obj->coefficient();
     1469      if (xColumn!=yColumn) {
     1470        iColumn[n]=xColumn;
     1471        jColumn[n]=yColumn;
     1472        element[n++]=coefficient;
     1473        iColumn[n]=yColumn;
     1474        jColumn[n]=xColumn;
     1475        element[n++]=coefficient;
     1476      } else {
     1477        iColumn[n]=xColumn;
     1478        jColumn[n]=xColumn;
     1479        element[n++]=coefficient;
     1480      }
     1481    }
     1482    // First sort in column order
     1483    CoinSort_3(iColumn,iColumn+n,jColumn,element);
     1484    // marker at end
     1485    iColumn[n]=numberColumns;
     1486    int lastI=iColumn[0];
     1487    // compute starts
     1488    start[0]=0;
     1489    for (i=1;i<n+1;i++) {
     1490      if (iColumn[i]!=lastI) {
     1491        while (lastI<iColumn[i]) {
     1492          start[lastI+1]=i;
     1493          lastI++;
     1494        }
     1495        lastI=iColumn[i];
     1496      }
     1497    }
     1498    // -1 unknown, 0 convex, 1 nonconvex
     1499    int status=-1;
     1500    int statusNegative=-1;
     1501    int numberLong=0; // number with >2 elements
     1502    for (int k=0;k<numberColumns;k++) {
     1503      int first = start[k];
     1504      int last = start[k+1];
     1505      if (last>first) {
     1506        int j;
     1507        double diagonal = 0.0;
     1508        int whichK=-1;
     1509        for (j=first;j<last;j++) {
     1510          if (jColumn[j]==k) {
     1511            diagonal = element[j];
     1512            status = diagonal >0 ? 0 : 1;
     1513            statusNegative = diagonal <0 ? 0 : 1;
     1514            whichK = (j==first) ? j+1 :j-1;
     1515            break;
     1516          }
     1517        }
     1518        if (last==first+1) {
     1519          // just one entry
     1520          if (!diagonal) {
     1521            // one off diagonal - not positive semi definite
     1522            status=1;
     1523            statusNegative=1;
     1524          }
     1525        } else if (diagonal) {
     1526          if (last==first+2) {
     1527            // other column and element
     1528            double otherElement=element[whichK];;
     1529            int otherColumn = jColumn[whichK];
     1530            double otherDiagonal=0.0;
     1531            // check 2x2 determinant - unless past and 2 long
     1532            if (otherColumn>i||start[otherColumn+1]>start[otherColumn]+2) {
     1533              for (j=start[otherColumn];j<start[otherColumn+1];j++) {
     1534                if (jColumn[j]==otherColumn) {
     1535                  otherDiagonal = element[j];
     1536                  break;
     1537                }
     1538              }
     1539              // determinant
     1540              double determinant = diagonal*otherDiagonal-otherElement*otherElement;
     1541              if (determinant<-1.0e-12) {
     1542                // not positive semi definite
     1543                status=1;
     1544                statusNegative=1;
     1545              } else if (start[otherColumn+1]>start[otherColumn]+2&&determinant<1.0e-12) {
     1546                // not positive semi definite
     1547                status=1;
     1548                statusNegative=1;
     1549              }
     1550            }
     1551          } else {
     1552            numberLong++;
     1553          }
     1554        }
     1555      }
     1556    }
     1557    if ((status==0||statusNegative==0)&&numberLong) {
     1558      // need to do more work
     1559      printf("Needs more work\n");
     1560    }
     1561    assert (status>0||statusNegative>0);
     1562    if (!status) {
     1563      convex_[iNon]=1;
     1564      // equality may be ok
     1565      assert (rowUpper[iRow]<1.0e20);
     1566      specialOptions2_ |= 8;
     1567    } else if (!statusNegative) {
     1568      convex_[iNon]=-1;
     1569      // equality may be ok
     1570      assert (rowLower[iRow]>-1.0e20);
     1571      specialOptions2_ |= 8;
     1572    } else {
     1573      convex_[iNon]=0;
     1574    }
     1575    printf("Convexity of row %d is %d\n",iRow,convex_[iNon]);
     1576    delete [] iColumn;
     1577    delete [] jColumn;
     1578    delete [] element;
     1579  }
     1580  delete [] start;
     1581}
    13141582//-------------------------------------------------------------------
    13151583// Clone
     
    13661634    delete [] bestSolution_;
    13671635    delete quadraticModel_;
     1636    delete [] startNonLinear_;
     1637    delete [] rowNonLinear_;
     1638    delete [] convex_;
     1639    delete [] whichNonLinear_;
    13681640  }
    13691641  matrix_ = NULL;
    13701642  originalRowCopy_ = NULL;
    13711643  quadraticModel_ = NULL;
     1644  numberNonLinearRows_=0;
     1645  startNonLinear_ = NULL;
     1646  rowNonLinear_ = NULL;
     1647  convex_ = NULL;
     1648  whichNonLinear_ = NULL;
    13721649  cbcModel_ = NULL;
    13731650  info_ = NULL;
     
    13881665  coinModel_ = rhs.coinModel_;
    13891666  numberVariables_ = rhs.numberVariables_;
     1667  numberNonLinearRows_ = rhs.numberNonLinearRows_;
    13901668  specialOptions2_ = rhs.specialOptions2_;
    13911669  objectiveRow_=rhs.objectiveRow_;
     
    14151693      bestSolution_=NULL;
    14161694    }
     1695  }
     1696  if (numberNonLinearRows_) {
     1697    startNonLinear_ = CoinCopyOfArray(rhs.startNonLinear_,numberNonLinearRows_+1);
     1698    rowNonLinear_ = CoinCopyOfArray(rhs.rowNonLinear_,numberNonLinearRows_);
     1699    convex_ = CoinCopyOfArray(rhs.convex_,numberNonLinearRows_);
     1700    int numberEntries = startNonLinear_[numberNonLinearRows_];
     1701    whichNonLinear_ = CoinCopyOfArray(rhs.whichNonLinear_,numberEntries);
    14171702  }
    14181703  if (rhs.quadraticModel_) {
  • branches/devel/Cbc/src/CbcLinked.hpp

    r500 r501  
    9191  /// Update coefficients - returns number updated if in updating mode
    9292  int updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix);
     93  /// Analyze constraints to see which are convex (quadratic)
     94  void analyzeObjects();
    9395  /// Objective value of best solution found internally
    9496  inline double bestObjectiveValue() const
     
    181183  /// Pointer back to CbcModel
    182184  CbcModel * cbcModel_;
     185  /// Number of rows with nonLinearities
     186  int numberNonLinearRows_;
     187  /// Starts of lists
     188  int * startNonLinear_;
     189  /// Row number for a list
     190  int * rowNonLinear_;
     191  /** Indicator whether is convex, concave or neither
     192      -1 concave, 0 neither, +1 convex
     193  */
     194  int * convex_;
     195  /// Indices in a list/row
     196  int * whichNonLinear_;
    183197  /// Model in CoinModel format
    184198  CoinModel coinModel_;
     
    189203  /**
    190204     0 bit (1) - don't do mini B&B
    191      1 bit (2) - quadratic only in objective
     205     1 bit (2) - quadratic only in objective (add OA cuts)
    192206     2 bit (4) - convex
     207     4 bit (8) - try adding OA cuts
    193208  */
    194209  int specialOptions2_;
Note: See TracChangeset for help on using the changeset viewer.