Changeset 1878 for trunk/Clp/src/ClpPresolve.cpp
 Timestamp:
 Aug 30, 2012 11:43:19 AM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpPresolve.cpp
r1834 r1878 13 13 14 14 #include "CoinHelperFunctions.hpp" 15 #if CLP_HAS_ABC 16 #include "CoinAbcCommon.hpp" 17 #endif 15 18 16 19 #include "CoinPackedMatrix.hpp" … … 55 58 nrows_(0), 56 59 nelems_(0), 60 #ifdef ABC_INHERIT 61 numberPasses_(20), 62 #else 57 63 numberPasses_(5), 64 #endif 58 65 substitution_(3), 59 66 #ifndef CLP_NO_STD … … 274 281 #ifndef CLP_NO_STD 275 282 } 276 #endif 283 #endif 277 284 // put back duals 278 285 CoinMemcpyN(prob.rowduals_, nrows_, originalModel_>dualRowSolution()); … … 444 451 } 445 452 #endif 453 static int tightenDoubletons2(CoinPresolveMatrix * prob) 454 { 455 // columnmajor representation 456 const int ncols = prob>ncols_ ; 457 const CoinBigIndex *const mcstrt = prob>mcstrt_ ; 458 const int *const hincol = prob>hincol_ ; 459 const int *const hrow = prob>hrow_ ; 460 double * colels = prob>colels_ ; 461 double * cost = prob>cost_ ; 462 463 // column type, bounds, solution, and status 464 const unsigned char *const integerType = prob>integerType_ ; 465 double * clo = prob>clo_ ; 466 double * cup = prob>cup_ ; 467 // rowmajor representation 468 //const int nrows = prob>nrows_ ; 469 const CoinBigIndex *const mrstrt = prob>mrstrt_ ; 470 const int *const hinrow = prob>hinrow_ ; 471 const int *const hcol = prob>hcol_ ; 472 double * rowels = prob>rowels_ ; 473 474 // row bounds 475 double *const rlo = prob>rlo_ ; 476 double *const rup = prob>rup_ ; 477 478 // tolerances 479 //const double ekkinf2 = PRESOLVE_SMALL_INF ; 480 //const double ekkinf = ekkinf2*1.0e8 ; 481 //const double ztolcbarj = prob>ztoldj_ ; 482 //const CoinRelFltEq relEq(prob>ztolzb_) ; 483 int numberChanged=0; 484 double bound[2]; 485 double alpha[2]={0.0,0.0}; 486 double offset=0.0; 487 488 for (int icol=0;icol<ncols;icol++) { 489 if (hincol[icol]==2) { 490 CoinBigIndex start=mcstrt[icol]; 491 int row0 = hrow[start]; 492 if (hinrow[row0]!=2) 493 continue; 494 int row1 = hrow[start+1]; 495 if (hinrow[row1]!=2) 496 continue; 497 double element0 = colels[start]; 498 double rowUpper0=rup[row0]; 499 bool swapSigns0=false; 500 if (rlo[row0]>1.0e30) { 501 if (rup[row0]>1.0e30) { 502 swapSigns0=true; 503 rowUpper0=rlo[row0]; 504 element0=element0; 505 } else { 506 // range or equality 507 continue; 508 } 509 } else if (rup[row0]>1.0e30) { 510 // free 511 continue; 512 } 513 #if 0 514 // skip here for speed 515 // skip if no cost (should be able to get rid of) 516 if (!cost[icol]) { 517 printf("should be able to get rid of %d with no cost\n",icol); 518 continue; 519 } 520 // skip if negative cost for now 521 if (cost[icol]<0.0) { 522 printf("code for negative cost\n"); 523 continue; 524 } 525 #endif 526 double element1 = colels[start+1]; 527 double rowUpper1=rup[row1]; 528 bool swapSigns1=false; 529 if (rlo[row1]>1.0e30) { 530 if (rup[row1]>1.0e30) { 531 swapSigns1=true; 532 rowUpper1=rlo[row1]; 533 element1=element1; 534 } else { 535 // range or equality 536 continue; 537 } 538 } else if (rup[row1]>1.0e30) { 539 // free 540 continue; 541 } 542 double lowerX=clo[icol]; 543 double upperX=cup[icol]; 544 int otherCol=1; 545 CoinBigIndex startRow=mrstrt[row0]; 546 for (CoinBigIndex j=startRow;j<startRow+2;j++) { 547 int jcol=hcol[j]; 548 if (jcol!=icol) { 549 alpha[0]=swapSigns0 ? rowels[j] :rowels[j]; 550 otherCol=jcol; 551 } 552 } 553 startRow=mrstrt[row1]; 554 bool possible=true; 555 for (CoinBigIndex j=startRow;j<startRow+2;j++) { 556 int jcol=hcol[j]; 557 if (jcol!=icol) { 558 if (jcol==otherCol) { 559 alpha[1]=swapSigns1 ? rowels[j] :rowels[j]; 560 } else { 561 possible=false; 562 } 563 } 564 } 565 if (possible) { 566 // skip if no cost (should be able to get rid of) 567 if (!cost[icol]) { 568 printf("should be able to get rid of %d with no cost\n",icol); 569 continue; 570 } 571 // skip if negative cost for now 572 if (cost[icol]<0.0) { 573 printf("code for negative cost\n"); 574 continue; 575 } 576 bound[0]=clo[otherCol]; 577 bound[1]=cup[otherCol]; 578 double lowestLowest=COIN_DBL_MAX; 579 double highestLowest=COIN_DBL_MAX; 580 double lowestHighest=COIN_DBL_MAX; 581 double highestHighest=COIN_DBL_MAX; 582 int binding0=0; 583 int binding1=0; 584 for (int k=0;k<2;k++) { 585 bool infLow0=false; 586 bool infLow1=false; 587 double sum0=0.0; 588 double sum1=0.0; 589 double value=bound[k]; 590 if (fabs(value)<1.0e30) { 591 sum0+=alpha[0]*value; 592 sum1+=alpha[1]*value; 593 } else { 594 if (alpha[0]>0.0) { 595 if (value<0.0) 596 infLow0 =true; 597 } else if (alpha[0]<0.0) { 598 if (value>0.0) 599 infLow0 =true; 600 } 601 if (alpha[1]>0.0) { 602 if (value<0.0) 603 infLow1 =true; 604 } else if (alpha[1]<0.0) { 605 if (value>0.0) 606 infLow1 =true; 607 } 608 } 609 /* Got sums 610 */ 611 double thisLowest0=COIN_DBL_MAX; 612 double thisHighest0=COIN_DBL_MAX; 613 if (element0>0.0) { 614 // upper bound unless inf&2 !=0 615 if (!infLow0) 616 thisHighest0 = (rowUpper0sum0)/element0; 617 } else { 618 // lower bound unless inf&2 !=0 619 if (!infLow0) 620 thisLowest0 = (rowUpper0sum0)/element0; 621 } 622 double thisLowest1=COIN_DBL_MAX; 623 double thisHighest1=COIN_DBL_MAX; 624 if (element1>0.0) { 625 // upper bound unless inf&2 !=0 626 if (!infLow1) 627 thisHighest1 = (rowUpper1sum1)/element1; 628 } else { 629 // lower bound unless inf&2 !=0 630 if (!infLow1) 631 thisLowest1 = (rowUpper1sum1)/element1; 632 } 633 if (thisLowest0>thisLowest1+1.0e12) { 634 if (thisLowest0>lowerX+1.0e12) 635 binding0= 1<<k; 636 } else if (thisLowest1>thisLowest0+1.0e12) { 637 if (thisLowest1>lowerX+1.0e12) 638 binding1= 1<<k; 639 thisLowest0=thisLowest1; 640 } 641 if (thisHighest0<thisHighest11.0e12) { 642 if (thisHighest0<upperX1.0e12) 643 binding0= 1<<k; 644 } else if (thisHighest1<thisHighest01.0e12) { 645 if (thisHighest1<upperX1.0e12) 646 binding1= 1<<k; 647 thisHighest0=thisHighest1; 648 } 649 lowestLowest=CoinMin(lowestLowest,thisLowest0); 650 highestHighest=CoinMax(highestHighest,thisHighest0); 651 lowestHighest=CoinMin(lowestHighest,thisHighest0); 652 highestLowest=CoinMax(highestLowest,thisLowest0); 653 } 654 // see if any good 655 //#define PRINT_VALUES 656 if (!binding0!binding1) { 657 printf("Row redundant for column %d\n",icol); 658 } else { 659 #ifdef PRINT_VALUES 660 printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n", 661 icol,lowerX,upperX,lowestLowest,lowestHighest, 662 highestLowest,highestHighest); 663 #endif 664 // if integer adjust 665 if (integerType[icol]) { 666 lowestLowest=ceil(lowestLowest1.0e5); 667 highestLowest=ceil(highestLowest1.0e5); 668 lowestHighest=floor(lowestHighest+1.0e5); 669 highestHighest=floor(highestHighest+1.0e5); 670 } 671 // if costed may be able to adjust 672 if (cost[icol]>=0.0) { 673 if (highestLowest<upperX&&highestLowest>=lowerX&&highestHighest<1.0e30) { 674 highestHighest=CoinMin(highestHighest,highestLowest); 675 } 676 } 677 if (cost[icol]<=0.0) { 678 if (lowestHighest>lowerX&&lowestHighest<=upperX&&lowestHighest>1.0e30) { 679 lowestLowest=CoinMax(lowestLowest,lowestHighest); 680 } 681 } 682 #if 1 683 if (lowestLowest>lowerX+1.0e8) { 684 #ifdef PRINT_VALUES 685 printf("Can increase lower bound on %d from %g to %g\n", 686 icol,lowerX,lowestLowest); 687 #endif 688 lowerX=lowestLowest; 689 } 690 if (highestHighest<upperX1.0e8) { 691 #ifdef PRINT_VALUES 692 printf("Can decrease upper bound on %d from %g to %g\n", 693 icol,upperX,highestHighest); 694 #endif 695 upperX=highestHighest; 696 697 } 698 #endif 699 // see if we can move costs 700 double xValue; 701 double yValue0; 702 double yValue1; 703 double newLower=COIN_DBL_MAX; 704 double newUpper=COIN_DBL_MAX; 705 double ranges0[2]; 706 double ranges1[2]; 707 double costEqual; 708 double slope[2]; 709 assert (binding0+binding1==3); 710 // get where equal 711 xValue=(rowUpper0*element1rowUpper1*element0)/(alpha[0]*element1alpha[1]*element0); 712 yValue0=(rowUpper0xValue*alpha[0])/element0; 713 yValue1=(rowUpper1xValue*alpha[1])/element1; 714 newLower=CoinMin(newLower,CoinMax(yValue0,yValue1)); 715 newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1)); 716 double xValueEqual=xValue; 717 double yValueEqual=yValue0; 718 costEqual = xValue*cost[otherCol]+yValueEqual*cost[icol]; 719 if (binding0==1) { 720 ranges0[0]=bound[0]; 721 ranges0[1]=yValue0; 722 ranges1[0]=yValue0; 723 ranges1[1]=bound[1]; 724 // take x 1.0 down 725 double x=xValue1.0; 726 double y=(rowUpper0x*alpha[0])/element0; 727 double costTotal = x*cost[otherCol]+y*cost[icol]; 728 slope[0] = costEqualcostTotal; 729 // take x 1.0 up 730 x=xValue+1.0; 731 y=(rowUpper1x*alpha[1])/element0; 732 costTotal = x*cost[otherCol]+y*cost[icol]; 733 slope[1] = costTotalcostEqual; 734 } else { 735 ranges1[0]=bound[0]; 736 ranges1[1]=yValue0; 737 ranges0[0]=yValue0; 738 ranges0[1]=bound[1]; 739 // take x 1.0 down 740 double x=xValue1.0; 741 double y=(rowUpper1x*alpha[1])/element0; 742 double costTotal = x*cost[otherCol]+y*cost[icol]; 743 slope[1] = costEqualcostTotal; 744 // take x 1.0 up 745 x=xValue+1.0; 746 y=(rowUpper0x*alpha[0])/element0; 747 costTotal = x*cost[otherCol]+y*cost[icol]; 748 slope[0] = costTotalcostEqual; 749 } 750 #ifdef PRINT_VALUES 751 printf("equal value of %d is %g, value of %d is max(%g,%g)  %g\n", 752 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1)); 753 printf("Cost at equality %g for constraint 0 ranges %g > %g slope %g for constraint 1 ranges %g > %g slope %g\n", 754 costEqual,ranges0[0],ranges0[1],slope[0],ranges1[0],ranges1[1],slope[1]); 755 #endif 756 xValue=bound[0]; 757 yValue0=(rowUpper0xValue*alpha[0])/element0; 758 yValue1=(rowUpper1xValue*alpha[1])/element1; 759 #ifdef PRINT_VALUES 760 printf("value of %d is %g, value of %d is max(%g,%g)  %g\n", 761 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1)); 762 #endif 763 newLower=CoinMin(newLower,CoinMax(yValue0,yValue1)); 764 // cost>0 so will be at lower 765 //double yValueAtBound0=newLower; 766 newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1)); 767 xValue=bound[1]; 768 yValue0=(rowUpper0xValue*alpha[0])/element0; 769 yValue1=(rowUpper1xValue*alpha[1])/element1; 770 #ifdef PRINT_VALUES 771 printf("value of %d is %g, value of %d is max(%g,%g)  %g\n", 772 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1)); 773 #endif 774 newLower=CoinMin(newLower,CoinMax(yValue0,yValue1)); 775 // cost>0 so will be at lower 776 //double yValueAtBound1=newLower; 777 newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1)); 778 lowerX=CoinMax(lowerX,newLower1.0e12*fabs(newLower)); 779 upperX=CoinMin(upperX,newUpper+1.0e12*fabs(newUpper)); 780 // Now make duplicate row 781 // keep row 0 so need to adjust costs so same 782 #ifdef PRINT_VALUES 783 printf("Costs for x %g,%g,%g are %g,%g,%g\n", 784 xValueEqual1.0,xValueEqual,xValueEqual+1.0, 785 costEqualslope[0],costEqual,costEqual+slope[1]); 786 #endif 787 double costOther=cost[otherCol]+slope[1]; 788 double costThis=cost[icol]+slope[1]*(element0/alpha[0]); 789 xValue=xValueEqual; 790 yValue0=CoinMax((rowUpper0xValue*alpha[0])/element0,lowerX); 791 double thisOffset=costEqual(costOther*xValue+costThis*yValue0); 792 offset += thisOffset; 793 #ifdef PRINT_VALUES 794 printf("new cost at equal %g\n",costOther*xValue+costThis*yValue0+thisOffset); 795 #endif 796 xValue=xValueEqual1.0; 797 yValue0=CoinMax((rowUpper0xValue*alpha[0])/element0,lowerX); 798 #ifdef PRINT_VALUES 799 printf("new cost at 1 %g\n",costOther*xValue+costThis*yValue0+thisOffset); 800 #endif 801 assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)(costEqualslope[0]))<1.0e5); 802 xValue=xValueEqual+1.0; 803 yValue0=CoinMax((rowUpper0xValue*alpha[0])/element0,lowerX); 804 #ifdef PRINT_VALUES 805 printf("new cost at +1 %g\n",costOther*xValue+costThis*yValue0+thisOffset); 806 #endif 807 assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)(costEqual+slope[1]))<1.0e5); 808 numberChanged++; 809 // continue; 810 cost[otherCol] = costOther; 811 cost[icol] = costThis; 812 clo[icol]=lowerX; 813 cup[icol]=upperX; 814 int startCol[2]; 815 int endCol[2]; 816 startCol[0]=mcstrt[icol]; 817 endCol[0]=startCol[0]+2; 818 startCol[1]=mcstrt[otherCol]; 819 endCol[1]=startCol[1]+hincol[otherCol]; 820 double values[2]={0.0,0.0}; 821 for (int k=0;k<2;k++) { 822 for (CoinBigIndex i=startCol[k];i<endCol[k];i++) { 823 if (hrow[i]==row0) 824 values[k]=colels[i]; 825 } 826 for (CoinBigIndex i=startCol[k];i<endCol[k];i++) { 827 if (hrow[i]==row1) 828 colels[i]=values[k]; 829 } 830 } 831 for (CoinBigIndex i=mrstrt[row1];i<mrstrt[row1]+2;i++) { 832 if (hcol[i]==icol) 833 rowels[i]=values[0]; 834 else 835 rowels[i]=values[1]; 836 } 837 } 838 } 839 } 840 } 841 #if ABC_NORMAL_DEBUG>0 842 if (offset) 843 printf("Cost offset %g\n",offset); 844 #endif 845 return numberChanged; 846 } 446 847 //#define COIN_PRESOLVE_BUG 447 848 #ifdef COIN_PRESOLVE_BUG … … 557 958 } 558 959 960 // transfer costs (may want to do it in OsiPresolve) 961 // need a transfer back at end of postsolve transferCosts(prob); 559 962 560 963 int iLoop; … … 569 972 paction_ = dupcol_action::presolve(prob, paction_); 570 973 } 974 #ifdef ABC_INHERIT 975 if (doTwoxTwo()) { 976 possibleSkip; 977 paction_ = twoxtwo_action::presolve(prob, paction_); 978 } 979 #endif 571 980 if (duprow) { 572 981 possibleSkip; 573 paction_ = duprow_action::presolve(prob, paction_); 982 int nTightened=tightenDoubletons2(prob); 983 if (nTightened) 984 printf("%d doubletons tightened\n",nTightened); 985 paction_ = duprow_action::presolve(prob, paction_); 574 986 } 575 987 if (doGubrow()) { … … 583 995 int lastDropped = 0; 584 996 prob>pass_ = 0; 997 #ifdef ABC_INHERIT 998 int numberRowsStart=nrows_prob>countEmptyRows(); 999 int numberColumnsStart=ncols_prob>countEmptyCols(); 1000 int numberRowsLeft=numberRowsStart; 1001 int numberColumnsLeft=numberColumnsStart; 1002 bool lastPassWasGood=true; 1003 #if ABC_NORMAL_DEBUG 1004 printf("Original rows,columns %d,%d starting first pass with %d,%d\n", 1005 nrows_,ncols_,numberRowsLeft,numberColumnsLeft); 1006 #endif 1007 #endif 585 1008 if (numberPasses_<=5) 586 1009 prob>presolveOptions_ = 0x10000; // say more lightweight … … 867 1290 if (paction_ == paction0  stopLoop) 868 1291 break; 869 1292 #ifdef ABC_INHERIT 1293 // see whether to stop anyway 1294 int numberRowsNow=nrows_prob>countEmptyRows(); 1295 int numberColumnsNow=ncols_prob>countEmptyCols(); 1296 #if ABC_NORMAL_DEBUG 1297 printf("Original rows,columns %d,%d  last %d,%d end of pass %d has %d,%d\n", 1298 nrows_,ncols_,numberRowsLeft,numberColumnsLeft,iLoop+1,numberRowsNow, 1299 numberColumnsNow); 1300 #endif 1301 int rowsDeleted=numberRowsLeftnumberRowsNow; 1302 int columnsDeleted=numberColumnsLeftnumberColumnsNow; 1303 if (iLoop>15) { 1304 if (rowsDeleted*100<numberRowsStart&& 1305 columnsDeleted*100<numberColumnsStart) 1306 break; 1307 lastPassWasGood=true; 1308 } else if (rowsDeleted*100<numberRowsStart&&rowsDeleted<500&& 1309 columnsDeleted*100<numberColumnsStart&&columnsDeleted<500) { 1310 if (!lastPassWasGood) 1311 break; 1312 else 1313 lastPassWasGood=false; 1314 } else { 1315 lastPassWasGood=true; 1316 } 1317 numberRowsLeft=numberRowsNow; 1318 numberColumnsLeft=numberColumnsNow; 1319 #endif 870 1320 } 871 1321 } … … 943 1393 } 944 1394 } 1395 } 1396 if (prob.maxmin_<0) { 1397 //for (int i=0;i<presolvedModel_>numberRows();i++) 1398 //prob.rowduals_[i]=prob.rowduals_[i]; 1399 for (int i=0;i<ncols_;i++) { 1400 prob.cost_[i]=prob.cost_[i]; 1401 //prob.rcosts_[i]=prob.rcosts_[i]; 1402 } 1403 prob.maxmin_=1.0; 945 1404 } 946 1405 const CoinPresolveAction *paction = paction_; … … 1245 1704 mcstrt_[0] = 0; 1246 1705 ClpDisjointCopyN(m>getVectorLengths(), ncols_, hincol_); 1706 if (si>getObjSense() < 0.0) { 1707 for (int i=0;i<ncols_;i++) 1708 cost_[i]=cost_[i]; 1709 maxmin_=1.0; 1710 } 1247 1711 for (icol = 0; icol < ncols_; icol++) { 1248 1712 CoinBigIndex j; … … 1416 1880 #endif 1417 1881 { 1882 if (si>getObjSense() < 0.0) { 1883 for (int i=0;i<ncols_;i++) 1884 cost_[i]=cost_[i]; 1885 dobias_=dobias_; 1886 } 1418 1887 si>loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_, 1419 1888 clo_, cup_, cost_, rlo_, rup_); … … 1436 1905 #endif 1437 1906 si>setDblParam(ClpObjOffset, originalOffset_  dobias_); 1907 if (si>getObjSense() < 0.0) { 1908 // put back 1909 for (int i=0;i<ncols_;i++) 1910 cost_[i]=cost_[i]; 1911 dobias_=dobias_; 1912 maxmin_=1.0; 1913 } 1438 1914 1439 1915 } … … 1904 2380 if (rowObjective_) { 1905 2381 int iRow; 2382 #ifndef NDEBUG 1906 2383 int k = 1; 2384 #endif 1907 2385 int nObj = 0; 1908 2386 for (iRow = 0; iRow < nrowsNow; iRow++) { 1909 2387 int kRow = originalRow_[iRow]; 2388 #ifndef NDEBUG 1910 2389 assert (kRow > k); 1911 2390 k = kRow; 2391 #endif 1912 2392 rowObjective_[iRow] = rowObjective_[kRow]; 1913 2393 if (rowObjective_[iRow])
Note: See TracChangeset
for help on using the changeset viewer.