#include <CglGomory.hpp>
Inheritance diagram for CglGomory:
Public Member Functions | |
Generate Cuts | |
virtual void | generateCuts (const OsiSolverInterface &si, OsiCuts &cs) const |
int | generateCuts (const OsiRowCutDebugger *debugger, OsiCuts &cs, const CoinPackedMatrix &columnCopy, const double *objective, const double *colsol, const double *colLower, const double *colUpper, const double *rowLower, const double *rowUpper, const char *intVar, const CoinWarmStartBasis *warm) const |
Change limit on how many variables in cut (default 50) | |
void | setLimit (int limit) |
Set. | |
int | getLimit () const |
Get. | |
Change criterion on which variables to look at. All ones | |
more than "away" away from integrality will be investigated (default 0.05) | |
void | setAway (double value) |
Set. | |
double | getAway () const |
Get. | |
Constructors and destructors | |
CglGomory () | |
Default constructor. | |
CglGomory (const CglGomory &) | |
Copy constructor. | |
CglGomory & | operator= (const CglGomory &rhs) |
Assignment operator. | |
virtual | ~CglGomory () |
Destructor. | |
Private Attributes | |
Private member data | |
double | away_ |
Only investigate if more than this away from integrality. | |
int | limit_ |
Limit - only generate if fewer than this in cut. | |
Friends | |
void | CglGomoryUnitTest (const OsiSolverInterface *siP, const std::string mpdDir) |
Definition at line 12 of file CglGomory.hpp.
|
Generates cuts given matrix and solution etc, returns number of cuts generated Definition at line 198 of file CglGomory.cpp. References CoinIndexedVector::add(), away_, CoinIndexedVector::checkClear(), CoinIndexedVector::clear(), CoinIndexedVector::denseVector(), CoinWarmStartBasis::getArtifStatus(), CoinIndexedVector::getIndices(), CoinIndexedVector::getNumElements(), CoinWarmStartBasis::getStructStatus(), limit_, CoinIndexedVector::reserve(), CoinIndexedVector::setNumElements(), and CoinIndexedVector::setVector().
00206 { 00207 int numberRows=columnCopy.getNumRows(); 00208 int numberColumns=columnCopy.getNumCols(); 00209 00210 // Start of code to create a factorization from warm start (A) ==== 00211 // check factorization is okay 00212 CoinFactorization factorization; 00213 // We can either set increasing rows so ...IsBasic gives pivot row 00214 // or we can just increment iBasic one by one 00215 // for now let ...iBasic give pivot row 00216 factorization.increasingRows(2); 00217 int status=-100; 00218 // probably could use pivotVariables from OsiSimplexModel 00219 int * rowIsBasic = new int[numberRows]; 00220 int * columnIsBasic = new int[numberColumns]; 00221 int i; 00222 int numberBasic=0; 00223 for (i=0;i<numberRows;i++) { 00224 if (warm->getArtifStatus(i) == CoinWarmStartBasis::basic) { 00225 rowIsBasic[i]=1; 00226 numberBasic++; 00227 } else { 00228 rowIsBasic[i]=-1; 00229 } 00230 } 00231 for (i=0;i<numberColumns;i++) { 00232 if (warm->getStructStatus(i) == CoinWarmStartBasis::basic) { 00233 columnIsBasic[i]=1; 00234 numberBasic++; 00235 } else { 00236 columnIsBasic[i]=-1; 00237 } 00238 } 00239 //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */ 00240 while (status<-98) { 00241 status=factorization.factorize(columnCopy, 00242 rowIsBasic, columnIsBasic); 00243 if (status==-99) factorization.areaFactor(factorization.areaFactor() * 2.0); 00244 } 00245 if (status) { 00246 std::cout<<"Bad factorization of basis - status "<<status<<std::endl; 00247 return -1; 00248 } 00249 // End of creation of factorization (A) ==== 00250 00251 double bounds[2]={-DBL_MAX,0.0}; 00252 int iColumn,iRow; 00253 00254 // get row copy of matrix 00255 CoinPackedMatrix rowCopy = columnCopy; 00256 rowCopy.reverseOrdering(); 00257 const int * column = rowCopy.getIndices(); 00258 const int * rowStart = rowCopy.getVectorStarts(); 00259 const int * rowLength = rowCopy.getVectorLengths(); 00260 const double * rowElements = rowCopy.getElements(); 00261 const int * row = columnCopy.getIndices(); 00262 const int * columnStart = columnCopy.getVectorStarts(); 00263 const int * columnLength = columnCopy.getVectorLengths(); 00264 const double * columnElements = columnCopy.getElements(); 00265 00266 // we need to do book-keeping for variables at ub 00267 double tolerance = 1.0e-7; 00268 bool * swap= new bool [numberColumns]; 00269 for (iColumn=0;iColumn<numberColumns;iColumn++) { 00270 if (columnIsBasic[iColumn]<0&& 00271 colUpper[iColumn]-colsol[iColumn]<=tolerance) { 00272 swap[iColumn]=true; 00273 } else { 00274 swap[iColumn]=false; 00275 } 00276 } 00277 00278 // get row activities (could use solver but lets do here ) 00279 double * rowActivity = new double [numberRows]; 00280 CoinFillN(rowActivity,numberRows,0.0); 00281 for (iColumn=0;iColumn<numberColumns;iColumn++) { 00282 double value = colsol[iColumn]; 00283 int k; 00284 for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];k++) { 00285 iRow = row[k]; 00286 rowActivity[iRow] += columnElements[k]*value; 00287 } 00288 } 00289 /* we need to mark rows: 00290 0) equality 00291 1) slack at lb (activity at ub) 00292 2) slack at ub (activity at lb) 00293 and 4 bit is set if slack must be integer 00294 00295 */ 00296 int * rowType = new int[numberRows]; 00297 for (iRow=0;iRow<numberRows;iRow++) { 00298 if (rowIsBasic[iRow]<0&&rowUpper[iRow]>rowLower[iRow]+1.0e-7) { 00299 int type=0; 00300 double rhs=0.0; 00301 if (rowActivity[iRow]>=rowUpper[iRow]-1.0e-7) { 00302 type=1; 00303 rhs=rowUpper[iRow]; 00304 } else if (rowActivity[iRow]<=rowLower[iRow]+1.0e-7) { 00305 type=2; 00306 rhs=rowLower[iRow]; 00307 } else { 00308 // wrong - but probably large rhs 00309 rowType[iRow]=type; 00310 #ifdef CGL_DEBUG 00311 assert (min(rowUpper[iRow]-rowActivity[iRow], 00312 rowActivity[iRow]-rowUpper[iRow])<1.0e-7); 00313 abort(); 00314 #else 00315 continue; 00316 #endif 00317 } 00318 if (above_integer(rhs)<1.0e-10) { 00319 // could be integer slack 00320 bool allInteger=true; 00321 int k; 00322 for (k=rowStart[iRow]; 00323 k<rowStart[iRow]+rowLength[iRow];k++) { 00324 int iColumn=column[k]; 00325 if (!intVar[iColumn]||above_integer(rowElements[k])>1.0e-10) { 00326 // not integer slacks 00327 allInteger=false; 00328 break; 00329 } 00330 } 00331 if (allInteger) { 00332 type |= 4; 00333 } 00334 } 00335 rowType[iRow]=type; 00336 } else { 00337 // row is equality or basic 00338 rowType[iRow]=0; 00339 } 00340 } 00341 00342 // Start of code to create work arrays for factorization (B) ==== 00343 // two vectors for updating (one is work) 00344 CoinIndexedVector work; 00345 CoinIndexedVector array; 00346 // make sure large enough 00347 work.reserve(numberRows); 00348 array.reserve(numberRows); 00349 int * arrayRows = array.getIndices(); 00350 double * arrayElements = array.denseVector(); 00351 // End of code to create work arrays (B) ==== 00352 00353 int numberAdded=0; 00354 // we also need somewhere to accumulate cut 00355 CoinIndexedVector cutVector; 00356 cutVector.reserve(numberColumns+1); 00357 int * cutIndex = cutVector.getIndices(); 00358 double * cutElement = cutVector.denseVector(); 00359 // and for packed form (as not necessarily in order) 00360 double * packed = new double[numberColumns+1]; 00361 00362 for (iColumn=0;iColumn<numberColumns;iColumn++) { 00363 double reducedValue=above_integer(colsol[iColumn]);; 00364 // This returns pivot row for columns or -1 if not basic (C) ==== 00365 int iBasic=columnIsBasic[iColumn]; 00366 double ratio=reducedValue/(1.0-reducedValue); 00367 if (iBasic>=0) { 00368 // Debug code below computes tableau column of basic ==== 00369 int j; 00370 #ifdef CGL_DEBUG 00371 { 00372 // put column into array 00373 array.setVector(columnLength[iColumn],row+columnStart[iColumn], 00374 columnElements+columnStart[iColumn]); 00375 // get column in tableau 00376 factorization.updateColumn ( &work, &array ); 00377 int nn=0; 00378 int numberInArray=array.getNumElements(); 00379 for (j=0;j<numberInArray;j++) { 00380 int indexValue=arrayRows[j]; 00381 double value=arrayElements[indexValue]; 00382 if (fabs(value)>1.0e-6) { 00383 assert (fabs(value-1.0)<1.0e-7); 00384 assert (indexValue==iBasic); 00385 nn++; 00386 } 00387 } 00388 assert (nn==1); 00389 } 00390 #endif 00391 if(intVar[iColumn]&&reducedValue<1.0-away_&&reducedValue>away_) { 00392 #ifdef CGL_DEBUG 00393 cutVector.checkClear(); 00394 #endif 00395 // get row of tableau 00396 double one =1.0; 00397 array.setVector(1,&iBasic,&one); 00398 int numberNonInteger=0; 00399 //Code below computes tableau row ==== 00400 // get pi 00401 factorization.updateColumnTranspose ( &work, &array ); 00402 int numberInArray=array.getNumElements(); 00403 // check pivot on iColumn 00404 { 00405 double value=0.0; 00406 int k; 00407 // add in row of tableau 00408 for (k=columnStart[iColumn]; 00409 k<columnStart[iColumn]+columnLength[iColumn];k++) { 00410 iRow = row[k]; 00411 value += columnElements[k]*arrayElements[iRow]; 00412 } 00413 // should be 1 00414 #ifdef CGL_DEBUG 00415 assert (fabs(value-1.0) < 1.0e-7); 00416 #endif 00417 } 00418 //reducedValue=colsol[iColumn]; 00419 // coding from pg 130 of Wolsey 00420 // adjustment to rhs 00421 double rhs=0.0; 00422 int number=0; 00423 // columns 00424 for (j=0;j<numberColumns;j++) { 00425 double value=0.0; 00426 if (columnIsBasic[j]<0&& 00427 colUpper[j]>colLower[j]+1.0e-8) { 00428 int k; 00429 // add in row of tableau 00430 for (k=columnStart[j];k<columnStart[j]+columnLength[j];k++) { 00431 iRow = row[k]; 00432 value += columnElements[k]*arrayElements[iRow]; 00433 } 00434 // value is entry in tableau row end (C) ==== 00435 if (fabs(value)<1.0e-16) { 00436 // small value 00437 continue; 00438 } else { 00439 #if CGL_DEBUG>1 00440 if (iColumn==314) printf("for basic %d, column %d has alpha %g, colsol %g\n", 00441 iColumn,j,value,colsol[j]); 00442 #endif 00443 // deal with bounds 00444 if (swap[j]) { 00445 //reducedValue -= value*colUpper[j]; 00446 // negate 00447 value = - value; 00448 } else { 00449 //reducedValue -= value*colLower[j]; 00450 } 00451 #if CGL_DEBUG>1 00452 if (iColumn==348) printf("%d value %g reduced %g int %d rhs %g swap %d\n", 00453 j,value,reducedValue,intVar[j],rhs,swap[j]); 00454 #endif 00455 double coefficient; 00456 if (intVar[j]) { 00457 // integer 00458 coefficient = above_integer(value); 00459 if (coefficient > reducedValue) { 00460 coefficient = ratio * (1.0-coefficient); 00461 } 00462 } else { 00463 // continuous 00464 numberNonInteger++; 00465 if (value > 0.0) { 00466 coefficient = value; 00467 } else { 00468 //??? sign wrong in book 00469 coefficient = -ratio*value; 00470 } 00471 } 00472 if (swap[j]) { 00473 // negate 00474 coefficient = - coefficient; 00475 rhs += colUpper[j]*coefficient; 00476 } else { 00477 rhs += colLower[j]*coefficient; 00478 } 00479 if (fabs(coefficient)>= COIN_INDEXED_TINY_ELEMENT) { 00480 cutElement[j] = coefficient; 00481 cutIndex[number++]=j; 00482 // If too many - break from loop 00483 if (number>limit_) 00484 break; 00485 } 00486 } 00487 } else { 00488 // basic 00489 continue; 00490 } 00491 } 00492 cutVector.setNumElements(number); 00493 // If too many - just clear vector and skip 00494 if (number>limit_) { 00495 cutVector.clear(); 00496 continue; 00497 } 00498 //check will be cut 00499 //reducedValue=above_integer(reducedValue); 00500 rhs += reducedValue; 00501 double violation = reducedValue; 00502 #ifdef CGL_DEBUG 00503 std::cout<<"cut has violation of "<<violation 00504 <<" value "<<colsol[iColumn]<<std::endl; 00505 #endif 00506 // now do slacks part 00507 for (j=0;j<numberInArray;j++) { 00508 iRow=arrayRows[j]; 00509 double value = arrayElements[iRow]; 00510 int type=rowType[iRow]; 00511 if (type&&fabs(value)>=1.0e-16) { 00512 if ((type&1)==0) { 00513 // negate to get correct coefficient 00514 value = - value; 00515 } 00516 double coefficient; 00517 if ((type&4)!=0) { 00518 // integer 00519 coefficient = above_integer(value); 00520 if (coefficient > reducedValue) { 00521 coefficient = ratio * (1.0-coefficient); 00522 } 00523 } else { 00524 numberNonInteger++; 00525 // continuous 00526 if (value > 0.0) { 00527 coefficient = value; 00528 } else { 00529 coefficient = -ratio*value; 00530 } 00531 } 00532 if ((type&1)!=0) { 00533 // slack at ub - treat as +1.0 00534 rhs -= coefficient*rowUpper[iRow]; 00535 } else { 00536 // negate yet again ? 00537 coefficient = - coefficient; 00538 rhs -= coefficient*rowLower[iRow]; 00539 } 00540 int k; 00541 for (k=rowStart[iRow]; 00542 k<rowStart[iRow]+rowLength[iRow];k++) { 00543 int jColumn=column[k]; 00544 double value=rowElements[k]; 00545 cutVector.add(jColumn,-coefficient*value); 00546 } 00547 } 00548 } 00549 //check again and pack down 00550 // also change signs 00551 // also zero cutElement 00552 double sum=0.0; 00553 rhs = - rhs; 00554 int n = cutVector.getNumElements(); 00555 number=0; 00556 for (j=0;j<n;j++) { 00557 int jColumn =cutIndex[j]; 00558 double value=-cutElement[jColumn]; 00559 cutElement[jColumn]=0.0; 00560 if (fabs(value)>1.0e-13) { 00561 sum+=value*colsol[jColumn]; 00562 packed[number]=value; 00563 cutIndex[number++]=jColumn; 00564 } 00565 } 00566 // say zeroed out 00567 cutVector.setNumElements(0); 00568 if (sum >rhs+0.9*away_&& 00569 fabs((sum-rhs)-violation)<1.0e-6) { 00570 //#ifdef CGL_DEBUG 00571 #ifdef CGL_DEBUG 00572 #if CGL_DEBUG<=1 00573 if (number<=-10) { 00574 #endif 00575 for (j=0;j<number;j++) { 00576 std::cout<<" ["<<cutIndex[j]<<","<<packed[j]<<"]"; 00577 } 00578 std::cout<<" <= "<<rhs<<std::endl; 00579 #if CGL_DEBUG<=1 00580 } 00581 #endif 00582 #endif 00583 if (!numberNonInteger&&number) { 00584 #ifdef CGL_DEBUG 00585 assert (sizeof(Rational)==sizeof(double)); 00586 #endif 00587 Rational * cleaned = (Rational *) cutElement; 00588 int * xInt = (int *) cutElement; 00589 // cut should have an integer slack so try and simplify 00590 // add in rhs and put in cutElements (remember to zero later) 00591 cutIndex[number]=numberColumns+1; 00592 packed[number]=rhs; 00593 int numberNonSmall=0; 00594 int lcm = 1; 00595 00596 for (j=0;j<number+1;j++) { 00597 double value=above_integer(fabs(packed[j])); 00598 if (fabs(value)<1.0e-4) { 00599 // too small 00600 continue; 00601 } else { 00602 numberNonSmall++; 00603 } 00604 00605 cleaned[j]=nearestRational(value,100000); 00606 if (cleaned[j].denominator<0) { 00607 // bad rational 00608 lcm=-1; 00609 break; 00610 } 00611 int thisGcd = gcd(lcm,cleaned[j].denominator); 00612 // may need to check for overflow 00613 lcm /= thisGcd; 00614 lcm *= cleaned[j].denominator; 00615 } 00616 if (lcm>0&&numberNonSmall) { 00617 double multiplier = lcm; 00618 int nOverflow = 0; 00619 for (j=0; j<number+1;j++) { 00620 double value = fabs(packed[j]); 00621 double dxInt = value*multiplier; 00622 xInt[j]= (int) (dxInt+0.5); 00623 #if CGL_DEBUG>1 00624 printf("%g => %g \n",value,dxInt); 00625 #endif 00626 if (dxInt>1.0e9||fabs(dxInt-xInt[j])> 1.0e-8) { 00627 nOverflow++; 00628 break; 00629 } 00630 } 00631 00632 if (nOverflow){ 00633 #ifdef CGL_DEBUG 00634 printf("Gomory Scaling: Warning: Overflow detected \n"); 00635 #endif 00636 numberNonInteger=-1; 00637 } else { 00638 00639 // find greatest common divisor of the elements 00640 j=0; 00641 while (!xInt[j]) 00642 j++; // skip zeros 00643 int thisGcd = gcd(xInt[j],xInt[j+1]); 00644 j++; 00645 for (;j<number+1;j++) { 00646 if (xInt[j]) 00647 thisGcd = gcd(thisGcd,xInt[j]); 00648 } 00649 #if 0 00650 // Check nothing too illegal - FIX this 00651 for (j=0;j<number+1;j++) { 00652 double old = lcm*packed[j]; 00653 int newOne; 00654 if (old>0.0) 00655 newOne=xInt[j]/thisGcd; 00656 else 00657 newOne=-xInt[j]/thisGcd; 00658 if (fabs(((double) newOne)-old)> 00659 1.0e-10*(fabs(newOne)+1.0)) { 00660 // say no good - first see if happens 00661 printf("Fix this test 456 - just skip\n"); 00662 abort(); 00663 } 00664 } 00665 #endif 00666 #if CGL_DEBUG>1 00667 printf("The gcd of xInt is %i\n",thisGcd); 00668 #endif 00669 00670 // construct new cut by dividing through by gcd and 00671 double minMultiplier=1.0e100; 00672 double maxMultiplier=0.0; 00673 for (j=0; j<number+1;j++) { 00674 double old=packed[j]; 00675 if (old>0.0) { 00676 packed[j]=xInt[j]/thisGcd; 00677 } else { 00678 packed[j]=-xInt[j]/thisGcd; 00679 } 00680 #if CGL_DEBUG>1 00681 printf("%g => %g \n",old,packed[j]); 00682 #endif 00683 if (packed[j]) { 00684 if (fabs(packed[j])>maxMultiplier*fabs(old)) 00685 maxMultiplier = packed[j]/old; 00686 if (fabs(packed[j])<minMultiplier*fabs(old)) 00687 minMultiplier = packed[j]/old; 00688 } 00689 } 00690 rhs = packed[number]; 00691 double ratio=maxMultiplier/minMultiplier; 00692 #ifdef CGL_DEBUG 00693 printf("min, max multipliers - %g, %g\n", 00694 minMultiplier,maxMultiplier); 00695 #endif 00696 assert(ratio>0.9999&&ratio<1.0001); 00697 } 00698 } 00699 // erase cutElement 00700 CoinFillN(cutElement,number+1,0.0); 00701 } else { 00702 // relax rhs a tiny bit 00703 rhs += 1.0e-8; 00704 // relax if lots of elements for mixed gomory 00705 if (number>=20) { 00706 rhs += 1.0e-7*((double) (number/20)); 00707 } 00708 } 00709 // Take off tiny elements 00710 // for first pass reject 00711 #define TINY_ELEMENT 1.0e-7 00712 { 00713 int i,number2=number; 00714 number=0; 00715 double largest=0.0; 00716 double smallest=1.0e30; 00717 for (i=0;i<number2;i++) { 00718 double value=fabs(packed[i]); 00719 if (value<TINY_ELEMENT) { 00720 int iColumn = cutIndex[i]; 00721 if (colUpper[iColumn]-colLower[iColumn]<10.0) { 00722 // weaken cut 00723 if (packed[i]>0.0) 00724 rhs -= value*colLower[iColumn]; 00725 else 00726 rhs += value*colUpper[iColumn]; 00727 } else { 00728 // throw away 00729 number=limit_+1; 00730 numberNonInteger=1; 00731 break; 00732 } 00733 } else { 00734 largest=max(largest,value); 00735 smallest=min(smallest,value); 00736 cutIndex[number]=cutIndex[i]; 00737 packed[number++]=packed[i]; 00738 } 00739 } 00740 if (largest>1.0e9*smallest) { 00741 number=limit_+1; //reject 00742 numberNonInteger=1; 00743 } 00744 } 00745 if (number<limit_||!numberNonInteger) { 00746 bounds[1]=rhs; 00747 if (number>50&&numberNonInteger) 00748 bounds[1] += 1.0e-6; // weaken 00749 { 00750 OsiRowCut rc; 00751 rc.setRow(number,cutIndex,packed); 00752 rc.setLb(bounds[0]); 00753 rc.setUb(bounds[1]); 00754 cs.insert(rc); 00755 #ifdef CGL_DEBUG 00756 if (!number) 00757 std::cout<<"******* Empty cut - infeasible"<<std::endl; 00758 if (debugger) 00759 assert(!debugger->invalidCut(rc)); 00760 #endif 00761 } 00762 numberAdded++; 00763 } else { 00764 #ifdef CGL_DEBUG 00765 std::cout<<"cut has "<<number<<" entries - skipped" 00766 <<std::endl; 00767 if (!number) 00768 std::cout<<"******* Empty cut - infeasible"<<std::endl; 00769 #endif 00770 } 00771 } else { 00772 // why dropped? 00773 #ifdef CGL_DEBUG 00774 std::cout<<"first violation "<<violation<<" now " 00775 <<sum-rhs<<" why?, rhs= "<<rhs<<std::endl; 00776 00777 for (j=0;j<number;j++) { 00778 int jColumn =cutIndex[j]; 00779 double value=packed[j]; 00780 std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn] 00781 <<") "; 00782 } 00783 std::cout<<std::endl; 00784 abort(); 00785 #endif 00786 } 00787 } 00788 } else { 00789 // not basic 00790 #if CGL_DEBUG>1 00791 { 00792 // put column into array 00793 array.setVector(columnLength[iColumn],row+columnStart[iColumn], 00794 columnElements+columnStart[iColumn]); 00795 // get column in tableau 00796 factorization.updateColumn ( &work, &array ); 00797 int numberInArray=array.getNumElements(); 00798 printf("non-basic %d\n",iColumn); 00799 for (j=0;j<numberInArray;j++) { 00800 int indexValue=arrayRows[j]; 00801 double value=arrayElements[indexValue]; 00802 if (fabs(value)>1.0e-6) { 00803 printf("%d %g\n",indexValue,value); 00804 } 00805 } 00806 } 00807 #endif 00808 } 00809 } 00810 00811 delete [] rowActivity; 00812 delete [] swap; 00813 delete [] rowType; 00814 delete [] packed; 00815 delete [] rowIsBasic; 00816 delete [] columnIsBasic; 00817 return numberAdded; 00818 } |
|
Generate Mixed Integer Gomory cuts for the model of the solver interface, si. Insert the generated cuts into OsiCut, cs. There is a limit option, which will only generate cuts with less than this number of entries. We can also only look at 0-1 variables a certain distance from integer. Implements CglCutGenerator. Definition at line 26 of file CglGomory.cpp.
00028 { 00029 // Get basic problem information 00030 int numberColumns=si.getNumCols(); 00031 00032 // get integer variables and basis 00033 char * intVar = new char[numberColumns]; 00034 int i; 00035 CoinWarmStart * warmstart = si.getWarmStart(); 00036 const CoinWarmStartBasis* warm = 00037 dynamic_cast<const CoinWarmStartBasis*>(warmstart); 00038 const double * colUpper = si.getColUpper(); 00039 const double * colLower = si.getColLower(); 00040 for (i=0;i<numberColumns;i++) { 00041 if (si.isInteger(i)) { 00042 if (colUpper[i]>colLower[i]+0.5) { 00043 if (fabs(colUpper[i]-1.0)<1.0e-12&&fabs(colLower[i])<1.0e-12) { 00044 intVar[i]=1; //0-1 00045 } else if (colLower[i]>=0.0) { 00046 intVar[i] = 2; // other 00047 } else { 00048 // negative bounds - I am not sure works 00049 intVar[i] = 0; 00050 } 00051 } else { 00052 intVar[i] = 0; 00053 } 00054 } else { 00055 intVar[i]=0; 00056 } 00057 } 00058 #ifdef CGL_DEBUG 00059 const OsiRowCutDebugger * debugger = si.getRowCutDebugger(); 00060 if (debugger&&!debugger->onOptimalPath(si)) 00061 debugger = NULL; 00062 #else 00063 const OsiRowCutDebugger * debugger = NULL; 00064 #endif 00065 00066 generateCuts(debugger, cs, *si.getMatrixByCol(), 00067 si.getObjCoefficients(), si.getColSolution(), 00068 si.getColLower(), si.getColUpper(), 00069 si.getRowLower(), si.getRowUpper(), 00070 intVar,warm); 00071 00072 delete warmstart; 00073 delete [] intVar; 00074 } |
|
A function that tests the methods in the CglGomory class. The only reason for it not to be a member method is that this way it doesn't have to be compiled into the library. And that's a gain, because the library should be compiled with optimization on, but this method should be compiled with debugging. Definition at line 25 of file CglGomoryTest.cpp.
00028 { 00029 CoinRelFltEq eq(0.000001); 00030 00031 // Test default constructor 00032 { 00033 CglGomory aGenerator; 00034 assert (aGenerator.getLimit()==50); 00035 assert (aGenerator.getAway()==0.05); 00036 } 00037 00038 // Test copy & assignment etc 00039 { 00040 CglGomory rhs; 00041 { 00042 CglGomory bGenerator; 00043 bGenerator.setLimit(99); 00044 bGenerator.setAway(0.2); 00045 CglGomory cGenerator(bGenerator); 00046 rhs=bGenerator; 00047 assert (rhs.getLimit()==99); 00048 assert (rhs.getAway()==0.2); 00049 } 00050 } 00051 00052 // Test explicit form - all integer (pg 125 Wolsey) 00053 if (1) { 00054 OsiCuts osicuts; 00055 CglGomory test1; 00056 int i; 00057 int nOldCuts=0,nRowCuts; 00058 00059 // matrix data 00060 //deliberate hiccup of 2 between 0 and 1 00061 int start[5]={0,4,7,8,9}; 00062 int length[5]={2,3,1,1,1}; 00063 int rows[11]={0,2,-1,-1,0,1,2,0,1,2}; 00064 double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1}; 00065 CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length); 00066 00067 // rim data (objective not used just yet) 00068 double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0}; 00069 double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10}; 00070 double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10}; 00071 double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0}; 00072 double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0}; 00073 00074 // integer 00075 char intVar[7]={2,2,2,2,2,2,2}; 00076 00077 // basis 1 00078 int rowBasis1[3]={-1,-1,-1}; 00079 int colBasis1[5]={1,1,-1,-1,1}; 00080 CoinWarmStartBasis warm; 00081 warm.setSize(5,3); 00082 for (i=0;i<3;i++) { 00083 if (rowBasis1[i]<0) { 00084 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00085 } else { 00086 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00087 } 00088 } 00089 for (i=0;i<5;i++) { 00090 if (colBasis1[i]<0) { 00091 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00092 } else { 00093 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00094 } 00095 } 00096 00097 // solution 1 00098 double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0}; 00099 test1.generateCuts(NULL, osicuts, matrix, 00100 objective, colsol1, 00101 colLower, colUpper, 00102 rowLower, rowUpper, intVar, &warm); 00103 nRowCuts = osicuts.sizeRowCuts(); 00104 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00105 assert (nRowCuts==2); 00106 // cuts always <= 00107 int testCut=0; // test first cut as stronger 00108 double rhs=-6.0; 00109 double testCut1[5]={0.0,0.0,-1.0,-2.0,0.0}; 00110 double * cut = testCut1; 00111 double * colsol = colsol1; 00112 for (i=nOldCuts; i<nRowCuts; i++){ 00113 OsiRowCut rcut; 00114 CoinPackedVector rpv; 00115 rcut = osicuts.rowCut(i); 00116 rpv = rcut.row(); 00117 const int n = rpv.getNumElements(); 00118 const int * indices = rpv.getIndices(); 00119 double* elements = rpv.getElements(); 00120 double sum2=0.0; 00121 int k=0; 00122 double lb=rcut.lb(); 00123 double ub=rcut.ub(); 00124 for (k=0; k<n; k++){ 00125 int column=indices[k]; 00126 sum2 += colsol[column]*elements[k]; 00127 } 00128 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00129 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00130 for (k=0; k<n; k++){ 00131 int column=indices[k]; 00132 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00133 colsol[column]<<") "; 00134 } 00135 std::cout <<std::endl; 00136 } 00137 if (i-nOldCuts==testCut) { 00138 assert( eq(rhs,ub)); 00139 assert(n==2); 00140 for (k=0; k<n; k++){ 00141 int column=indices[k]; 00142 assert (eq(cut[column],elements[k])); 00143 } 00144 // add cut 00145 // explicit slack 00146 matrix.setDimensions(-1,6); 00147 rpv.insert(5,1.0*7.0); // to get cut in book 00148 rowLower[3]=ub; 00149 rowUpper[3]=ub; 00150 matrix.appendRow(rpv); 00151 } 00152 } 00153 nOldCuts=nRowCuts; 00154 // basis 2 00155 int rowBasis2[4]={-1,-1,-1,-1}; 00156 int colBasis2[6]={1,1,1,1,-1,-1}; 00157 warm.setSize(6,4); 00158 for (i=0;i<4;i++) { 00159 if (rowBasis2[i]<0) { 00160 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00161 } else { 00162 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00163 } 00164 } 00165 for (i=0;i<6;i++) { 00166 if (colBasis2[i]<0) { 00167 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00168 } else { 00169 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00170 } 00171 } 00172 00173 // solution 2 00174 double colsol2[6]={2.0,0.5,1.0,2.5,0.0,0.0}; 00175 test1.generateCuts(NULL, osicuts, matrix, 00176 objective, colsol2, 00177 colLower, colUpper, 00178 rowLower, rowUpper, intVar, &warm); 00179 nRowCuts = osicuts.sizeRowCuts(); 00180 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00181 assert (nRowCuts-nOldCuts==2); 00182 // cuts always <= 00183 testCut=0; // test first cut as stronger 00184 rhs=-1.0; 00185 double testCut2[6]={0.0,0.0,0.0,0.0,-1.0,0.0}; 00186 cut = testCut2; 00187 colsol = colsol2; 00188 for (i=nOldCuts; i<nRowCuts; i++){ 00189 OsiRowCut rcut; 00190 CoinPackedVector rpv; 00191 rcut = osicuts.rowCut(i); 00192 rpv = rcut.row(); 00193 const int n = rpv.getNumElements(); 00194 const int * indices = rpv.getIndices(); 00195 double* elements = rpv.getElements(); 00196 double sum2=0.0; 00197 int k=0; 00198 double lb=rcut.lb(); 00199 double ub=rcut.ub(); 00200 for (k=0; k<n; k++){ 00201 int column=indices[k]; 00202 sum2 += colsol[column]*elements[k]; 00203 } 00204 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00205 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00206 for (k=0; k<n; k++){ 00207 int column=indices[k]; 00208 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00209 colsol[column]<<") "; 00210 } 00211 std::cout <<std::endl; 00212 } 00213 if (i-nOldCuts==testCut) { 00214 assert( eq(rhs,ub)); 00215 assert(n==1); 00216 for (k=0; k<n; k++){ 00217 int column=indices[k]; 00218 assert (eq(cut[column],elements[k])); 00219 } 00220 // add cut 00221 // explicit slack 00222 matrix.setDimensions(-1,7); 00223 rpv.insert(6,1.0); 00224 rowLower[4]=ub; 00225 rowUpper[4]=ub; 00226 matrix.appendRow(rpv); 00227 } 00228 } 00229 nOldCuts=nRowCuts; 00230 // basis 3 00231 int rowBasis3[5]={-1,-1,-1,-1,-1}; 00232 int colBasis3[7]={1,1,1,1,1,-1,-1}; 00233 warm.setSize(7,5); 00234 for (i=0;i<5;i++) { 00235 if (rowBasis3[i]<0) { 00236 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00237 } else { 00238 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00239 } 00240 } 00241 for (i=0;i<7;i++) { 00242 if (colBasis3[i]<0) { 00243 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00244 } else { 00245 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00246 } 00247 } 00248 00249 // solution 3 00250 double colsol3[7]={2.0,1.0,2.0,2.0,1.0,0.0,0.0}; 00251 test1.generateCuts(NULL, osicuts, matrix, 00252 objective, colsol3, 00253 colLower, colUpper, 00254 rowLower, rowUpper, intVar, &warm); 00255 nRowCuts = osicuts.sizeRowCuts(); 00256 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00257 assert (nRowCuts==nOldCuts); 00258 00259 } 00260 // Test explicit form - this time with x4 flipped 00261 if (1) { 00262 OsiCuts osicuts; 00263 CglGomory test1; 00264 int i; 00265 int nOldCuts=0,nRowCuts; 00266 00267 // matrix data 00268 //deliberate hiccup of 2 between 0 and 1 00269 int start[5]={0,4,7,8,9}; 00270 int length[5]={2,3,1,1,1}; 00271 int rows[11]={0,2,-1,-1,0,1,2,0,1,2}; 00272 double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,-1,1}; 00273 CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length); 00274 00275 // rim data (objective not used just yet) 00276 double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0}; 00277 double rowLower[5]={14.0,-5.0,3.0,1.0e10,1.0e10}; 00278 double rowUpper[5]={14.0,-5.0,3.0,-1.0e10,-1.0e10}; 00279 double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0}; 00280 double colUpper[7]={100.0,100.0,100.0,8.0,100.0,100.0,100.0}; 00281 00282 // integer 00283 char intVar[7]={2,2,2,2,2,2,2}; 00284 00285 // basis 1 00286 int rowBasis1[3]={-1,-1,-1}; 00287 int colBasis1[5]={1,1,-1,-1,1}; 00288 CoinWarmStartBasis warm; 00289 warm.setSize(5,3); 00290 for (i=0;i<3;i++) { 00291 if (rowBasis1[i]<0) { 00292 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00293 } else { 00294 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00295 } 00296 } 00297 for (i=0;i<5;i++) { 00298 if (colBasis1[i]<0) { 00299 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00300 } else { 00301 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00302 } 00303 } 00304 00305 // solution 1 00306 double colsol1[5]={20.0/7.0,3.0,0.0,8.0,23.0/7.0}; 00307 test1.generateCuts(NULL, osicuts, matrix, 00308 objective, colsol1, 00309 colLower, colUpper, 00310 rowLower, rowUpper, intVar, &warm); 00311 nRowCuts = osicuts.sizeRowCuts(); 00312 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00313 assert (nRowCuts==2); 00314 // cuts always <= 00315 int testCut=0; // test first cut as stronger 00316 double rhs=10.0; 00317 double testCut1[5]={0.0,0.0,-1.0,2.0,0.0}; 00318 double * cut = testCut1; 00319 double * colsol = colsol1; 00320 for (i=nOldCuts; i<nRowCuts; i++){ 00321 OsiRowCut rcut; 00322 CoinPackedVector rpv; 00323 rcut = osicuts.rowCut(i); 00324 rpv = rcut.row(); 00325 const int n = rpv.getNumElements(); 00326 const int * indices = rpv.getIndices(); 00327 double* elements = rpv.getElements(); 00328 double sum2=0.0; 00329 int k=0; 00330 double lb=rcut.lb(); 00331 double ub=rcut.ub(); 00332 for (k=0; k<n; k++){ 00333 int column=indices[k]; 00334 sum2 += colsol[column]*elements[k]; 00335 } 00336 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00337 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00338 for (k=0; k<n; k++){ 00339 int column=indices[k]; 00340 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00341 colsol[column]<<") "; 00342 } 00343 std::cout <<std::endl; 00344 } 00345 if (i-nOldCuts==testCut) { 00346 assert( eq(rhs,ub)); 00347 assert(n==2); 00348 for (k=0; k<n; k++){ 00349 int column=indices[k]; 00350 assert (eq(cut[column],elements[k])); 00351 } 00352 // add cut 00353 // explicit slack 00354 matrix.setDimensions(-1,6); 00355 rpv.insert(5,1.0*7.0); // to get cut in book 00356 rowLower[3]=ub; 00357 rowUpper[3]=ub; 00358 matrix.appendRow(rpv); 00359 } 00360 } 00361 nOldCuts=nRowCuts; 00362 // basis 2 00363 int rowBasis2[4]={-1,-1,-1,-1}; 00364 int colBasis2[6]={1,1,1,1,-1,-1}; 00365 warm.setSize(6,4); 00366 for (i=0;i<4;i++) { 00367 if (rowBasis2[i]<0) { 00368 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00369 } else { 00370 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00371 } 00372 } 00373 for (i=0;i<6;i++) { 00374 if (colBasis2[i]<0) { 00375 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00376 } else { 00377 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00378 } 00379 } 00380 00381 // solution 2 00382 double colsol2[6]={2.0,0.5,1.0,5.5,0.0,0.0}; 00383 test1.generateCuts(NULL, osicuts, matrix, 00384 objective, colsol2, 00385 colLower, colUpper, 00386 rowLower, rowUpper, intVar, &warm); 00387 nRowCuts = osicuts.sizeRowCuts(); 00388 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00389 assert (nRowCuts-nOldCuts==2); 00390 // cuts always <= 00391 testCut=0; // test first cut as stronger 00392 rhs=-1.0; 00393 double testCut2[6]={0.0,0.0,0.0,0.0,-1.0,0.0}; 00394 cut = testCut2; 00395 colsol = colsol2; 00396 for (i=nOldCuts; i<nRowCuts; i++){ 00397 OsiRowCut rcut; 00398 CoinPackedVector rpv; 00399 rcut = osicuts.rowCut(i); 00400 rpv = rcut.row(); 00401 const int n = rpv.getNumElements(); 00402 const int * indices = rpv.getIndices(); 00403 double* elements = rpv.getElements(); 00404 double sum2=0.0; 00405 int k=0; 00406 double lb=rcut.lb(); 00407 double ub=rcut.ub(); 00408 for (k=0; k<n; k++){ 00409 int column=indices[k]; 00410 sum2 += colsol[column]*elements[k]; 00411 } 00412 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00413 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00414 for (k=0; k<n; k++){ 00415 int column=indices[k]; 00416 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00417 colsol[column]<<") "; 00418 } 00419 std::cout <<std::endl; 00420 } 00421 if (i-nOldCuts==testCut) { 00422 assert( eq(rhs,ub)); 00423 assert(n==1); 00424 for (k=0; k<n; k++){ 00425 int column=indices[k]; 00426 assert (eq(cut[column],elements[k])); 00427 } 00428 // add cut 00429 // explicit slack 00430 matrix.setDimensions(-1,7); 00431 rpv.insert(6,1.0); 00432 rowLower[4]=ub; 00433 rowUpper[4]=ub; 00434 matrix.appendRow(rpv); 00435 } 00436 } 00437 nOldCuts=nRowCuts; 00438 // basis 3 00439 int rowBasis3[5]={-1,-1,-1,-1,-1}; 00440 int colBasis3[7]={1,1,1,1,1,-1,-1}; 00441 warm.setSize(7,5); 00442 for (i=0;i<5;i++) { 00443 if (rowBasis3[i]<0) { 00444 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00445 } else { 00446 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00447 } 00448 } 00449 for (i=0;i<7;i++) { 00450 if (colBasis3[i]<0) { 00451 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00452 } else { 00453 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00454 } 00455 } 00456 00457 // solution 3 00458 double colsol3[7]={2.0,1.0,2.0,6.0,1.0,0.0,0.0}; 00459 test1.generateCuts(NULL, osicuts, matrix, 00460 objective, colsol3, 00461 colLower, colUpper, 00462 rowLower, rowUpper, intVar, &warm); 00463 nRowCuts = osicuts.sizeRowCuts(); 00464 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00465 assert (nRowCuts==nOldCuts); 00466 00467 } 00468 // Test with slacks 00469 if (1) { 00470 OsiCuts osicuts; 00471 CglGomory test1; 00472 int i; 00473 int nOldCuts=0,nRowCuts; 00474 00475 // matrix data 00476 //deliberate hiccup of 2 between 0 and 1 00477 int start[5]={0,4}; 00478 int length[5]={2,3}; 00479 int rows[11]={0,2,-1,-1,0,1,2}; 00480 double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0}; 00481 CoinPackedMatrix matrix(true,3,2,5,elements,rows,start,length); 00482 00483 // rim data (objective not used just yet) 00484 double objective[2]={-4.0,1.0}; 00485 double rowLower[5]={-1.0e10,-1.0e10,-1.0e10,1.0e10,1.0e10}; 00486 double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10}; 00487 double colLower[2]={0.0,0.0}; 00488 double colUpper[2]={100.0,100.0}; 00489 00490 // integer 00491 char intVar[2]={2,2}; 00492 00493 // basis 1 00494 int rowBasis1[3]={-1,-1,1}; 00495 int colBasis1[2]={1,1}; 00496 CoinWarmStartBasis warm; 00497 warm.setSize(2,3); 00498 for (i=0;i<3;i++) { 00499 if (rowBasis1[i]<0) { 00500 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00501 } else { 00502 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00503 } 00504 } 00505 for (i=0;i<2;i++) { 00506 if (colBasis1[i]<0) { 00507 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00508 } else { 00509 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00510 } 00511 } 00512 00513 // solution 1 00514 double colsol1[2]={20.0/7.0,3.0}; 00515 test1.generateCuts(NULL, osicuts, matrix, 00516 objective, colsol1, 00517 colLower, colUpper, 00518 rowLower, rowUpper, intVar, &warm); 00519 nRowCuts = osicuts.sizeRowCuts(); 00520 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00521 assert (nRowCuts==1); 00522 // cuts always <= 00523 int testCut=0; // test first cut as stronger 00524 double rhs=2.0; 00525 double testCut1[2]={1.0,0.0}; 00526 double * cut = testCut1; 00527 double * colsol = colsol1; 00528 for (i=nOldCuts; i<nRowCuts; i++){ 00529 OsiRowCut rcut; 00530 CoinPackedVector rpv; 00531 rcut = osicuts.rowCut(i); 00532 rpv = rcut.row(); 00533 const int n = rpv.getNumElements(); 00534 const int * indices = rpv.getIndices(); 00535 double* elements = rpv.getElements(); 00536 double sum2=0.0; 00537 int k=0; 00538 double lb=rcut.lb(); 00539 double ub=rcut.ub(); 00540 for (k=0; k<n; k++){ 00541 int column=indices[k]; 00542 sum2 += colsol[column]*elements[k]; 00543 } 00544 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00545 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00546 for (k=0; k<n; k++){ 00547 int column=indices[k]; 00548 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00549 colsol[column]<<") "; 00550 } 00551 std::cout <<std::endl; 00552 } 00553 if (i-nOldCuts==testCut) { 00554 assert( eq(rhs,ub)); 00555 assert(n==1); 00556 for (k=0; k<n; k++){ 00557 int column=indices[k]; 00558 assert (eq(cut[column],elements[k])); 00559 } 00560 // add cut 00561 rowLower[3]=-1.0e100; 00562 rowUpper[3]=ub; 00563 matrix.appendRow(rpv); 00564 } 00565 } 00566 nOldCuts=nRowCuts; 00567 // basis 2 00568 int rowBasis2[4]={1,1,-1,-1}; 00569 int colBasis2[2]={1,1}; 00570 warm.setSize(2,4); 00571 for (i=0;i<4;i++) { 00572 if (rowBasis2[i]<0) { 00573 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00574 } else { 00575 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00576 } 00577 } 00578 for (i=0;i<2;i++) { 00579 if (colBasis2[i]<0) { 00580 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00581 } else { 00582 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00583 } 00584 } 00585 00586 // solution 2 00587 double colsol2[2]={2.0,0.5}; 00588 test1.generateCuts(NULL, osicuts, matrix, 00589 objective, colsol2, 00590 colLower, colUpper, 00591 rowLower, rowUpper, intVar, &warm); 00592 nRowCuts = osicuts.sizeRowCuts(); 00593 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00594 assert (nRowCuts-nOldCuts==1); 00595 // cuts always <= 00596 testCut=0; // test first cut as stronger 00597 rhs=1.0; 00598 double testCut2[2]={1.0,-1.0}; 00599 cut = testCut2; 00600 colsol = colsol2; 00601 for (i=nOldCuts; i<nRowCuts; i++){ 00602 OsiRowCut rcut; 00603 CoinPackedVector rpv; 00604 rcut = osicuts.rowCut(i); 00605 rpv = rcut.row(); 00606 const int n = rpv.getNumElements(); 00607 const int * indices = rpv.getIndices(); 00608 double* elements = rpv.getElements(); 00609 double sum2=0.0; 00610 int k=0; 00611 double lb=rcut.lb(); 00612 double ub=rcut.ub(); 00613 for (k=0; k<n; k++){ 00614 int column=indices[k]; 00615 sum2 += colsol[column]*elements[k]; 00616 } 00617 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00618 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00619 for (k=0; k<n; k++){ 00620 int column=indices[k]; 00621 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00622 colsol[column]<<") "; 00623 } 00624 std::cout <<std::endl; 00625 } 00626 if (i-nOldCuts==testCut) { 00627 assert( eq(rhs,ub)); 00628 assert(n==2); 00629 for (k=0; k<n; k++){ 00630 int column=indices[k]; 00631 assert (eq(cut[column],elements[k])); 00632 } 00633 // add cut 00634 rowLower[4]=-1.0e100; 00635 rowUpper[4]=ub; 00636 matrix.appendRow(rpv); 00637 } 00638 } 00639 nOldCuts=nRowCuts; 00640 // basis 3 00641 int rowBasis3[5]={1,1,1,-1,-1}; 00642 int colBasis3[2]={1,1}; 00643 warm.setSize(2,5); 00644 for (i=0;i<5;i++) { 00645 if (rowBasis3[i]<0) { 00646 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00647 } else { 00648 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00649 } 00650 } 00651 for (i=0;i<2;i++) { 00652 if (colBasis3[i]<0) { 00653 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00654 } else { 00655 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00656 } 00657 } 00658 00659 // solution 3 00660 double colsol3[2]={2.0,1.0}; 00661 test1.generateCuts(NULL, osicuts, matrix, 00662 objective, colsol3, 00663 colLower, colUpper, 00664 rowLower, rowUpper, intVar, &warm); 00665 nRowCuts = osicuts.sizeRowCuts(); 00666 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00667 assert (nRowCuts==nOldCuts); 00668 00669 } 00670 // swap some rows to G 00671 if (1) { 00672 OsiCuts osicuts; 00673 CglGomory test1; 00674 int i; 00675 int nOldCuts=0,nRowCuts; 00676 00677 // matrix data 00678 //deliberate hiccup of 2 between 0 and 1 00679 int start[5]={0,4}; 00680 int length[5]={2,3}; 00681 int rows[11]={0,2,-1,-1,0,1,2}; 00682 double elements[11]={-7.0,-2.0,1.0e10,1.0e10,+2.0,1.0,+2.0}; 00683 CoinPackedMatrix matrix(true,3,2,5,elements,rows,start,length); 00684 00685 // rim data (objective not used just yet) 00686 double objective[2]={-4.0,1.0}; 00687 double rowUpper[5]={1.0e10,3.0,1.0e10,-1.0e10,-1.0e10}; 00688 double rowLower[5]={-14.0,-1.0e10,-3.0,1.0e10,1.0e10}; 00689 double colLower[2]={0.0,0.0}; 00690 double colUpper[2]={100.0,100.0}; 00691 00692 // integer 00693 char intVar[2]={2,2}; 00694 00695 // basis 1 00696 int rowBasis1[3]={-1,-1,1}; 00697 int colBasis1[2]={1,1}; 00698 CoinWarmStartBasis warm; 00699 warm.setSize(2,3); 00700 for (i=0;i<3;i++) { 00701 if (rowBasis1[i]<0) { 00702 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00703 } else { 00704 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00705 } 00706 } 00707 for (i=0;i<2;i++) { 00708 if (colBasis1[i]<0) { 00709 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00710 } else { 00711 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00712 } 00713 } 00714 00715 // solution 1 00716 double colsol1[2]={20.0/7.0,3.0}; 00717 test1.generateCuts(NULL, osicuts, matrix, 00718 objective, colsol1, 00719 colLower, colUpper, 00720 rowLower, rowUpper, intVar, &warm); 00721 nRowCuts = osicuts.sizeRowCuts(); 00722 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00723 assert (nRowCuts==1); 00724 // cuts always <= 00725 int testCut=0; // test first cut as stronger 00726 double rhs=2.0; 00727 double testCut1[2]={1.0,0.0}; 00728 double * cut = testCut1; 00729 double * colsol = colsol1; 00730 for (i=nOldCuts; i<nRowCuts; i++){ 00731 OsiRowCut rcut; 00732 CoinPackedVector rpv; 00733 rcut = osicuts.rowCut(i); 00734 rpv = rcut.row(); 00735 const int n = rpv.getNumElements(); 00736 const int * indices = rpv.getIndices(); 00737 double* elements = rpv.getElements(); 00738 double sum2=0.0; 00739 int k=0; 00740 double lb=rcut.lb(); 00741 double ub=rcut.ub(); 00742 for (k=0; k<n; k++){ 00743 int column=indices[k]; 00744 sum2 += colsol[column]*elements[k]; 00745 } 00746 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00747 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00748 for (k=0; k<n; k++){ 00749 int column=indices[k]; 00750 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00751 colsol[column]<<") "; 00752 } 00753 std::cout <<std::endl; 00754 } 00755 if (i-nOldCuts==testCut) { 00756 assert( eq(rhs,ub)); 00757 assert(n==1); 00758 for (k=0; k<n; k++){ 00759 int column=indices[k]; 00760 assert (eq(cut[column],elements[k])); 00761 } 00762 // add cut 00763 rowLower[3]=-1.0e100; 00764 rowUpper[3]=ub; 00765 matrix.appendRow(rpv); 00766 } 00767 } 00768 nOldCuts=nRowCuts; 00769 // basis 2 00770 int rowBasis2[4]={1,1,-1,-1}; 00771 int colBasis2[2]={1,1}; 00772 warm.setSize(2,4); 00773 for (i=0;i<4;i++) { 00774 if (rowBasis2[i]<0) { 00775 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00776 } else { 00777 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00778 } 00779 } 00780 for (i=0;i<2;i++) { 00781 if (colBasis2[i]<0) { 00782 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00783 } else { 00784 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00785 } 00786 } 00787 00788 // solution 2 00789 double colsol2[2]={2.0,0.5}; 00790 test1.generateCuts(NULL, osicuts, matrix, 00791 objective, colsol2, 00792 colLower, colUpper, 00793 rowLower, rowUpper, intVar, &warm); 00794 nRowCuts = osicuts.sizeRowCuts(); 00795 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00796 assert (nRowCuts-nOldCuts==1); 00797 // cuts always <= 00798 testCut=0; // test first cut as stronger 00799 rhs=1.0; 00800 double testCut2[2]={1.0,-1.0}; 00801 cut = testCut2; 00802 colsol = colsol2; 00803 for (i=nOldCuts; i<nRowCuts; i++){ 00804 OsiRowCut rcut; 00805 CoinPackedVector rpv; 00806 rcut = osicuts.rowCut(i); 00807 rpv = rcut.row(); 00808 const int n = rpv.getNumElements(); 00809 const int * indices = rpv.getIndices(); 00810 double* elements = rpv.getElements(); 00811 double sum2=0.0; 00812 int k=0; 00813 double lb=rcut.lb(); 00814 double ub=rcut.ub(); 00815 for (k=0; k<n; k++){ 00816 int column=indices[k]; 00817 sum2 += colsol[column]*elements[k]; 00818 } 00819 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00820 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00821 for (k=0; k<n; k++){ 00822 int column=indices[k]; 00823 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00824 colsol[column]<<") "; 00825 } 00826 std::cout <<std::endl; 00827 } 00828 if (i-nOldCuts==testCut) { 00829 assert( eq(rhs,ub)); 00830 assert(n==2); 00831 for (k=0; k<n; k++){ 00832 int column=indices[k]; 00833 assert (eq(cut[column],elements[k])); 00834 } 00835 // add cut 00836 rowLower[4]=-1.0e100; 00837 rowUpper[4]=ub; 00838 matrix.appendRow(rpv); 00839 } 00840 } 00841 nOldCuts=nRowCuts; 00842 // basis 3 00843 int rowBasis3[5]={1,1,1,-1,-1}; 00844 int colBasis3[2]={1,1}; 00845 warm.setSize(2,5); 00846 for (i=0;i<5;i++) { 00847 if (rowBasis3[i]<0) { 00848 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00849 } else { 00850 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00851 } 00852 } 00853 for (i=0;i<2;i++) { 00854 if (colBasis3[i]<0) { 00855 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00856 } else { 00857 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00858 } 00859 } 00860 00861 // solution 3 00862 double colsol3[2]={2.0,1.0}; 00863 test1.generateCuts(NULL, osicuts, matrix, 00864 objective, colsol3, 00865 colLower, colUpper, 00866 rowLower, rowUpper, intVar, &warm); 00867 nRowCuts = osicuts.sizeRowCuts(); 00868 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00869 assert (nRowCuts==nOldCuts); 00870 00871 } 00872 00873 00874 // NOW mixed integer gomory cuts 00875 00876 // Test explicit form - (pg 130 Wolsey) 00877 // Some arrays left same size as previously although not used in full 00878 if (1) { 00879 OsiCuts osicuts; 00880 CglGomory test1; 00881 int i; 00882 int nOldCuts=0,nRowCuts; 00883 00884 // matrix data 00885 //deliberate hiccup of 2 between 0 and 1 00886 int start[5]={0,4,7,8,9}; 00887 int length[5]={2,3,1,1,1}; 00888 int rows[11]={0,2,-1,-1,0,1,2,0,1,2}; 00889 double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1}; 00890 CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length); 00891 00892 // rim data (objective not used just yet) 00893 double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0}; 00894 double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10}; 00895 double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10}; 00896 double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0}; 00897 double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0}; 00898 00899 // integer 00900 char intVar[7]={2,0,0,0,0,0,0}; 00901 00902 // basis 1 00903 int rowBasis1[3]={-1,-1,-1}; 00904 int colBasis1[5]={1,1,-1,-1,1}; 00905 CoinWarmStartBasis warm; 00906 warm.setSize(5,3); 00907 for (i=0;i<3;i++) { 00908 if (rowBasis1[i]<0) { 00909 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00910 } else { 00911 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00912 } 00913 } 00914 for (i=0;i<5;i++) { 00915 if (colBasis1[i]<0) { 00916 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00917 } else { 00918 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00919 } 00920 } 00921 00922 // solution 1 00923 double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0}; 00924 test1.generateCuts(NULL, osicuts, matrix, 00925 objective, colsol1, 00926 colLower, colUpper, 00927 rowLower, rowUpper, intVar, &warm); 00928 nRowCuts = osicuts.sizeRowCuts(); 00929 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 00930 assert (nRowCuts==1); 00931 // cuts always <= 00932 int testCut=0; // test first cut as stronger 00933 double rhs=-6.0/7.0; 00934 double testCut1[5]={0.0,0.0,-1.0/7.0,-2.0/7.0,0.0}; 00935 double * cut = testCut1; 00936 double * colsol = colsol1; 00937 for (i=nOldCuts; i<nRowCuts; i++){ 00938 OsiRowCut rcut; 00939 CoinPackedVector rpv; 00940 rcut = osicuts.rowCut(i); 00941 rpv = rcut.row(); 00942 const int n = rpv.getNumElements(); 00943 const int * indices = rpv.getIndices(); 00944 double* elements = rpv.getElements(); 00945 double sum2=0.0; 00946 int k=0; 00947 double lb=rcut.lb(); 00948 double ub=rcut.ub(); 00949 for (k=0; k<n; k++){ 00950 int column=indices[k]; 00951 sum2 += colsol[column]*elements[k]; 00952 } 00953 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 00954 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 00955 for (k=0; k<n; k++){ 00956 int column=indices[k]; 00957 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 00958 colsol[column]<<") "; 00959 } 00960 std::cout <<std::endl; 00961 } 00962 if (i-nOldCuts==testCut) { 00963 assert( eq(rhs,ub)); 00964 assert(n==2); 00965 for (k=0; k<n; k++){ 00966 int column=indices[k]; 00967 assert (eq(cut[column],elements[k])); 00968 } 00969 // add cut 00970 // explicit slack 00971 matrix.setDimensions(-1,6); 00972 rpv.insert(5,1.0); // to get cut in book 00973 rowLower[3]=ub; 00974 rowUpper[3]=ub; 00975 matrix.appendRow(rpv); 00976 } 00977 } 00978 nOldCuts=nRowCuts; 00979 // basis 2 00980 int rowBasis2[4]={-1,-1,-1,-1}; 00981 int colBasis2[6]={1,1,1,1,-1,-1}; 00982 warm.setSize(6,4); 00983 for (i=0;i<4;i++) { 00984 if (rowBasis2[i]<0) { 00985 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 00986 } else { 00987 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 00988 } 00989 } 00990 for (i=0;i<6;i++) { 00991 if (colBasis2[i]<0) { 00992 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 00993 } else { 00994 warm.setStructStatus(i,CoinWarmStartBasis::basic); 00995 } 00996 } 00997 00998 // solution 2 00999 double colsol2[6]={2.0,0.5,1.0,2.5,0.0,0.0}; 01000 test1.generateCuts(NULL, osicuts, matrix, 01001 objective, colsol2, 01002 colLower, colUpper, 01003 rowLower, rowUpper, intVar, &warm); 01004 nRowCuts = osicuts.sizeRowCuts(); 01005 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 01006 assert (nRowCuts==nOldCuts); 01007 01008 } 01009 // Test explicit form - this time with x4 flipped 01010 if (1) { 01011 OsiCuts osicuts; 01012 CglGomory test1; 01013 int i; 01014 int nOldCuts=0,nRowCuts; 01015 01016 // matrix data 01017 //deliberate hiccup of 2 between 0 and 1 01018 int start[5]={0,4,7,8,9}; 01019 int length[5]={2,3,1,1,1}; 01020 int rows[11]={0,2,-1,-1,0,1,2,0,1,2}; 01021 double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,-1,1}; 01022 CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length); 01023 01024 // rim data (objective not used just yet) 01025 double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0}; 01026 double rowLower[5]={14.0,-5.0,3.0,1.0e10,1.0e10}; 01027 double rowUpper[5]={14.0,-5.0,3.0,-1.0e10,-1.0e10}; 01028 double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0}; 01029 double colUpper[7]={100.0,100.0,100.0,8.0,100.0,100.0,100.0}; 01030 01031 // integer 01032 char intVar[7]={2,0,0,0,0,0,0}; 01033 01034 // basis 1 01035 int rowBasis1[3]={-1,-1,-1}; 01036 int colBasis1[5]={1,1,-1,-1,1}; 01037 CoinWarmStartBasis warm; 01038 warm.setSize(5,3); 01039 for (i=0;i<3;i++) { 01040 if (rowBasis1[i]<0) { 01041 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 01042 } else { 01043 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 01044 } 01045 } 01046 for (i=0;i<5;i++) { 01047 if (colBasis1[i]<0) { 01048 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 01049 } else { 01050 warm.setStructStatus(i,CoinWarmStartBasis::basic); 01051 } 01052 } 01053 01054 // solution 1 01055 double colsol1[5]={20.0/7.0,3.0,0.0,8.0,23.0/7.0}; 01056 test1.generateCuts(NULL, osicuts, matrix, 01057 objective, colsol1, 01058 colLower, colUpper, 01059 rowLower, rowUpper, intVar, &warm); 01060 nRowCuts = osicuts.sizeRowCuts(); 01061 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 01062 assert (nRowCuts==1); 01063 // cuts always <= 01064 int testCut=0; 01065 double rhs=10.0/7.0; 01066 double testCut1[5]={0.0,0.0,-1.0/7.0,2.0/7.0,0.0}; 01067 double * cut = testCut1; 01068 double * colsol = colsol1; 01069 for (i=nOldCuts; i<nRowCuts; i++){ 01070 OsiRowCut rcut; 01071 CoinPackedVector rpv; 01072 rcut = osicuts.rowCut(i); 01073 rpv = rcut.row(); 01074 const int n = rpv.getNumElements(); 01075 const int * indices = rpv.getIndices(); 01076 double* elements = rpv.getElements(); 01077 double sum2=0.0; 01078 int k=0; 01079 double lb=rcut.lb(); 01080 double ub=rcut.ub(); 01081 for (k=0; k<n; k++){ 01082 int column=indices[k]; 01083 sum2 += colsol[column]*elements[k]; 01084 } 01085 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 01086 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 01087 for (k=0; k<n; k++){ 01088 int column=indices[k]; 01089 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 01090 colsol[column]<<") "; 01091 } 01092 std::cout <<std::endl; 01093 } 01094 if (i-nOldCuts==testCut) { 01095 assert( eq(rhs,ub)); 01096 assert(n==2); 01097 for (k=0; k<n; k++){ 01098 int column=indices[k]; 01099 assert (eq(cut[column],elements[k])); 01100 } 01101 // add cut 01102 // explicit slack 01103 matrix.setDimensions(-1,6); 01104 rpv.insert(5,1.0); // to get cut in book 01105 rowLower[3]=ub; 01106 rowUpper[3]=ub; 01107 matrix.appendRow(rpv); 01108 } 01109 } 01110 nOldCuts=nRowCuts; 01111 // basis 2 01112 int rowBasis2[4]={-1,-1,-1,-1}; 01113 int colBasis2[6]={1,1,1,1,-1,-1}; 01114 warm.setSize(6,4); 01115 for (i=0;i<4;i++) { 01116 if (rowBasis2[i]<0) { 01117 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 01118 } else { 01119 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 01120 } 01121 } 01122 for (i=0;i<6;i++) { 01123 if (colBasis2[i]<0) { 01124 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 01125 } else { 01126 warm.setStructStatus(i,CoinWarmStartBasis::basic); 01127 } 01128 } 01129 01130 // solution 2 01131 double colsol2[6]={2.0,0.5,1.0,5.5,0.0,0.0}; 01132 test1.generateCuts(NULL, osicuts, matrix, 01133 objective, colsol2, 01134 colLower, colUpper, 01135 rowLower, rowUpper, intVar, &warm); 01136 nRowCuts = osicuts.sizeRowCuts(); 01137 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 01138 assert (nRowCuts==nOldCuts); 01139 01140 } 01141 // Test with slacks 01142 if (1) { 01143 OsiCuts osicuts; 01144 CglGomory test1; 01145 int i; 01146 int nOldCuts=0,nRowCuts; 01147 01148 // matrix data 01149 //deliberate hiccup of 2 between 0 and 1 01150 int start[5]={0,4}; 01151 int length[5]={2,3}; 01152 int rows[11]={0,2,-1,-1,0,1,2}; 01153 double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0}; 01154 CoinPackedMatrix matrix(true,3,2,5,elements,rows,start,length); 01155 01156 // rim data (objective not used just yet) 01157 double objective[2]={-4.0,1.0}; 01158 double rowLower[5]={-1.0e10,-1.0e10,-1.0e10,1.0e10,1.0e10}; 01159 double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10}; 01160 double colLower[2]={0.0,0.0}; 01161 double colUpper[2]={100.0,100.0}; 01162 01163 // integer 01164 char intVar[2]={2,0}; 01165 01166 // basis 1 01167 int rowBasis1[3]={-1,-1,1}; 01168 int colBasis1[2]={1,1}; 01169 CoinWarmStartBasis warm; 01170 warm.setSize(2,3); 01171 for (i=0;i<3;i++) { 01172 if (rowBasis1[i]<0) { 01173 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 01174 } else { 01175 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 01176 } 01177 } 01178 for (i=0;i<2;i++) { 01179 if (colBasis1[i]<0) { 01180 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 01181 } else { 01182 warm.setStructStatus(i,CoinWarmStartBasis::basic); 01183 } 01184 } 01185 01186 // solution 1 01187 double colsol1[2]={20.0/7.0,3.0}; 01188 test1.generateCuts(NULL, osicuts, matrix, 01189 objective, colsol1, 01190 colLower, colUpper, 01191 rowLower, rowUpper, intVar, &warm); 01192 nRowCuts = osicuts.sizeRowCuts(); 01193 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 01194 assert (nRowCuts==1); 01195 // cuts always <= 01196 int testCut=0; // test first cut as stronger 01197 double rhs=2.0; 01198 double testCut1[2]={1.0,0.0}; 01199 double * cut = testCut1; 01200 double * colsol = colsol1; 01201 for (i=nOldCuts; i<nRowCuts; i++){ 01202 OsiRowCut rcut; 01203 CoinPackedVector rpv; 01204 rcut = osicuts.rowCut(i); 01205 rpv = rcut.row(); 01206 const int n = rpv.getNumElements(); 01207 const int * indices = rpv.getIndices(); 01208 double* elements = rpv.getElements(); 01209 double sum2=0.0; 01210 int k=0; 01211 double lb=rcut.lb(); 01212 double ub=rcut.ub(); 01213 for (k=0; k<n; k++){ 01214 int column=indices[k]; 01215 sum2 += colsol[column]*elements[k]; 01216 } 01217 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 01218 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 01219 for (k=0; k<n; k++){ 01220 int column=indices[k]; 01221 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 01222 colsol[column]<<") "; 01223 } 01224 std::cout <<std::endl; 01225 } 01226 if (i-nOldCuts==testCut) { 01227 assert( eq(rhs,ub)); 01228 assert(n==1); 01229 for (k=0; k<n; k++){ 01230 int column=indices[k]; 01231 assert (eq(cut[column],elements[k])); 01232 } 01233 // add cut 01234 rowLower[3]=-1.0e100; 01235 rowUpper[3]=ub; 01236 matrix.appendRow(rpv); 01237 } 01238 } 01239 nOldCuts=nRowCuts; 01240 // basis 2 01241 int rowBasis2[4]={1,1,-1,-1}; 01242 int colBasis2[2]={1,1}; 01243 warm.setSize(2,4); 01244 for (i=0;i<4;i++) { 01245 if (rowBasis2[i]<0) { 01246 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 01247 } else { 01248 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 01249 } 01250 } 01251 for (i=0;i<2;i++) { 01252 if (colBasis2[i]<0) { 01253 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 01254 } else { 01255 warm.setStructStatus(i,CoinWarmStartBasis::basic); 01256 } 01257 } 01258 01259 // solution 2 01260 double colsol2[2]={2.0,0.5}; 01261 test1.generateCuts(NULL, osicuts, matrix, 01262 objective, colsol2, 01263 colLower, colUpper, 01264 rowLower, rowUpper, intVar, &warm); 01265 nRowCuts = osicuts.sizeRowCuts(); 01266 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 01267 assert (nRowCuts==nOldCuts); 01268 01269 } 01270 // swap some rows to G 01271 if (1) { 01272 OsiCuts osicuts; 01273 CglGomory test1; 01274 int i; 01275 int nOldCuts=0,nRowCuts; 01276 01277 // matrix data 01278 //deliberate hiccup of 2 between 0 and 1 01279 int start[5]={0,4}; 01280 int length[5]={2,3}; 01281 int rows[11]={0,2,-1,-1,0,1,2}; 01282 double elements[11]={-7.0,-2.0,1.0e10,1.0e10,+2.0,1.0,+2.0}; 01283 CoinPackedMatrix matrix(true,3,2,5,elements,rows,start,length); 01284 01285 // rim data (objective not used just yet) 01286 double objective[2]={-4.0,1.0}; 01287 double rowUpper[5]={1.0e10,3.0,1.0e10,-1.0e10,-1.0e10}; 01288 double rowLower[5]={-14.0,-1.0e10,-3.0,1.0e10,1.0e10}; 01289 double colLower[2]={0.0,0.0}; 01290 double colUpper[2]={100.0,100.0}; 01291 01292 // integer 01293 char intVar[2]={2,0}; 01294 01295 // basis 1 01296 int rowBasis1[3]={-1,-1,1}; 01297 int colBasis1[2]={1,1}; 01298 CoinWarmStartBasis warm; 01299 warm.setSize(2,3); 01300 for (i=0;i<3;i++) { 01301 if (rowBasis1[i]<0) { 01302 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 01303 } else { 01304 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 01305 } 01306 } 01307 for (i=0;i<2;i++) { 01308 if (colBasis1[i]<0) { 01309 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 01310 } else { 01311 warm.setStructStatus(i,CoinWarmStartBasis::basic); 01312 } 01313 } 01314 01315 // solution 1 01316 double colsol1[2]={20.0/7.0,3.0}; 01317 test1.generateCuts(NULL, osicuts, matrix, 01318 objective, colsol1, 01319 colLower, colUpper, 01320 rowLower, rowUpper, intVar, &warm); 01321 nRowCuts = osicuts.sizeRowCuts(); 01322 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 01323 assert (nRowCuts==1); 01324 // cuts always <= 01325 int testCut=0; // test first cut as stronger 01326 double rhs=2.0; 01327 double testCut1[2]={1.0,0.0}; 01328 double * cut = testCut1; 01329 double * colsol = colsol1; 01330 for (i=nOldCuts; i<nRowCuts; i++){ 01331 OsiRowCut rcut; 01332 CoinPackedVector rpv; 01333 rcut = osicuts.rowCut(i); 01334 rpv = rcut.row(); 01335 const int n = rpv.getNumElements(); 01336 const int * indices = rpv.getIndices(); 01337 double* elements = rpv.getElements(); 01338 double sum2=0.0; 01339 int k=0; 01340 double lb=rcut.lb(); 01341 double ub=rcut.ub(); 01342 for (k=0; k<n; k++){ 01343 int column=indices[k]; 01344 sum2 += colsol[column]*elements[k]; 01345 } 01346 if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) { 01347 std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl; 01348 for (k=0; k<n; k++){ 01349 int column=indices[k]; 01350 std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<< 01351 colsol[column]<<") "; 01352 } 01353 std::cout <<std::endl; 01354 } 01355 if (i-nOldCuts==testCut) { 01356 assert( eq(rhs,ub)); 01357 assert(n==1); 01358 for (k=0; k<n; k++){ 01359 int column=indices[k]; 01360 assert (eq(cut[column],elements[k])); 01361 } 01362 // add cut 01363 rowLower[3]=-1.0e100; 01364 rowUpper[3]=ub; 01365 matrix.appendRow(rpv); 01366 } 01367 } 01368 nOldCuts=nRowCuts; 01369 // basis 2 01370 int rowBasis2[4]={1,1,-1,-1}; 01371 int colBasis2[2]={1,1}; 01372 warm.setSize(2,4); 01373 for (i=0;i<4;i++) { 01374 if (rowBasis2[i]<0) { 01375 warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound); 01376 } else { 01377 warm.setArtifStatus(i,CoinWarmStartBasis::basic); 01378 } 01379 } 01380 for (i=0;i<2;i++) { 01381 if (colBasis2[i]<0) { 01382 warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound); 01383 } else { 01384 warm.setStructStatus(i,CoinWarmStartBasis::basic); 01385 } 01386 } 01387 01388 // solution 2 01389 double colsol2[2]={2.0,0.5}; 01390 test1.generateCuts(NULL, osicuts, matrix, 01391 objective, colsol2, 01392 colLower, colUpper, 01393 rowLower, rowUpper, intVar, &warm); 01394 nRowCuts = osicuts.sizeRowCuts(); 01395 std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl; 01396 assert (nRowCuts==nOldCuts); 01397 01398 } 01399 01400 // Miplib3 problem p0033 01401 if (1) { 01402 // Setup 01403 OsiSolverInterface * siP = baseSiP->clone(); 01404 std::string fn(mpsDir+"p0033"); 01405 siP->readMps(fn.c_str(),"mps"); 01406 siP->activateRowCutDebugger("p0033"); 01407 CglGomory test; 01408 01409 // Solve the LP relaxation of the model and 01410 // print out ofv for sake of comparison 01411 siP->initialSolve(); 01412 double lpRelaxBefore=siP->getObjValue(); 01413 assert( eq(lpRelaxBefore, 2520.5717391304347) ); 01414 double mycs[] = {0, 1, 0, 0, -2.0837010502455788e-19, 1, 0, 0, 1, 01415 0.021739130434782594, 0.35652173913043478, 01416 -6.7220534694101275e-18, 5.3125906451789717e-18, 01417 1, 0, 1.9298798670241979e-17, 0, 0, 0, 01418 7.8875708048320448e-18, 0.5, 0, 01419 0.85999999999999999, 1, 1, 0.57999999999999996, 01420 1, 0, 1, 0, 0.25, 0, 0.67500000000000004}; 01421 siP->setColSolution(mycs); 01422 #ifdef CGL_DEBUG 01423 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 01424 #endif 01425 01426 OsiCuts cuts; 01427 01428 // Test generateCuts method 01429 test.generateCuts(*siP,cuts); 01430 OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts); 01431 01432 siP->resolve(); 01433 double lpRelaxAfter=siP->getObjValue(); 01434 //assert( eq(lpRelaxAfter, 2592.1908295194507) ); 01435 assert( eq(lpRelaxAfter, 2582.6167554453768) ); 01436 #ifdef CGL_DEBUG 01437 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 01438 printf("\n\nFinal LP min=%f\n",lpRelaxAfter); 01439 #endif 01440 assert( lpRelaxBefore < lpRelaxAfter ); 01441 01442 delete siP; 01443 } 01444 //exit(0); 01445 01446 } |