Changeset 678 for branches/devel/Cbc/src/CbcLinked.cpp
 Timestamp:
 Jul 9, 2007 4:52:34 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/devel/Cbc/src/CbcLinked.cpp
r667 r678 181 181 //sprintf(temp,"cc%d",iPass); 182 182 //writeMps(temp); 183 //writeMps("tight"); 184 //exit(33); 183 185 //printf("wrote cc%d\n",iPass); 184 186 OsiClpSolverInterface::initialSolve(); … … 1497 1499 } 1498 1500 } 1501 delete [] which; 1502 #if 1 1503 addTighterConstraints(); 1504 #endif 1505 } 1506 // Add reformulated bilinear constraints 1507 void 1508 OsiSolverLink::addTighterConstraints() 1509 { 1510 // This is first attempt  for now get working on trimloss 1511 int numberW=0; 1512 int * xW = new int[numberObjects_]; 1513 int * yW = new int[numberObjects_]; 1514 // Points to firstlambda 1515 int * wW = new int[numberObjects_]; 1516 // Coefficient 1517 double * alphaW = new double[numberObjects_]; 1518 // Objects 1519 OsiBiLinear ** objW = new OsiBiLinear * [numberObjects_]; 1520 int numberColumns = getNumCols(); 1521 int firstLambda=numberColumns; 1522 // set up list (better to rethink and do properly as column ordered) 1523 int * list = new int[numberColumns]; 1524 memset(list,0,numberColumns*sizeof(int)); 1525 int i; 1526 for ( i =0;i<numberObjects_;i++) { 1527 OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]); 1528 if (obj) { 1529 //obj>setBranchingStrategy(4); // ***** temp 1530 objW[numberW]=obj; 1531 xW[numberW]=obj>xColumn(); 1532 yW[numberW]=obj>yColumn(); 1533 list[xW[numberW]]=1; 1534 list[yW[numberW]]=1; 1535 wW[numberW]=obj>firstLambda(); 1536 firstLambda = CoinMin(firstLambda,obj>firstLambda()); 1537 alphaW[numberW]=obj>coefficient(); 1538 //assert (alphaW[numberW]==1.0); // fix when occurs 1539 numberW++; 1540 } 1541 } 1542 int nList = 0; 1543 for (i=0;i<numberColumns;i++) { 1544 if (list[i]) 1545 list[nList++]=i; 1546 } 1547 // set up mark array 1548 char * mark = new char [firstLambda*firstLambda]; 1549 memset(mark,0,firstLambda*firstLambda); 1550 for (i=0;i<numberW;i++) { 1551 int x = xW[i]; 1552 int y = yW[i]; 1553 mark[x*firstLambda+y]=1; 1554 mark[y*firstLambda+x]=1; 1555 } 1556 int numberRows2 = originalRowCopy_>getNumRows(); 1557 int * addColumn = new int [numberColumns]; 1558 double * addElement = new double [numberColumns]; 1559 assert (objectiveRow_<0); // fix when occurs 1560 for (int iRow=0;iRow<numberRows2;iRow++) { 1561 for (int iList=0;iList<nList;iList++) { 1562 int kColumn = list[iList]; 1563 const double * columnLower = getColLower(); 1564 //const double * columnUpper = getColUpper(); 1565 const double * rowLower = getRowLower(); 1566 const double * rowUpper = getRowUpper(); 1567 const CoinPackedMatrix * rowCopy = getMatrixByRow(); 1568 const double * element = rowCopy>getElements(); 1569 const int * column = rowCopy>getIndices(); 1570 const CoinBigIndex * rowStart = rowCopy>getVectorStarts(); 1571 const int * rowLength = rowCopy>getVectorLengths(); 1572 CoinBigIndex j; 1573 int numberElements = rowLength[iRow]; 1574 int n=0; 1575 for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) { 1576 int iColumn = column[j]; 1577 if (iColumn>=firstLambda) { 1578 // no good 1579 n=1; 1580 break; 1581 } 1582 if (mark[iColumn*firstLambda+kColumn]) 1583 n++; 1584 } 1585 if (n==numberElements) { 1586 printf("can add row %d\n",iRow); 1587 assert (columnLower[kColumn]>=0); // might be able to fix 1588 int xColumn=kColumn; 1589 n=0; 1590 for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) { 1591 int yColumn = column[j]; 1592 int k; 1593 for (k=0;k<numberW;k++) { 1594 if ((xW[k]==yColumn&&yW[k]==xColumn) 1595 (yW[k]==yColumn&&xW[k]==xColumn)) 1596 break; 1597 } 1598 assert (k<numberW); 1599 if (xW[k]!=xColumn) { 1600 int temp=xColumn; 1601 xColumn=yColumn; 1602 yColumn=temp; 1603 } 1604 int start = wW[k]; 1605 double value = element[j]; 1606 for (int kk=0;kk<4;kk++) { 1607 // Dummy value 1608 addElement[n]= 1.0e19; 1609 addColumn[n++]=start+kk; 1610 } 1611 // Tell object about this 1612 objW[k]>addExtraRow(matrix_>getNumRows(),value); 1613 } 1614 addColumn[n++] = kColumn; 1615 double lo = rowLower[iRow]; 1616 double up = rowUpper[iRow]; 1617 if (lo>1.0e20) { 1618 addElement[n1]=lo; 1619 if (lo==up) 1620 addRow(n,addColumn,addElement,0.0,0.0); 1621 else 1622 addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX); 1623 matrix_>appendRow(n,addColumn,addElement); 1624 } 1625 if (up<1.0e20&&up>lo) { 1626 addElement[n1]=up; 1627 addRow(n,addColumn,addElement,COIN_DBL_MAX,0.0); 1628 matrix_>appendRow(n,addColumn,addElement); 1629 } 1630 } 1631 } 1632 } 1499 1633 #if 0 1500 //int fake[10]={3,4,2,5,5,3,1,11,2,0}; 1501 int fake[10]={11,5,5,4,3,3,2,2,1,0}; 1502 for (int kk=0;kk<10;kk++) { 1503 setColUpper(kk,fake[kk]); 1504 setColLower(kk,fake[kk]); 1634 // possibly do bounds 1635 for (int iColumn=0;iColumn<firstLambda;iColumn++) { 1636 for (int iList=0;iList<nList;iList++) { 1637 int kColumn = list[iList]; 1638 const double * columnLower = getColLower(); 1639 const double * columnUpper = getColUpper(); 1640 if (mark[iColumn*firstLambda+kColumn]) { 1641 printf("can add column %d\n",iColumn); 1642 assert (columnLower[kColumn]>=0); // might be able to fix 1643 int xColumn=kColumn; 1644 int yColumn = iColumn; 1645 int k; 1646 for (k=0;k<numberW;k++) { 1647 if ((xW[k]==yColumn&&yW[k]==xColumn) 1648 (yW[k]==yColumn&&xW[k]==xColumn)) 1649 break; 1650 } 1651 assert (k<numberW); 1652 if (xW[k]!=xColumn) { 1653 int temp=xColumn; 1654 xColumn=yColumn; 1655 yColumn=temp; 1656 } 1657 int start = wW[k]; 1658 int n=0; 1659 for (int kk=0;kk<4;kk++) { 1660 // Dummy value 1661 addElement[n]= 1.0e19; 1662 addColumn[n++]=start+kk; 1663 } 1664 // Tell object about this 1665 objW[k]>addExtraRow(matrix_>getNumRows(),1.0); 1666 addColumn[n++] = kColumn; 1667 double lo = columnLower[iColumn]; 1668 double up = columnUpper[iColumn]; 1669 if (lo>1.0e20) { 1670 addElement[n1]=lo; 1671 if (lo==up) 1672 addRow(n,addColumn,addElement,0.0,0.0); 1673 else 1674 addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX); 1675 matrix_>appendRow(n,addColumn,addElement); 1676 } 1677 if (up<1.0e20&&up>lo) { 1678 addElement[n1]=up; 1679 addRow(n,addColumn,addElement,COIN_DBL_MAX,0.0); 1680 matrix_>appendRow(n,addColumn,addElement); 1681 } 1682 } 1683 } 1505 1684 } 1506 1685 #endif 1507 delete [] which; 1686 delete [] xW; 1687 delete [] yW; 1688 delete [] wW; 1689 delete [] alphaW; 1690 delete [] addColumn; 1691 delete [] addElement; 1692 delete [] mark; 1693 delete [] list; 1694 delete [] objW; 1508 1695 } 1509 1696 // Set all biLinear priorities on xx variables … … 4441 4628 xyRow_(1), 4442 4629 convexity_(1), 4630 numberExtraRows_(0), 4631 multiplier_(NULL), 4632 extraRow_(NULL), 4443 4633 chosen_(1) 4444 4634 { … … 4469 4659 xyRow_(xyRow), 4470 4660 convexity_(1), 4661 numberExtraRows_(0), 4662 multiplier_(NULL), 4663 extraRow_(NULL), 4471 4664 chosen_(1) 4472 4665 { … … 4680 4873 xyRow_(xyRow), 4681 4874 convexity_(1), 4875 numberExtraRows_(0), 4876 multiplier_(NULL), 4877 extraRow_(NULL), 4682 4878 chosen_(1) 4683 4879 { … … 4896 5092 xyRow_(rhs.xyRow_), 4897 5093 convexity_(rhs.convexity_), 5094 numberExtraRows_(rhs.numberExtraRows_), 5095 multiplier_(NULL), 5096 extraRow_(NULL), 4898 5097 chosen_(rhs.chosen_) 4899 5098 { 5099 if (numberExtraRows_) { 5100 multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_); 5101 extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_); 5102 } 4900 5103 } 4901 5104 … … 4931 5134 xyRow_ = rhs.xyRow_; 4932 5135 convexity_ = rhs.convexity_; 5136 numberExtraRows_ = rhs.numberExtraRows_; 5137 delete [] multiplier_; 5138 delete [] extraRow_; 5139 if (numberExtraRows_) { 5140 multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_); 5141 extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_); 5142 } else { 5143 multiplier_ = NULL; 5144 extraRow_ = NULL; 5145 } 4933 5146 chosen_ = rhs.chosen_; 4934 5147 } … … 4939 5152 OsiBiLinear::~OsiBiLinear () 4940 5153 { 5154 delete [] multiplier_; 5155 delete [] extraRow_; 5156 } 5157 // Adds in data for extra row with variable coefficients 5158 void 5159 OsiBiLinear::addExtraRow(int row, double multiplier) 5160 { 5161 int * tempI = new int [numberExtraRows_+1]; 5162 double * tempD = new double [numberExtraRows_+1]; 5163 memcpy(tempI,extraRow_,numberExtraRows_*sizeof(int)); 5164 memcpy(tempD,multiplier_,numberExtraRows_*sizeof(double)); 5165 tempI[numberExtraRows_]=row; 5166 tempD[numberExtraRows_]=multiplier; 5167 numberExtraRows_++; 5168 delete [] extraRow_; 5169 extraRow_ = tempI; 5170 delete [] multiplier_; 5171 multiplier_ = tempD; 4941 5172 } 4942 5173 static bool testCoarse=true; … … 5121 5352 } 5122 5353 xyLambda /= coefficient_; 5354 } 5355 if (0) { 5356 // only true with positive values 5357 // see if all convexification constraints OK with true 5358 assert (xyTrue+1.0e5>xB[0]*y+yB[0]*x  xB[0]*yB[0]); 5359 assert (xyTrue+1.0e5>xB[1]*y+yB[1]*x  xB[1]*yB[1]); 5360 assert (xyTrue1.0e5<xB[1]*y+yB[0]*x  xB[1]*yB[0]); 5361 assert (xyTrue1.0e5<xB[0]*y+yB[1]*x  xB[0]*yB[1]); 5362 // see if all convexification constraints OK with lambda version 5363 #if 1 5364 assert (xyLambda+1.0e5>xB[0]*y+yB[0]*x  xB[0]*yB[0]); 5365 assert (xyLambda+1.0e5>xB[1]*y+yB[1]*x  xB[1]*yB[1]); 5366 assert (xyLambda1.0e5<xB[1]*y+yB[0]*x  xB[1]*yB[0]); 5367 assert (xyLambda1.0e5<xB[0]*y+yB[1]*x  xB[0]*yB[1]); 5368 #endif 5369 // see if other bound stuff true 5370 assert (xyLambda+1.0e5>xB[0]*y); 5371 assert (xyLambda+1.0e5>yB[0]*x); 5372 assert (xyLambda1.0e5<xB[1]*y); 5373 assert (xyLambda1.0e5<yB[1]*x); 5374 #define SIZE 2 5375 if (yColumn_==xColumn_+SIZE) { 5376 #if SIZE==6 5377 double bMax = 2200.0; 5378 double bMin = bMax  100.0; 5379 double b[] = {330.0,360.0,380.0,430.0,490.0,530.0}; 5380 #elif SIZE==2 5381 double bMax = 1900.0; 5382 double bMin = bMax  200.0; 5383 double b[] = {460.0,570.0}; 5384 #else 5385 abort(); 5386 #endif 5387 double sum =0.0; 5388 double sum2 =0.0; 5389 int m=xColumn_; 5390 double x = info>solution_[m]; 5391 double xB[2]; 5392 double yB[2]; 5393 xB[0]=info>lower_[m]; 5394 xB[1]=info>upper_[m]; 5395 for (int i=0;i<SIZE*SIZE;i+=SIZE) { 5396 int n = i+SIZE+m; 5397 double y = info>solution_[n]; 5398 yB[0]=info>lower_[n]; 5399 yB[1]=info>upper_[n]; 5400 int firstLambda=SIZE*SIZE+2*SIZE+4*i+4*m; 5401 double xyLambda=0.0; 5402 for (int j=0;j<4;j++) { 5403 int iX = j>>1; 5404 int iY = j&1; 5405 xyLambda += xB[iX]*yB[iY]*info>solution_[firstLambda+j]; 5406 } 5407 sum += xyLambda*b[i/SIZE]; 5408 double xyTrue = x*y; 5409 sum2 += xyTrue*b[i/SIZE]; 5410 } 5411 if (sum>bMax*x+1.0e5sum<bMin*x1.0e5) { 5412 //if (sum<bMax*x+1.0e5&&sum>bMin*x1.0e5) { 5413 printf("bmin*x %g b*w %g bmax*x %g (true) %g\n",bMin*x,sum,bMax*x,sum2); 5414 printf("m %d lb %g value %g up %g\n", 5415 m,xB[0],x,xB[1]); 5416 sum=0.0; 5417 for (int i=0;i<SIZE*SIZE;i+=SIZE) { 5418 int n = i+SIZE+m; 5419 double y = info>solution_[n]; 5420 yB[0]=info>lower_[n]; 5421 yB[1]=info>upper_[n]; 5422 printf("n %d lb %g value %g up %g\n", 5423 n,yB[0],y,yB[1]); 5424 int firstLambda=SIZE*SIZE+2*SIZE+4*i+m*4; 5425 double xyLambda=0.0; 5426 for (int j=0;j<4;j++) { 5427 int iX = j>>1; 5428 int iY = j&1; 5429 xyLambda += xB[iX]*yB[iY]*info>solution_[firstLambda+j]; 5430 printf("j %d l %d new xylambda %g ",j,firstLambda+j,xyLambda); 5431 } 5432 sum += xyLambda*b[i/SIZE]; 5433 printf("  sum now %g\n",sum); 5434 } 5435 } 5436 if (sum2>bMax*x+1.0e5sum2<bMin*x1.0e5) { 5437 printf("bmin*x %g b*x*y %g bmax*x %g (estimate) %g\n",bMin*x,sum2,bMax*x,sum); 5438 printf("m %d lb %g value %g up %g\n", 5439 m,xB[0],x,xB[1]); 5440 sum2=0.0; 5441 for (int i=0;i<SIZE*SIZE;i+=SIZE) { 5442 int n = i+SIZE+m; 5443 double y = info>solution_[n]; 5444 yB[0]=info>lower_[n]; 5445 yB[1]=info>upper_[n]; 5446 printf("n %d lb %g value %g up %g\n", 5447 n,yB[0],y,yB[1]); 5448 double xyTrue = x*y; 5449 sum2 += xyTrue*b[i/SIZE]; 5450 printf("xyTrue %g  sum now %g\n",xyTrue,sum2); 5451 } 5452 } 5453 } 5123 5454 } 5124 5455 // If pseudo shadow prices then see what would happen … … 5841 6172 const int * row = matrix>getIndices(); 5842 6173 const CoinBigIndex * columnStart = matrix>getVectorStarts(); 5843 //const int * columnLength = matrix>getVectorLengths();6174 const int * columnLength = matrix>getVectorLengths(); 5844 6175 // order is LxLy, LxUy, UxLy and UxUy 5845 6176 double xB[2]; … … 5862 6193 double y = yB[iY]; 5863 6194 CoinBigIndex k = columnStart[j+firstLambda_]; 6195 CoinBigIndex last = k+columnLength[j+firstLambda_]; 5864 6196 double value; 5865 6197 // xy … … 5896 6228 element[k++]=value; 5897 6229 numberUpdated++; 6230 } 6231 // Do extra rows 6232 for (int i=0;i<numberExtraRows_;i++) { 6233 int iRow = extraRow_[i]; 6234 for (;k<last;k++) { 6235 if (row[k]==iRow) 6236 break; 6237 } 6238 assert (k<last); 6239 element[k++] = x*y*multiplier_[i]; 5898 6240 } 5899 6241 }
Note: See TracChangeset
for help on using the changeset viewer.