Changeset 501 for branches/devel/Cbc
 Timestamp:
 Dec 12, 2006 4:22:19 PM (14 years ago)
 Location:
 branches/devel/Cbc/src
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

branches/devel/Cbc/src/CbcLinked.cpp
r500 r501 557 557 } 558 558 // ???  try 559 // But skip if strong branching 560 CbcModel * cbcModel = (modelPtr_>maximumIterations()<10000) ? cbcModel_ : NULL; 559 561 if ((specialOptions2_&2)!=0) { 560 562 // If model has stored then add cut (if convex) 561 if (cbcModel _&&(specialOptions2_&4)!=0&&quadraticModel_) {563 if (cbcModel&&(specialOptions2_&4)!=0&&quadraticModel_) { 562 564 int numberGenerators = cbcModel_>numberCutGenerators(); 563 565 int iGenerator; … … 627 629 } 628 630 } 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.0e12) { 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.0e5) 702 gen2>addCut(COIN_DBL_MAX,offset+1.0e7,n,column,gradient); 703 else if (!convex&&rhs<offset1.0e5) 704 gen2>addCut(offset1.0e7,COIN_DBL_MAX,n,column,gradient); 705 } 706 delete [] gradient; 707 delete [] column; 708 break; 709 } 710 } 629 711 } 630 712 } … … 826 908 int * which = new int[numberColumns]; 827 909 numberVariables_=0; 828 specialOptions2_=0;910 //specialOptions2_=0; 829 911 int iColumn; 830 912 int numberErrors=0; … … 1024 1106 elementQuadratic = new double [numberQuadratic]; 1025 1107 numberQuadratic=0; 1108 } 1109 if (quadraticObjective((specialOptions2_&8)!=0&&saveNBi)) { 1026 1110 // add in objective as constraint 1027 1111 objectiveVariable_= numberColumns; … … 1245 1329 numberRows2,whichRows, 1246 1330 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; 1247 1378 delete [] whichRows; 1248 delete [] whichColumns;1379 analyzeObjects(); 1249 1380 } 1250 1381 } … … 1312 1443 delete [] which; 1313 1444 } 1445 // Analyze constraints to see which are convex (quadratic) 1446 void 1447 OsiSolverLink::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 :j1; 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>istart[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*otherDiagonalotherElement*otherElement; 1541 if (determinant<1.0e12) { 1542 // not positive semi definite 1543 status=1; 1544 statusNegative=1; 1545 } else if (start[otherColumn+1]>start[otherColumn]+2&&determinant<1.0e12) { 1546 // not positive semi definite 1547 status=1; 1548 statusNegative=1; 1549 } 1550 } 1551 } else { 1552 numberLong++; 1553 } 1554 } 1555 } 1556 } 1557 if ((status==0statusNegative==0)&&numberLong) { 1558 // need to do more work 1559 printf("Needs more work\n"); 1560 } 1561 assert (status>0statusNegative>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 } 1314 1582 // 1315 1583 // Clone … … 1366 1634 delete [] bestSolution_; 1367 1635 delete quadraticModel_; 1636 delete [] startNonLinear_; 1637 delete [] rowNonLinear_; 1638 delete [] convex_; 1639 delete [] whichNonLinear_; 1368 1640 } 1369 1641 matrix_ = NULL; 1370 1642 originalRowCopy_ = NULL; 1371 1643 quadraticModel_ = NULL; 1644 numberNonLinearRows_=0; 1645 startNonLinear_ = NULL; 1646 rowNonLinear_ = NULL; 1647 convex_ = NULL; 1648 whichNonLinear_ = NULL; 1372 1649 cbcModel_ = NULL; 1373 1650 info_ = NULL; … … 1388 1665 coinModel_ = rhs.coinModel_; 1389 1666 numberVariables_ = rhs.numberVariables_; 1667 numberNonLinearRows_ = rhs.numberNonLinearRows_; 1390 1668 specialOptions2_ = rhs.specialOptions2_; 1391 1669 objectiveRow_=rhs.objectiveRow_; … … 1415 1693 bestSolution_=NULL; 1416 1694 } 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); 1417 1702 } 1418 1703 if (rhs.quadraticModel_) { 
branches/devel/Cbc/src/CbcLinked.hpp
r500 r501 91 91 /// Update coefficients  returns number updated if in updating mode 92 92 int updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix); 93 /// Analyze constraints to see which are convex (quadratic) 94 void analyzeObjects(); 93 95 /// Objective value of best solution found internally 94 96 inline double bestObjectiveValue() const … … 181 183 /// Pointer back to CbcModel 182 184 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_; 183 197 /// Model in CoinModel format 184 198 CoinModel coinModel_; … … 189 203 /** 190 204 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) 192 206 2 bit (4)  convex 207 4 bit (8)  try adding OA cuts 193 208 */ 194 209 int specialOptions2_;
Note: See TracChangeset
for help on using the changeset viewer.