Changeset 530 for branches/devel
- Timestamp:
- Feb 8, 2007 12:45:15 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/devel/Cbc/src/CoinSolve.cpp
r525 r530 38 38 #include "OsiRowCutDebugger.hpp" 39 39 #include "OsiChooseVariable.hpp" 40 //#define CLP_MALLOC_STATISTICS 41 #ifdef CLP_MALLOC_STATISTICS 42 #include <exception> 43 #include <new> 44 static double malloc_times=0.0; 45 static double malloc_total=0.0; 46 static int malloc_amount[]={0,32,128,256,1024,4096,16384,65536,262144,INT_MAX}; 47 static int malloc_n=10; 48 double malloc_counts[10]={0,0,0,0,0,0,0,0,0,0}; 49 void * operator new (size_t size) 50 { 51 malloc_times ++; 52 malloc_total += size; 53 int i; 54 for (i=0;i<malloc_n;i++) { 55 if ((int) size<=malloc_amount[i]) { 56 malloc_counts[i]++; 57 break; 58 } 59 } 60 void * p =malloc(size); 61 return p; 62 } 63 void operator delete (void *p) 64 { 65 free(p); 66 } 67 static voif malloc_stats() 68 { 69 double average = malloc_total/malloc_times; 70 printf("count %g bytes %g - average %g\n",malloc_times,malloc_total,average); 71 for (int i=0;i<malloc_n;i++) 72 printf("%g ",malloc_counts[i]); 73 printf("\n"); 74 } 75 #endif 40 76 #ifdef DMALLOC 41 77 #include "dmalloc.h" … … 56 92 #include "OsiRowCut.hpp" 57 93 #include "OsiColCut.hpp" 94 #ifndef COIN_HAS_LINK 58 95 #ifdef COIN_HAS_ASL 59 96 #define COIN_HAS_LINK 60 97 #endif 61 #ifndef COIN_DEVELOP62 #undef COIN_HAS_LINK63 98 #endif 64 99 #ifdef COIN_HAS_LINK … … 426 461 } 427 462 } 463 /* Returns OsiSolverInterface (User should delete) 464 On entry numberKnapsack is maximum number of Total entries 465 */ 466 static OsiSolverInterface * 467 expandKnapsack(CoinModel & model, int * whichColumn, int * knapsackStart, 468 int * knapsackRow, int &numberKnapsack, 469 CglStored & stored, int logLevel, 470 int fixedPriority) 471 { 472 int maxTotal = numberKnapsack; 473 // load from coin model 474 OsiSolverLink *si = new OsiSolverLink(); 475 OsiSolverInterface * finalModel=NULL; 476 si->setDefaultMeshSize(0.001); 477 // need some relative granularity 478 si->setDefaultBound(100.0); 479 si->setDefaultMeshSize(0.01); 480 si->setDefaultBound(1000.0); 481 si->setIntegerPriority(1000); 482 si->setBiLinearPriority(10000); 483 si->load(model,true,logLevel); 484 // get priorities 485 const int * priorities=model.priorities(); 486 int numberColumns = model.numberColumns(); 487 int SOSPriority=10000; 488 if (priorities) { 489 OsiObject ** objects = si->objects(); 490 int numberObjects = si->numberObjects(); 491 for (int iObj = 0;iObj<numberObjects;iObj++) { 492 int iColumn = objects[iObj]->columnNumber(); 493 if (iColumn>=0&&iColumn<numberColumns) { 494 OsiSimpleInteger * obj = 495 dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ; 496 assert (obj); 497 int iPriority = priorities[iColumn]; 498 if (iPriority>0) 499 objects[iObj]->setPriority(iPriority); 500 } 501 } 502 if (fixedPriority>0) { 503 si->setFixedPriority(fixedPriority); 504 SOSPriority=fixedPriority+1; 505 } 506 } 507 CoinModel coinModel=*si->coinModel(); 508 assert(coinModel.numberRows()>0); 509 int numberRows = coinModel.numberRows(); 510 // Mark variables 511 int * whichKnapsack = new int [numberColumns]; 512 int iRow,iColumn; 513 for (iColumn=0;iColumn<numberColumns;iColumn++) 514 whichKnapsack[iColumn]=-1; 515 int kRow; 516 bool badModel=false; 517 // analyze 518 if (logLevel>1) { 519 for (iRow=0;iRow<numberRows;iRow++) { 520 /* Just obvious one at first 521 positive non unit coefficients 522 all integer 523 positive rowUpper 524 for now - linear (but further down in code may use nonlinear) 525 column bounds should be tight 526 */ 527 //double lower = coinModel.getRowLower(iRow); 528 double upper = coinModel.getRowUpper(iRow); 529 if (upper<1.0e10) { 530 CoinModelLink triple=coinModel.firstInRow(iRow); 531 bool possible=true; 532 int n=0; 533 int n1=0; 534 while (triple.column()>=0) { 535 int iColumn = triple.column(); 536 const char * el = coinModel.getElementAsString(iRow,iColumn); 537 if (!strcmp("Numeric",el)) { 538 if (coinModel.columnLower(iColumn)==coinModel.columnUpper(iColumn)) { 539 triple=coinModel.next(triple); 540 continue; // fixed 541 } 542 double value=coinModel.getElement(iRow,iColumn); 543 if (value<0.0) { 544 possible=false; 545 } else { 546 n++; 547 if (value==1.0) 548 n1++; 549 if (coinModel.columnLower(iColumn)<0.0) 550 possible=false; 551 if (!coinModel.isInteger(iColumn)) 552 possible=false; 553 if (whichKnapsack[iColumn]>=0) 554 possible=false; 555 } 556 } else { 557 possible=false; // non linear 558 } 559 triple=coinModel.next(triple); 560 } 561 if (n-n1>1&&possible) { 562 double lower = coinModel.getRowLower(iRow); 563 double upper = coinModel.getRowUpper(iRow); 564 CoinModelLink triple=coinModel.firstInRow(iRow); 565 while (triple.column()>=0) { 566 int iColumn = triple.column(); 567 lower -= coinModel.columnLower(iColumn)*triple.value(); 568 upper -= coinModel.columnLower(iColumn)*triple.value(); 569 triple=coinModel.next(triple); 570 } 571 printf("%d is possible %g <=",iRow,lower); 572 // print 573 triple=coinModel.firstInRow(iRow); 574 while (triple.column()>=0) { 575 int iColumn = triple.column(); 576 if (coinModel.columnLower(iColumn)!=coinModel.columnUpper(iColumn)) 577 printf(" (%d,el %g up %g)",iColumn,triple.value(), 578 coinModel.columnUpper(iColumn)-coinModel.columnLower(iColumn)); 579 triple=coinModel.next(triple); 580 } 581 printf(" <= %g\n",upper); 582 } 583 } 584 } 585 } 586 numberKnapsack=0; 587 for (kRow=0;kRow<numberRows;kRow++) { 588 iRow=kRow; 589 /* Just obvious one at first 590 positive non unit coefficients 591 all integer 592 positive rowUpper 593 for now - linear (but further down in code may use nonlinear) 594 column bounds should be tight 595 */ 596 //double lower = coinModel.getRowLower(iRow); 597 double upper = coinModel.getRowUpper(iRow); 598 if (upper<1.0e10) { 599 CoinModelLink triple=coinModel.firstInRow(iRow); 600 bool possible=true; 601 int n=0; 602 int n1=0; 603 while (triple.column()>=0) { 604 int iColumn = triple.column(); 605 const char * el = coinModel.getElementAsString(iRow,iColumn); 606 if (!strcmp("Numeric",el)) { 607 if (coinModel.columnLower(iColumn)==coinModel.columnUpper(iColumn)) { 608 triple=coinModel.next(triple); 609 continue; // fixed 610 } 611 double value=coinModel.getElement(iRow,iColumn); 612 if (value<0.0) { 613 possible=false; 614 } else { 615 n++; 616 if (value==1.0) 617 n1++; 618 if (coinModel.columnLower(iColumn)<0.0) 619 possible=false; 620 if (!coinModel.isInteger(iColumn)) 621 possible=false; 622 if (whichKnapsack[iColumn]>=0) 623 possible=false; 624 } 625 } else { 626 possible=false; // non linear 627 } 628 triple=coinModel.next(triple); 629 } 630 if (n-n1>1&&possible) { 631 // try 632 CoinModelLink triple=coinModel.firstInRow(iRow); 633 while (triple.column()>=0) { 634 int iColumn = triple.column(); 635 if (coinModel.columnLower(iColumn)!=coinModel.columnUpper(iColumn)) 636 whichKnapsack[iColumn]=numberKnapsack; 637 triple=coinModel.next(triple); 638 } 639 knapsackRow[numberKnapsack++]=iRow; 640 } 641 } 642 } 643 if (logLevel>0) 644 printf("%d out of %d candidate rows are possible\n",numberKnapsack,numberRows); 645 // Check whether we can get rid of nonlinearities 646 /* mark rows 647 -2 in knapsack and other variables 648 -1 not involved 649 n only in knapsack n 650 */ 651 int * markRow = new int [numberRows]; 652 for (iRow=0;iRow<numberRows;iRow++) 653 markRow[iRow]=-1; 654 int canDo=1; // OK and linear 655 for (iColumn=0;iColumn<numberColumns;iColumn++) { 656 CoinModelLink triple=coinModel.firstInColumn(iColumn); 657 int iKnapsack = whichKnapsack[iColumn]; 658 bool linear=true; 659 // See if quadratic objective 660 const char * expr = coinModel.getColumnObjectiveAsString(iColumn); 661 if (strcmp(expr,"Numeric")) { 662 linear=false; 663 } 664 while (triple.row()>=0) { 665 int iRow = triple.row(); 666 if (iKnapsack>=0) { 667 if (markRow[iRow]==-1) { 668 markRow[iRow]=iKnapsack; 669 } else if (markRow[iRow]!=iKnapsack) { 670 markRow[iRow]=-2; 671 } 672 } 673 const char * expr = coinModel.getElementAsString(iRow,iColumn); 674 if (strcmp(expr,"Numeric")) { 675 linear=false; 676 } 677 triple=coinModel.next(triple); 678 } 679 if (!linear) { 680 if (whichKnapsack[iColumn]<0) { 681 canDo=0; 682 break; 683 } else { 684 canDo=2; 685 } 686 } 687 } 688 int * markKnapsack = NULL; 689 double * coefficient = NULL; 690 double * linear = NULL; 691 int * whichRow = NULL; 692 int * lookupRow = NULL; 693 badModel=(canDo==0); 694 if (numberKnapsack&&canDo) { 695 /* double check - OK if 696 no nonlinear 697 nonlinear only on columns in knapsack 698 nonlinear only on columns in knapsack * ONE other - same for all in knapsack 699 AND that is only row connected to knapsack 700 (theoretically could split knapsack if two other and small numbers) 701 also ONE could be ONE expression - not just a variable 702 */ 703 int iKnapsack; 704 markKnapsack = new int [numberKnapsack]; 705 coefficient = new double [numberKnapsack]; 706 linear = new double [numberColumns]; 707 for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) 708 markKnapsack[iKnapsack]=-1; 709 if (canDo==2) { 710 for (iRow=-1;iRow<numberRows;iRow++) { 711 int numberOdd; 712 CoinPackedMatrix * row = coinModel.quadraticRow(iRow,linear,numberOdd); 713 if (row) { 714 // see if valid 715 const double * element = row->getElements(); 716 const int * column = row->getIndices(); 717 const CoinBigIndex * columnStart = row->getVectorStarts(); 718 const int * columnLength = row->getVectorLengths(); 719 int numberLook = row->getNumCols(); 720 for (int i=0;i<numberLook;i++) { 721 int iKnapsack=whichKnapsack[i]; 722 if (iKnapsack<0) { 723 // might be able to swap - but for now can't have knapsack in 724 for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) { 725 int iColumn = column[j]; 726 if (whichKnapsack[iColumn]>=0) { 727 canDo=0; // no good 728 badModel=true; 729 break; 730 } 731 } 732 } else { 733 // OK if in same knapsack - or maybe just one 734 int marked=markKnapsack[iKnapsack]; 735 for (int j=columnStart[i];j<columnStart[i]+columnLength[i];j++) { 736 int iColumn = column[j]; 737 if (whichKnapsack[iColumn]!=iKnapsack&&whichKnapsack[iColumn]>=0) { 738 canDo=0; // no good 739 badModel=true; 740 break; 741 } else if (marked==-1) { 742 markKnapsack[iKnapsack]=iColumn; 743 marked=iColumn; 744 coefficient[iKnapsack]=element[j]; 745 coinModel.associateElement(coinModel.columnName(iColumn),1.0); 746 } else if (marked!=iColumn) { 747 badModel=true; 748 canDo=0; // no good 749 break; 750 } else { 751 // could manage with different coefficients - but for now ... 752 assert(coefficient[iKnapsack]==element[j]); 753 } 754 } 755 } 756 } 757 delete row; 758 } 759 } 760 } 761 if (canDo) { 762 // for any rows which are cuts 763 whichRow = new int [numberRows]; 764 lookupRow = new int [numberRows]; 765 bool someNonlinear=false; 766 double maxCoefficient=1.0; 767 for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) { 768 if (markKnapsack[iKnapsack]>=0) { 769 someNonlinear=true; 770 int iColumn = markKnapsack[iKnapsack]; 771 maxCoefficient = CoinMax(maxCoefficient,fabs(coefficient[iKnapsack]*coinModel.columnUpper(iColumn))); 772 } 773 } 774 if (someNonlinear) { 775 // associate all columns to stop possible error messages 776 for (iColumn=0;iColumn<numberColumns;iColumn++) { 777 coinModel.associateElement(coinModel.columnName(iColumn),1.0); 778 } 779 } 780 ClpSimplex tempModel; 781 tempModel.loadProblem(coinModel); 782 // Create final model - first without knapsacks 783 int nCol=0; 784 int nRow=0; 785 for (iRow=0;iRow<numberRows;iRow++) { 786 if (markRow[iRow]<0) { 787 lookupRow[iRow]=nRow; 788 whichRow[nRow++]=iRow; 789 } else { 790 lookupRow[iRow]=-1; 791 } 792 } 793 for (iColumn=0;iColumn<numberColumns;iColumn++) { 794 if (whichKnapsack[iColumn]<0) 795 whichColumn[nCol++]=iColumn; 796 } 797 ClpSimplex finalModelX(&tempModel,nRow,whichRow,nCol,whichColumn,false,false,false); 798 OsiClpSolverInterface finalModelY(&finalModelX,true); 799 finalModel = finalModelY.clone(); 800 finalModelY.releaseClp(); 801 // Put back priorities 802 const int * priorities=model.priorities(); 803 if (priorities) { 804 finalModel->findIntegers(false); 805 OsiObject ** objects = finalModel->objects(); 806 int numberObjects = finalModel->numberObjects(); 807 for (int iObj = 0;iObj<numberObjects;iObj++) { 808 int iColumn = objects[iObj]->columnNumber(); 809 if (iColumn>=0&&iColumn<nCol) { 810 OsiSimpleInteger * obj = 811 dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ; 812 assert (obj); 813 int iPriority = priorities[whichColumn[iColumn]]; 814 if (iPriority>0) 815 objects[iObj]->setPriority(iPriority); 816 } 817 } 818 } 819 for (iRow=0;iRow<numberRows;iRow++) { 820 whichRow[iRow]=iRow; 821 } 822 int numberOther=finalModel->getNumCols(); 823 int nLargest=0; 824 int nelLargest=0; 825 int nTotal=0; 826 for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) { 827 iRow = knapsackRow[iKnapsack]; 828 int nCreate = maxTotal; 829 int nelCreate=coinModel.expandKnapsack(iRow,nCreate,NULL,NULL,NULL,NULL); 830 if (nelCreate<0) 831 badModel=true; 832 nTotal+=nCreate; 833 nLargest = CoinMax(nLargest,nCreate); 834 nelLargest = CoinMax(nelLargest,nelCreate); 835 } 836 if (nTotal>maxTotal) 837 badModel=true; 838 if (!badModel) { 839 // Now arrays for building 840 nelLargest = CoinMax(nelLargest,nLargest)+1; 841 double * buildObj = new double [nLargest]; 842 double * buildElement = new double [nelLargest]; 843 int * buildStart = new int[nLargest+1]; 844 int * buildRow = new int[nelLargest]; 845 OsiObject ** object = new OsiObject * [numberKnapsack]; 846 int nSOS=0; 847 for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) { 848 knapsackStart[iKnapsack]=finalModel->getNumCols(); 849 iRow = knapsackRow[iKnapsack]; 850 int nCreate = 10000; 851 coinModel.expandKnapsack(iRow,nCreate,buildObj,buildStart,buildRow,buildElement); 852 // Redo row numbers 853 for (iColumn=0;iColumn<nCreate;iColumn++) { 854 for (int j=buildStart[iColumn];j<buildStart[iColumn+1];j++) { 855 int jRow=buildRow[j]; 856 jRow=lookupRow[jRow]; 857 assert (jRow>=0&&jRow<nRow); 858 buildRow[j]=jRow; 859 } 860 } 861 finalModel->addCols(nCreate,buildStart,buildRow,buildElement,NULL,NULL,buildObj); 862 int numberFinal=finalModel->getNumCols(); 863 for (iColumn=numberOther;iColumn<numberFinal;iColumn++) { 864 if (markKnapsack[iKnapsack]<0) { 865 finalModel->setColUpper(iColumn,maxCoefficient); 866 finalModel->setInteger(iColumn); 867 } else { 868 finalModel->setColUpper(iColumn,maxCoefficient+1.0); 869 finalModel->setInteger(iColumn); 870 } 871 buildRow[iColumn-numberOther]=iColumn; 872 buildElement[iColumn-numberOther]=1.0; 873 } 874 if (markKnapsack[iKnapsack]<0) { 875 // convexity row 876 finalModel->addRow(numberFinal-numberOther,buildRow,buildElement,1.0,1.0); 877 } else { 878 int iColumn = markKnapsack[iKnapsack]; 879 int n=numberFinal-numberOther; 880 buildRow[n]=iColumn; 881 buildElement[n++]=coefficient[iKnapsack]; 882 // convexity row (sort of) 883 finalModel->addRow(n,buildRow,buildElement,0.0,0.0); 884 OsiSOS * sosObject = new OsiSOS(finalModel,n-1,buildRow,NULL,1); 885 sosObject->setPriority(iKnapsack+SOSPriority); 886 // Say not integral even if is (switch off heuristics) 887 sosObject->setIntegerValued(false); 888 object[nSOS++]=sosObject; 889 } 890 numberOther=numberFinal; 891 } 892 finalModel->addObjects(nSOS,object); 893 for (iKnapsack=0;iKnapsack<nSOS;iKnapsack++) 894 delete object[iKnapsack]; 895 delete [] object; 896 // Can we move any rows to cuts 897 const int * cutMarker = coinModel.cutMarker(); 898 if (cutMarker) { 899 // Row copy 900 const CoinPackedMatrix * matrixByRow = finalModel->getMatrixByRow(); 901 const double * elementByRow = matrixByRow->getElements(); 902 const int * column = matrixByRow->getIndices(); 903 const CoinBigIndex * rowStart = matrixByRow->getVectorStarts(); 904 const int * rowLength = matrixByRow->getVectorLengths(); 905 906 const double * rowLower = finalModel->getRowLower(); 907 const double * rowUpper = finalModel->getRowUpper(); 908 int nDelete=0; 909 for (iRow=0;iRow<numberRows;iRow++) { 910 if (cutMarker[iRow]&&lookupRow[iRow]>=0) { 911 int jRow=lookupRow[iRow]; 912 whichRow[nDelete++]=jRow; 913 int start = rowStart[jRow]; 914 stored.addCut(rowLower[jRow],rowUpper[jRow], 915 rowLength[jRow],column+start,elementByRow+start); 916 } 917 } 918 finalModel->deleteRows(nDelete,whichRow); 919 } 920 knapsackStart[numberKnapsack]=finalModel->getNumCols(); 921 delete [] buildObj; 922 delete [] buildElement; 923 delete [] buildStart; 924 delete [] buildRow; 925 finalModel->writeMps("full"); 926 } 927 } 928 } 929 delete [] whichKnapsack; 930 delete [] markRow; 931 delete [] markKnapsack; 932 delete [] coefficient; 933 delete [] linear; 934 delete [] whichRow; 935 delete [] lookupRow; 936 delete si; 937 si=NULL; 938 if (!badModel) { 939 return finalModel; 940 } else { 941 delete finalModel; 942 return NULL; 943 } 944 } 945 // Fills in original solution (coinModel length) 946 static void 947 afterKnapsack(const CoinModel & coinModel2, const int * whichColumn, const int * knapsackStart, 948 const int * knapsackRow, int numberKnapsack, 949 const double * knapsackSolution, double * solution, int logLevel) 950 { 951 CoinModel coinModel = coinModel2; 952 int numberColumns = coinModel.numberColumns(); 953 int iColumn; 954 // associate all columns to stop possible error messages 955 for (iColumn=0;iColumn<numberColumns;iColumn++) { 956 coinModel.associateElement(coinModel.columnName(iColumn),1.0); 957 } 958 CoinZeroN(solution,numberColumns); 959 int nCol=knapsackStart[0]; 960 for (iColumn=0;iColumn<nCol;iColumn++) { 961 int jColumn = whichColumn[iColumn]; 962 solution[jColumn]=knapsackSolution[iColumn]; 963 } 964 int * buildRow = new int [numberColumns]; // wild overkill 965 double * buildElement = new double [numberColumns]; 966 int iKnapsack; 967 for (iKnapsack=0;iKnapsack<numberKnapsack;iKnapsack++) { 968 int k=-1; 969 double value=0.0; 970 for (iColumn=knapsackStart[iKnapsack];iColumn<knapsackStart[iKnapsack+1];iColumn++) { 971 if (knapsackSolution[iColumn]>1.0e-5) { 972 if (k>=0) { 973 printf("Two nonzero values for knapsack %d at (%d,%g) and (%d,%g)\n",iKnapsack, 974 k,knapsackSolution[k],iColumn,knapsackSolution[iColumn]); 975 abort(); 976 } 977 k=iColumn; 978 value=floor(knapsackSolution[iColumn]+0.5); 979 assert (fabs(value-knapsackSolution[iColumn])<1.0e-5); 980 } 981 } 982 if (k>=0) { 983 int iRow = knapsackRow[iKnapsack]; 984 int nCreate = 10000; 985 int nel=coinModel.expandKnapsack(iRow,nCreate,NULL,NULL,buildRow,buildElement,k-knapsackStart[iKnapsack]); 986 assert (nel); 987 if (logLevel>0) 988 printf("expanded column %d in knapsack %d has %d nonzero entries:\n", 989 k-knapsackStart[iKnapsack],iKnapsack,nel); 990 for (int i=0;i<nel;i++) { 991 int jColumn = buildRow[i]; 992 double value = buildElement[i]; 993 if (logLevel>0) 994 printf("%d - original %d has value %g\n",i,jColumn,value); 995 solution[jColumn]=value; 996 } 997 } 998 } 999 delete [] buildRow; 1000 delete [] buildElement; 1001 #if 0 1002 for (iColumn=0;iColumn<numberColumns;iColumn++) { 1003 if (solution[iColumn]>1.0e-5&&coinModel.isInteger(iColumn)) 1004 printf("%d %g\n",iColumn,solution[iColumn]); 1005 } 1006 #endif 1007 } 428 1008 static int outDupRow(OsiSolverInterface * solver) 429 1009 { … … 606 1186 char * sosType = NULL; 607 1187 double * sosReference = NULL; 1188 int * cut=NULL; 608 1189 int * sosPriority=NULL; 609 1190 #ifdef COIN_HAS_ASL 610 1191 ampl_info info; 1192 CglStored storedAmpl; 611 1193 CoinModel * coinModel = NULL; 1194 CoinModel saveCoinModel; 1195 int * whichColumn = NULL; 1196 int * knapsackStart=NULL; 1197 int * knapsackRow=NULL; 1198 int numberKnapsack=0; 612 1199 memset(&info,0,sizeof(info)); 613 1200 if (argc>2&&!strcmp(argv[2],"-AMPL")) { … … 653 1240 info.columnLower,info.columnUpper,info.objective, 654 1241 info.rowLower,info.rowUpper); 1242 // take off cuts if ampl wants that 1243 if (info.cut) { 1244 int numberRows = info.numberRows; 1245 int * whichRow = new int [numberRows]; 1246 // Row copy 1247 const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow(); 1248 const double * elementByRow = matrixByRow->getElements(); 1249 const int * column = matrixByRow->getIndices(); 1250 const CoinBigIndex * rowStart = matrixByRow->getVectorStarts(); 1251 const int * rowLength = matrixByRow->getVectorLengths(); 1252 1253 const double * rowLower = solver->getRowLower(); 1254 const double * rowUpper = solver->getRowUpper(); 1255 int nDelete=0; 1256 for (int iRow=0;iRow<numberRows;iRow++) { 1257 if (info.cut[iRow]) { 1258 whichRow[nDelete++]=iRow; 1259 int start = rowStart[iRow]; 1260 storedAmpl.addCut(rowLower[iRow],rowUpper[iRow], 1261 rowLength[iRow],column+start,elementByRow+start); 1262 } 1263 } 1264 solver->deleteRows(nDelete,whichRow); 1265 delete [] whichRow; 1266 } 655 1267 #ifdef COIN_HAS_LINK 656 1268 } else { 1269 // save 1270 saveCoinModel = *coinModel; 657 1271 // load from coin model 658 1272 OsiSolverLink solver1; … … 1820 2434 } 1821 2435 if (!miplib) { 2436 if (!preSolve) { 2437 model.solver()->setHintParam(OsiDoPresolveInInitial,false,OsiHintTry); 2438 model.solver()->setHintParam(OsiDoPresolveInResolve,false,OsiHintTry); 2439 } 1822 2440 model.initialSolve(); 1823 2441 OsiSolverInterface * solver = model.solver(); … … 1970 2588 preProcess=0; 1971 2589 useStrategy=true; 2590 // empty out any cuts 2591 if (storedAmpl.sizeRowCuts()) { 2592 printf("Emptying ampl stored cuts as internal preprocessing\n"); 2593 CglStored temp; 2594 storedAmpl=temp; 2595 } 1972 2596 } 1973 2597 if (preProcess&&type==BAB) { … … 2127 2751 } 2128 2752 int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue(); 2753 #ifdef COIN_HAS_LINK 2754 // If linked then see if expansion wanted 2755 { 2756 OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver()); 2757 if (solver3) { 2758 int options = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()/10000; 2759 if (options) { 2760 /* 2761 1 - force mini branch and bound 2762 2 - set priorities high on continuous 2763 4 - try adding OA cuts 2764 8 - try doing quadratic linearization 2765 16 - try expanding knapsacks 2766 */ 2767 if ((options&16)) { 2768 int numberColumns = saveCoinModel.numberColumns(); 2769 int numberRows = saveCoinModel.numberRows(); 2770 whichColumn = new int[numberColumns]; 2771 knapsackStart=new int[numberRows+1]; 2772 knapsackRow=new int[numberRows]; 2773 numberKnapsack=10000; 2774 OsiSolverInterface * solver = expandKnapsack(saveCoinModel,whichColumn,knapsackStart, 2775 knapsackRow,numberKnapsack, 2776 storedAmpl,2,slpValue); 2777 if (solver) { 2778 clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver); 2779 assert (clpSolver); 2780 lpSolver = clpSolver->getModelPtr(); 2781 babModel->assignSolver(solver); 2782 testOsiOptions=0; 2783 } else { 2784 numberKnapsack=0; 2785 delete [] whichColumn; 2786 delete [] knapsackStart; 2787 delete [] knapsackRow; 2788 whichColumn=NULL; 2789 knapsackStart=NULL; 2790 knapsackRow=NULL; 2791 } 2792 } 2793 } 2794 } 2795 } 2796 #endif 2129 2797 if (useCosts&&testOsiOptions<0) { 2130 2798 int numberColumns = babModel->getNumCols(); … … 2835 3503 free(sosReference); 2836 3504 sosReference=NULL; 3505 free(cut); 3506 cut=NULL; 2837 3507 free(sosPriority); 2838 3508 sosPriority=NULL; … … 2918 3588 4 - try adding OA cuts 2919 3589 8 - try doing quadratic linearization 3590 16 - try expanding knapsacks 2920 3591 */ 2921 3592 if ((options&2)) … … 3011 3682 osiclp->setSolveOptions(clpSolve); 3012 3683 osiclp->setHintParam(OsiDoDualInResolve,false); 3684 // switch off row copy 3685 osiclp->getModelPtr()->setSpecialOptions(osiclp->getModelPtr()->specialOptions()|256); 3686 } 3687 if (storedAmpl.sizeRowCuts()) { 3688 if (preProcess) { 3689 const int * originalColumns = process.originalColumns(); 3690 int numberColumns = babModel->getNumCols(); 3691 int * newColumn = new int[numberOriginalColumns]; 3692 int i; 3693 for (i=0;i<numberOriginalColumns;i++) 3694 newColumn[i]=-1; 3695 for (i=0;i<numberColumns;i++) { 3696 int iColumn = originalColumns[i]; 3697 newColumn[iColumn]=i; 3698 } 3699 int * buildColumn = new int[numberColumns]; 3700 // Build up valid cuts 3701 int nBad=0; 3702 int nCuts = storedAmpl.sizeRowCuts(); 3703 CglStored newCuts; 3704 for (i=0;i<nCuts;i++) { 3705 const OsiRowCut * cut = storedAmpl.rowCutPointer(i); 3706 double lb = cut->lb(); 3707 double ub = cut->ub(); 3708 int n=cut->row().getNumElements(); 3709 const int * column = cut->row().getIndices(); 3710 const double * element = cut->row().getElements(); 3711 bool bad=false; 3712 for (int i=0;i<n;i++) { 3713 int iColumn = column[i]; 3714 iColumn = newColumn[iColumn]; 3715 if (iColumn>=0) { 3716 buildColumn[i]=iColumn; 3717 } else { 3718 bad=true; 3719 break; 3720 } 3721 } 3722 if (!bad) { 3723 newCuts.addCut(lb,ub,n,buildColumn,element); 3724 } else { 3725 nBad++; 3726 } 3727 } 3728 storedAmpl=newCuts; 3729 if (nBad) 3730 printf("%d cuts dropped\n",nBad); 3731 delete [] newColumn; 3732 delete [] buildColumn; 3733 } 3734 if (storedAmpl.sizeRowCuts()) 3735 babModel->addCutGenerator(&storedAmpl,1,"AmplStored"); 3013 3736 } 3014 3737 babModel->branchAndBound(statistics); … … 3076 3799 if (type==STRENGTHEN&&strengthenedModel) 3077 3800 clpSolver = dynamic_cast< OsiClpSolverInterface*> (strengthenedModel); 3801 else if (usingAmpl) 3802 clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver()); 3078 3803 lpSolver = clpSolver->getModelPtr(); 3079 3804 if (numberChanged) { … … 3106 3831 #ifdef COIN_HAS_ASL 3107 3832 if (usingAmpl) { 3833 clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver()); 3834 lpSolver = clpSolver->getModelPtr(); 3108 3835 double value = babModel->getObjValue()*lpSolver->getObjSense(); 3109 3836 char buf[300]; … … 3148 3875 if (bestSolution) { 3149 3876 free(info.primalSolution); 3150 info.primalSolution = (double *) malloc(n*sizeof(double)); 3151 CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution); 3152 int numberRows = lpSolver->numberRows(); 3153 free(info.dualSolution); 3154 info.dualSolution = (double *) malloc(numberRows*sizeof(double)); 3155 CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution); 3877 if (!numberKnapsack) { 3878 info.primalSolution = (double *) malloc(n*sizeof(double)); 3879 CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution); 3880 int numberRows = lpSolver->numberRows(); 3881 free(info.dualSolution); 3882 info.dualSolution = (double *) malloc(numberRows*sizeof(double)); 3883 CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution); 3884 } else { 3885 // expanded knapsack 3886 info.dualSolution=NULL; 3887 int numberColumns = saveCoinModel.numberColumns(); 3888 info.primalSolution = (double *) malloc(numberColumns*sizeof(double)); 3889 // Fills in original solution (coinModel length) 3890 afterKnapsack(saveCoinModel, whichColumn, knapsackStart, 3891 knapsackRow, numberKnapsack, 3892 lpSolver->primalColumnSolution(), info.primalSolution,1); 3893 } 3156 3894 } else { 3157 3895 info.primalSolution=NULL; … … 3198 3936 free(sosReference); 3199 3937 sosReference=NULL; 3938 free(cut); 3939 cut=NULL; 3200 3940 free(sosPriority); 3201 3941 sosPriority=NULL;
Note: See TracChangeset
for help on using the changeset viewer.