#include <CglKnapsackCover.hpp>
Inheritance diagram for CglKnapsackCover:
Public Member Functions | |
void | setTestedRowIndices (int num, const int *ind) |
Generate Cuts | |
virtual void | generateCuts (const OsiSolverInterface &si, OsiCuts &cs) const |
Constructors and destructors | |
CglKnapsackCover () | |
Default constructor. | |
CglKnapsackCover (const CglKnapsackCover &) | |
Copy constructor. | |
CglKnapsackCover & | operator= (const CglKnapsackCover &rhs) |
Assignment operator. | |
virtual | ~CglKnapsackCover () |
Destructor. | |
Sets and gets | |
void | setMaxInKnapsack (int value) |
Set limit on number in knapsack. | |
int | getMaxInKnapsack () const |
get limit on number in knapsack | |
Private Member Functions | |
Private methods | |
int | deriveAKnapsack (const OsiSolverInterface &si, OsiCuts &cs, CoinPackedVector &krow, bool treatAsLRow, double &b, int *complement, double *xstar, int rowIndex, int numberElements, const int *index, const double *element) const |
int | deriveAKnapsack (const OsiSolverInterface &si, OsiCuts &cs, CoinPackedVector &krow, double &b, int *complement, double *xstar, int rowIndex, const CoinPackedVectorBase &matrixRow) const |
int | findExactMostViolatedMinCover (int nCols, int row, CoinPackedVector &krow, double b, double *xstar, CoinPackedVector &cover, CoinPackedVector &remainder) const |
int | findLPMostViolatedMinCover (int nCols, int row, CoinPackedVector &krow, double &b, double *xstar, CoinPackedVector &cover, CoinPackedVector &remainder) const |
int | findGreedyCover (int row, CoinPackedVector &krow, double &b, double *xstar, CoinPackedVector &cover, CoinPackedVector &remainder) const |
find a minimum cover by a simple greedy approach | |
int | liftCoverCut (double &b, int nRowElem, CoinPackedVector &cover, CoinPackedVector &remainder, CoinPackedVector &cut) const |
lift the cover inequality | |
int | liftAndUncomplementAndAdd (double rowub, CoinPackedVector &krow, double &b, int *complement, int row, CoinPackedVector &cover, CoinPackedVector &remainder, OsiCuts &cs) const |
sequence-independent lift and uncomplement and add the resulting cut to the cut set | |
void | seqLiftAndUncomplementAndAdd (int nCols, double *xstar, int *complement, int row, int nRowElem, double &b, CoinPackedVector &cover, CoinPackedVector &remainder, OsiCuts &cs) const |
sequence-dependent lift, uncomplement and add the resulting cut to the cut set | |
void | liftUpDownAndUncomplementAndAdd (int nCols, double *xstar, int *complement, int row, int nRowElem, double &b, CoinPackedVector &fracCover, CoinPackedVector &atOne, CoinPackedVector &remainder, OsiCuts &cs) const |
sequence-dependent lift binary variables either up or down, uncomplement and add to the cut set | |
int | findPseudoJohnAndEllisCover (int row, CoinPackedVector &krow, double &b, double *xstar, CoinPackedVector &cover, CoinPackedVector &remainder) const |
find a cover using a variation of the logic found in OSL (w/o SOS) | |
int | findJohnAndEllisCover (int row, CoinPackedVector &krow, double &b, double *xstar, CoinPackedVector &fracCover, CoinPackedVector &atOnes, CoinPackedVector &remainder) const |
find a cover using the basic logic found in OSL (w/o SOS) | |
int | exactSolveKnapsack (int n, double c, double const *pp, double const *ww, double &z, int *x) const |
Private Attributes | |
Private member data | |
double | epsilon_ |
epsilon | |
double | epsilon2_ |
Tolerance to use for violation - bigger than epsilon_. | |
double | onetol_ |
1-epsilon | |
int | maxInKnapsack_ |
Maximum in knapsack. | |
int | numRowsToCheck_ |
int * | rowsToCheck_ |
Friends | |
void | CglKnapsackCoverUnitTest (const OsiSolverInterface *siP, const std::string mpdDir) |
Definition at line 11 of file CglKnapsackCover.hpp.
|
deriveAKnapsack returns 1 if it is able to derive a (canonical) knapsack inequality in binary variables of the form ax<=b from the rowIndex-th row in the model, returns 0 otherwise. Definition at line 738 of file CglKnapsackCover.cpp. References CoinPackedVector::clear(), epsilon_, CoinPackedVector::getElements(), CoinPackedVector::getIndices(), CoinPackedVector::getNumElements(), CoinPackedVector::insert(), and onetol_. Referenced by generateCuts().
00750 { 00751 int i; 00752 00753 krow.clear(); 00754 00755 // if the matrixRow represent a ge inequality, then 00756 // leMatrixRow == -matrixRow // otherwise 00757 // leMatrixRow == matrixRow. 00758 00759 CoinPackedVector leMatrixRow(numberElements,index,element); 00760 00761 double maxKrowElement = -DBL_MAX; 00762 double minKrowElement = DBL_MAX; 00763 00764 00765 if (treatAsLRow) { 00766 // treat as if L row 00767 } else { 00768 // treat as if G row 00769 b=-b; 00770 std::transform(leMatrixRow.getElements(), 00771 leMatrixRow.getElements() + leMatrixRow.getNumElements(), 00772 leMatrixRow.getElements(), 00773 std::negate<double>()); 00774 } 00775 00776 // nBinUnsat is a counter for the number of unsatisfied 00777 // (i.e. fractional) binary vars 00778 int nBinUnsat =0; 00779 const double * colupper = si.getColUpper(); 00780 const double * collower = si.getColLower(); 00781 00782 // At this point, leMatrixRow and b represent a le inequality in general 00783 // variables. 00784 // To derive a canonical knapsack inequality in free binary variable, 00785 // process out the continuous & non-binary integer & fixed binary variables. 00786 // If the non-free-binary variables can be appropriately bounded, 00787 // net them out of the constraint, otherwise abandon this row and return 0. 00788 00789 const int * indices = leMatrixRow.getIndices(); 00790 const double * elements = leMatrixRow.getElements(); 00791 // for every variable in the constraint 00792 for (i=0; i<leMatrixRow.getNumElements(); i++){ 00793 // if the variable is not a free binary var 00794 if ( !si.isFreeBinary(indices[i]) ) { 00795 // and the coefficient is strictly negative 00796 if(elements[i]<-epsilon_){ 00797 // and the variable has a finite upper bound 00798 if (colupper[indices[i]] < si.getInfinity()){ 00799 // then replace the variable with its upper bound. 00800 b=b-elements[i]*colupper[indices[i]]; 00801 } 00802 else { 00803 return 0; 00804 } 00805 } 00806 // if the coefficient is strictly positive 00807 else if(elements[i]>epsilon_){ 00808 // and the variable has a finite lower bound 00809 if (collower[indices[i]] > -si.getInfinity()){ 00810 // then replace the variable with its lower bound. 00811 b=b-elements[i]*collower[indices[i]]; 00812 } 00813 else { 00814 return 0; 00815 } 00816 } 00817 // note: if the coefficient is zero, the variable is not included in the 00818 // knapsack inequality. 00819 } 00820 // else the variable is a free binary var and is included in the knapsack 00821 // inequality. 00822 // note: the variable is included regardless of its solution value to the 00823 // lp relaxation. 00824 else{ 00825 krow.insert(indices[i], elements[i]); 00826 00827 // if the binary variable is unsatified (i.e. has fractional value), 00828 // increment the counter. 00829 if(xstar[indices[i]] > epsilon_ && xstar[indices[i]] < onetol_) 00830 nBinUnsat++; 00831 00832 // keep track of the largest and smallest elements in the knapsack 00833 // (the idea is if there is not a lot of variation in the knapsack 00834 // coefficients, it is unlikely we will find a violated minimal 00835 // cover from this knapsack so don't even bother trying) 00836 if (fabs(elements[i]) > maxKrowElement) 00837 maxKrowElement = fabs(elements[i]); 00838 if (fabs(elements[i]) < minKrowElement) 00839 minKrowElement = fabs(elements[i]); 00840 } 00841 } 00842 00843 // If there's little variation in the knapsack coefficients, return 0. 00844 // If there are no unsatisfied binary variables, return. 00845 // If there's only one binary, return. 00846 // ToDo: but why return if 2 binary? ...there was some assumption in the 00847 // findVioMinCover..(?) 00848 // Anyway probing will probably find something 00849 if (krow.getNumElements() < 3 || 00850 nBinUnsat == 0 || 00851 maxKrowElement-minKrowElement < 1.0e-3*maxKrowElement ) { 00852 return 0; 00853 } 00854 00855 // However if we do decide to do when count is two - look carefully 00856 if (krow.getNumElements()==2) { 00857 const int * indices = krow.getIndices(); 00858 double * elements = krow.getElements(); 00859 double sum=0.0; 00860 for(i=0; i<2; i++){ 00861 int iColumn = indices[i]; 00862 sum += elements[i]*xstar[iColumn]; 00863 } 00864 if (sum<b-1.0e-4) { 00865 return 0; 00866 } else { 00867 printf("*** Doubleton Row is "); 00868 for(i=0; i<2; i++){ 00869 int iColumn = indices[i]; 00870 sum += elements[i]*xstar[iColumn]; 00871 printf("%d (coeff = %g, value = %g} ",indices[i], 00872 elements[i],xstar[iColumn]); 00873 } 00874 printf("<= %g - go for it\n",b); 00875 } 00876 } 00877 00878 00879 // At this point krow and b represent a le inequality in binary variables. 00880 // To obtain an le inequality with all positive coefficients, complement 00881 // any variable with a negative coefficient by changing the sign of 00882 // the coefficient, adjusting the rhs, and adjusting xstar, the column 00883 // solution vector. 00884 { 00885 const int s = krow.getNumElements(); 00886 const int * indices = krow.getIndices(); 00887 double * elements = krow.getElements(); 00888 for(i=0; i<s; i++){ 00889 if (elements[i] < -epsilon_) { 00890 complement[indices[i]]= 1; 00891 elements[i] *= -1; 00892 b+=elements[i]; 00893 xstar[indices[i]]=1.0-xstar[indices[i]]; 00894 } 00895 } 00896 } 00897 00898 // Quick feasibility check. 00899 // If the problem is infeasible, add an infeasible col cut to cut set 00900 // e.g. one that has lb > ub. 00901 // TODO: test this scenario in BCP 00902 if (b < 0 ){ 00903 OsiColCut cc; 00904 int index = krow.getIndices()[0]; 00905 const double fakeLb = colupper[index] + 1.0;; // yes, colupper. 00906 const double fakeUb = collower[index]; 00907 assert( fakeUb < fakeLb ); 00908 cc.setLbs( 1, &index, &fakeLb); 00909 cc.setUbs( 1, &index, &fakeLb); 00910 cc.setEffectiveness(DBL_MAX); 00911 cs.insert(cc); 00912 #ifdef PRINT_DEBUG 00913 printf("Cgl: Problem is infeasible\n"); 00914 #endif 00915 } 00916 00917 // At this point, krow and b represent a le inequality with postive 00918 // coefficients. 00919 // If any coefficient a_j > b, then x_j = 0, return 0 00920 // If any complemented var has coef a_j > b, then x_j = 1, return 0 00921 int fixed = 0; 00922 CoinPackedVector fixedBnd; 00923 for(i=0; i<krow.getNumElements(); i++){ 00924 if (krow.getElements()[i]> b){ 00925 fixedBnd.insert(krow.getIndices()[i],complement[krow.getIndices()[i]]); 00926 #ifdef PRINT_DEBUG 00927 printf("Variable %i being fixed to %i due to row %i.\n", 00928 krow.getIndices()[i],complement[krow.getIndices()[i]],rowIndex); 00929 #endif 00930 fixed = 1; 00931 } 00932 } 00933 00934 // After all possible variables are fixed by adding a column cut with 00935 // equivalent lower and upper bounds, return 00936 if (fixed) { 00937 OsiColCut cc; 00938 cc.setLbs(fixedBnd); 00939 cc.setUbs(fixedBnd); 00940 cc.setEffectiveness(DBL_MAX); 00941 return 0; 00942 } 00943 00944 return 1; 00945 } |
|
A C-style implementation of the Horowitz-Sahni exact solution procedure for solving knapsack problem. (ToDo: implement the more efficient dynamic programming approach) (Reference: Martello and Toth, Knapsack Problems, Wiley, 1990, p30.) Definition at line 2549 of file CglKnapsackCover.cpp. References epsilon2_. Referenced by findExactMostViolatedMinCover(), liftUpDownAndUncomplementAndAdd(), and seqLiftAndUncomplementAndAdd().
02556 { 02557 // The knapsack problem is to find: 02558 02559 // max {sum(j=1,n) p_j*x_j st. sum (j=1,n)w_j*x_j <= c, x binary} 02560 02561 // Notation: 02562 // xhat : current solution vector 02563 // zhat : current solution value = sum (j=1,n) p_j*xhat_j 02564 // chat : current residual capacity = c - sum (j=1,n) w_j*xhat_j 02565 // x : best solution so far, n-vector. 02566 // z : value of the best solution so far = sum (j=1,n) p_j*x_j 02567 02568 02569 // Input: n, the number of variables; 02570 // c, the rhs; 02571 // p, n-vector of objective func. coefficients; 02572 // w, n-vector of the row coeff. 02573 02574 // Output: z, the optimal objective function value; 02575 // x, the optimal (binary) solution n-vector; 02576 02577 // Assumes items are sorted p_1/w_1 >= p_2/w_2 >= ... >= p_n/w_n 02578 02579 memset(x, 0, (n)*sizeof(int)); 02580 int * xhat = new int[n+1]; 02581 memset(xhat, 0, (n+1)*sizeof(int)); 02582 int j; 02583 02584 // set up: adding the extra element and 02585 // accounting for the FORTRAN vs C difference in indexing arrays. 02586 double * p = new double[n+2]; 02587 double * w = new double[n+2]; 02588 int ii; 02589 for (ii=1; ii<n+1; ii++){ 02590 p[ii]=pp[ii-1]; 02591 w[ii]=ww[ii-1]; 02592 } 02593 02594 // 1. initialize 02595 double zhat = 0.0; 02596 z = 0.0; 02597 double chat = c+epsilon2_; 02598 p[n+1] = 0.0; 02599 w[n+1]= DBL_MAX; 02600 j=1; 02601 02602 while (1){ 02603 // 2. compute upper bound u 02604 // "find r = min {i: sum k=j,i w_k>chat};" 02605 ii=j; 02606 double wSemiSum = w[j]; 02607 double pSemiSum = p[j]; 02608 while (wSemiSum <= chat && ii<n+2){ 02609 ii++; 02610 wSemiSum+=w[ii]; 02611 pSemiSum+=p[ii]; 02612 } 02613 if (ii==n+2){ 02614 printf("Exceeded iterator limit. Aborting...\n"); 02615 abort(); 02616 } 02617 // r = ii at this point 02618 wSemiSum -= w[ii]; 02619 pSemiSum -= p[ii]; 02620 double u = pSemiSum + floor((chat - wSemiSum)*p[ii]/w[ii]); 02621 02622 // "if (z >= zhat + u) goto 5: backtrack;" 02623 if (!(z >= zhat + u)) { 02624 do { 02625 // 3. perform a forward step 02626 while (w[j] <= chat){ 02627 chat = chat - w[j]; 02628 zhat = zhat + p[j]; 02629 xhat[j] = 1; 02630 j+=1; 02631 } 02632 if (j<=n) { 02633 xhat[j]= 0; 02634 j+=1; 02635 } 02636 } while(j==n); 02637 02638 // "if (j<n) goto 2: compute_ub;" 02639 if (j<n) 02640 continue; 02641 02642 // 4. up date the best solution so far 02643 if (zhat > z) { 02644 z=zhat; 02645 int k; 02646 for (k=0; k<n; k++){ 02647 x[k]=xhat[k+1]; 02648 } 02649 } 02650 j=n; 02651 if (xhat[n] == 1){ 02652 chat = chat+ w[n]; 02653 zhat = zhat-p[n]; 02654 xhat[n]=0; 02655 } 02656 } 02657 // 5. backtrack 02658 // "find i=max{k<j:xhat[k]=1};" 02659 int i=j-1; 02660 while (!(xhat[i]==1) && i>0){ 02661 i--; 02662 } 02663 02664 // "if (no such i exists) return;" 02665 if (i==0){ 02666 delete [] p; 02667 delete [] w; 02668 delete [] xhat; 02669 return 1; 02670 } 02671 02672 chat = chat + w[i]; 02673 zhat=zhat -p[i]; 02674 xhat[i]=0; 02675 j=i+1; 02676 // "goto 2: compute_ub;" 02677 } 02678 } |
|
Find a violated minimal cover from a canonical form knapsack inequality by solving the -most- violated cover problem and postprocess to ensure minimality Definition at line 1196 of file CglKnapsackCover.cpp. References epsilon_, exactSolveKnapsack(), CoinPackedVector::getElements(), CoinPackedVector::getIndices(), CoinPackedVector::getNumElements(), CoinPackedVector::insert(), CoinPackedVector::reserve(), CoinPackedVector::sort(), CoinPackedVector::sortDecrElement(), CoinPackedVectorBase::sum(), and CoinPackedVector::truncate(). Referenced by generateCuts().
01204 { 01205 01206 // assumes the row is in canonical knapsack form 01207 01208 // A violated min.cover inequality exists if the 01209 // opt obj func value (oofv) < 1: 01210 // oofv = min sum (1-xstar_j)z_j 01211 // s.t. sum a_jz_j > b 01212 // x binary 01213 01214 // The vector z is the incidence vector 01215 // defines the set R and the cover inequality. 01216 // (sum j in R) x_j <= |R|-1 01217 01218 // This is the min version of the knapsack problem. 01219 // (note that strict inequalty...bleck) 01220 01221 // To obtain the max version, complement the z's, z_j=1-y_j and 01222 // adjust the constraint. 01223 01224 // oofv = (sum (1-xstar_j))- max sum (1-xstar)y_j 01225 // s.t. sum a_jy_j <= (sum over j)a_j - b (- EPSILON)] 01226 // y binary 01227 01228 // If oofv < 1, violated min cover inequality has 01229 // incidence vector z=1-y and rhs= num of nonzero's in z, i.e. 01230 // the number 0 in y. 01231 01232 // We solve the 0-1 knapsack problem by explicit ennumeration 01233 01234 double elementSum = krow.sum(); 01235 01236 // Redundant/useless adjusted rows should have been trapped in 01237 // transformation to canonical form of knapsack inequality 01238 if (elementSum < b + epsilon_) { 01239 #ifdef PRINT_DEBUG 01240 printf("Redundant/useless adjusted row\n"); 01241 #endif 01242 return -1; 01243 } 01244 01245 // Order krow in nonincreasing order of coefObj_j/a_j. 01246 // (1-xstar_1)/a_1 >= (1-xstar_2)/a_2 >= ... >= (1-xstar_n)/a_n 01247 // by defining this full-storage array "ratio" to be the external sort key. 01248 double * ratio= new double[nCols]; 01249 memset(ratio, 0, (nCols*sizeof(double))); 01250 01251 int i; 01252 { 01253 const int * indices = krow.getIndices(); 01254 const double * elements = krow.getElements(); 01255 for (i=0; i<krow.getNumElements(); i++){ 01256 if (fabs(elements[i])> epsilon_ ){ 01257 ratio[indices[i]]= (1.0-xstar[indices[i]]) / elements[i]; 01258 } 01259 else { 01260 ratio[indices[i]] = 0.0; 01261 } 01262 } 01263 } 01264 01265 // ToDo: would be nice to have sortkey NOT be full-storage vector 01266 CoinDecrSolutionOrdered dso(ratio); 01267 krow.sort(dso); 01268 01269 #ifdef CGL_DEBUG 01270 // sanity check 01271 for ( i=1; i<krow.getNumElements(); i++ ) { 01272 double ratioim1 = ratio[krow.getIndices()[i-1]]; 01273 double ratioi= ratio[krow.getIndices()[i]]; 01274 assert( ratioim1 >= ratioi ); 01275 } 01276 #endif 01277 01278 // Recall: 01279 // oofv = (sum (1-xstar_j))- max sum (1-xstar)y_j 01280 // s.t. sum a_jy_j <= (sum over j)a_j - b (- epsilon_)] 01281 // y binary 01282 01283 double objConst = 0.0; 01284 double exactOptVal = -1.0; 01285 int * exactOptSol = new int[krow.getNumElements()]; 01286 double * p = new double[krow.getNumElements()]; 01287 double * w = new double[krow.getNumElements()]; 01288 int kk; 01289 for (kk=0; kk<krow.getNumElements(); kk++){ 01290 p[kk]=1.0-xstar[krow.getIndices()[kk]]; 01291 w[kk]=krow.getElements()[kk]; 01292 objConst+=p[kk]; 01293 } 01294 01295 // vectors are indexed in ratioSortIndex order 01296 exactSolveKnapsack(krow.getNumElements(), (elementSum-b-epsilon_), p, w, 01297 exactOptVal, exactOptSol); 01298 01299 if(objConst-exactOptVal < 1){ 01300 cover.reserve(krow.getNumElements()); 01301 remainder.reserve(krow.getNumElements()); 01302 01303 // Partition the krow into the cover and the remainder. 01304 // The cover is complement of solution. 01305 double coverElementSum = 0; 01306 for(kk=0;kk<krow.getNumElements();kk++){ 01307 if(exactOptSol[kk]==0){ 01308 cover.insert(krow.getIndices()[kk],krow.getElements()[kk]); 01309 coverElementSum += krow.getElements()[kk]; 01310 } 01311 else { 01312 remainder.insert(krow.getIndices()[kk],krow.getElements()[kk]); 01313 } 01314 } 01315 01316 cover.sortDecrElement(); 01317 01318 // We have a violated cover inequality. 01319 // Construct a -minimal- violated cover 01320 // by testing and potentially tossing smallest 01321 // elements 01322 double oneLessCoverElementSum = 01323 coverElementSum - cover.getElements()[cover.getNumElements()-1]; 01324 while (oneLessCoverElementSum > b){ 01325 // move the excess cover member into the set of remainders 01326 remainder.insert(cover.getIndices()[cover.getNumElements()-1], 01327 cover.getElements()[cover.getNumElements()-1]); 01328 cover.truncate(cover.getNumElements()-1); 01329 oneLessCoverElementSum -= cover.getElements()[cover.getNumElements()-1]; 01330 } 01331 01332 #ifdef PRINT_DEBUG 01333 printf("Exact Most Violated Cover: row %i has cover of size %i.\n", 01334 row,cover.getNumElements()); 01335 //double sumCover = 0.0; 01336 for (i=0; i<cover.getNumElements(); i++){ 01337 printf("index %i, element %g, xstar value % g \n", 01338 cover.getIndices()[i], cover.getElements()[i], 01339 xstar[cover.getIndices()[i]]); 01340 //sumCover += cover.getElements()[i]; 01341 } 01342 printf("The b = %g, and the cover sum is %g\n\n", b, cover.sum() ); 01343 #endif 01344 01345 // local clean up 01346 delete [] exactOptSol; 01347 delete [] p; 01348 delete [] w; 01349 delete [] ratio; 01350 01351 return 1; // found an exact one 01352 } 01353 01354 // local clean up 01355 delete [] exactOptSol; 01356 delete [] p; 01357 delete [] w; 01358 delete [] ratio; 01359 01360 return 0; // didn't find an exact one 01361 01362 } |
|
Find the most violate minimum cover by solving the lp-relaxation of the most-violate-min-cover problem Definition at line 992 of file CglKnapsackCover.cpp. References epsilon_, CoinPackedVector::getElements(), CoinPackedVector::getIndices(), CoinPackedVector::getNumElements(), CoinPackedVector::insert(), CoinPackedVector::reserve(), CoinPackedVector::sort(), CoinPackedVector::sortDecrElement(), CoinPackedVectorBase::sum(), and CoinPackedVector::truncate(). Referenced by generateCuts().
01000 { 01001 01002 // Assumes krow and b describe a knapsack inequality in canonical form 01003 01004 // Given a knapsack inequality sum a_jx_j <= b, and optimal lp solution 01005 // xstart, a violated minimal cover inequality exists if the following 0-1 01006 // programming problem has an optimal objective function value (oofv) < 1 01007 // oofv = min sum (1-xstar_j)z_j 01008 // s.t. sum a_jz_j > b 01009 // z binary 01010 01011 // The vector z is an incidence vector, defining the cover R with the 01012 // associated cover inequality: 01013 // (sum j in R) x_j <= |R|-1 01014 01015 // This problem is itself a (min version of the) knapsack problem 01016 // but with a unsightly strict inequalty. 01017 01018 // To transform transform it into a max version, 01019 // complement the z's, z_j=1-y_j. 01020 // To compensate for the strict inequality, subtract epsilon from the rhs. 01021 01022 // oofv = (sum (1-xstar_j))- max sum (1-xstar)y_j 01023 // s.t. sum a_jy_j <= (sum over j)a_j - b (- EPSILON) 01024 // y binary 01025 01026 // If oofv < 1, then a violated min cover inequality has 01027 // incidence vector z with elements z_j=1-y_j and rhs= num of nonzero's in 01028 // z, i.e. the number of 0's in y. 01029 01030 // If the size of the knapsack is "small", we solve the problem exactly. 01031 // If the size of the knapsack is large, we solve the (simpler) lp relaxation 01032 // of the knapsack problem and postprocess to ensure the construction of a 01033 // minimimal cover. 01034 01035 // We also assume that testing/probing/fixing based on the knapsack structure 01036 // is done elsewhere. Only convenient-to-do sanity checks are done here. 01037 // (We do not assume that data is integer.) 01038 01039 double elementSum = krow.sum(); 01040 01041 // Redundant/useless adjusted rows should have been trapped in the 01042 // transformation to the canonical form of knapsack inequality 01043 if (elementSum < b + epsilon_) { 01044 return -1; 01045 } 01046 01047 // Order krow in nonincreasing order of coefObj_j/a_j. 01048 // (1-xstar_1)/a_1 >= (1-xstar_2)/a_2 >= ... >= (1-xstar_n)/a_n 01049 // by defining this full-storage array "ratio" to be the external sort key. 01050 double * ratio= new double[nCols]; 01051 memset(ratio, 0, (nCols*sizeof(double))); 01052 01053 int i; 01054 for (i=0; i<krow.getNumElements(); i++){ 01055 if (fabs(krow.getElements()[i])> epsilon_ ){ 01056 ratio[krow.getIndices()[i]]= 01057 (1.0-xstar[krow.getIndices()[i]]) / (krow.getElements()[i]); 01058 } 01059 else { 01060 ratio[krow.getIndices()[i]] = 0.0; 01061 } 01062 } 01063 01064 // ToDo: would be nice to have sortkey NOT be full-storage vector 01065 CoinDecrSolutionOrdered dso(ratio); 01066 krow.sort(dso); 01067 01068 // Find the "critical" element index "r" in the knapsack lp solution 01069 int r = 0; 01070 double sum = krow.getElements()[0]; 01071 while ( sum <= (elementSum - b - epsilon_ ) ){ 01072 r++; 01073 sum += krow.getElements()[r]; 01074 } 01075 01076 // Note: It is possible the r=0, and you get a violated minimal cover 01077 // if (r=0), then you've got a var with a really large coeff. compared 01078 // to the rest of the row. 01079 // r=0 says trivially that the 01080 // sum of ALL the binary vars in the row <= (cardinality of all the set -1) 01081 // Note: The cover may not be minimal if there are alternate optimals to the 01082 // maximization problem, so the cover must be post-processed to ensure 01083 // minimality. 01084 01085 // "r" is the critical element 01086 // The lp relaxation column solution is: 01087 // y_j = 1 for j=0,...,(r-1) 01088 // y_r = (elementSum - b - sum + krow.element()[r])/krow.element()[r] 01089 // y_j = 0 for j=r+1,...,krow.getNumElements() 01090 01091 // The number of nonzeros in the lp column solution is r+1 01092 01093 // if oofv to the lp knap >= 1, then no violated min cover is possible 01094 int nCover; 01095 01096 double lpoofv=0.0; 01097 for (i=r+1; i<krow.getNumElements(); i++){ 01098 lpoofv += (1-xstar[krow.getIndices()[i]]); 01099 } 01100 double ipofv = lpoofv + (1-xstar[krow.getIndices()[r]]); 01101 01102 // Couldn't find an lp violated min cover inequality 01103 if (ipofv > 1.0 - epsilon_){ 01104 delete [] ratio; 01105 return -1; 01106 } 01107 else { 01108 // Partition knapsack into cover and noncover (i.e. remainder) 01109 // pieces 01110 nCover = krow.getNumElements() - r; 01111 double coverSum =0.0; 01112 cover.reserve(nCover); 01113 remainder.reserve(r); 01114 01115 for (i=r; i<krow.getNumElements(); i++){ 01116 cover.insert(krow.getIndices()[i],krow.getElements()[i]); 01117 coverSum += krow.getElements()[i]; 01118 } 01119 for (i=0; i<r; i++){ 01120 remainder.insert(krow.getIndices()[i],krow.getElements()[i]); 01121 } 01122 01123 #ifdef PRINT_DEBUG 01124 if (coverSum <= b){ 01125 printf("The identified cover is NOT a cover\n"); 01126 abort(); 01127 } 01128 #endif 01129 01130 // Sort cover in terms of knapsack row coefficients 01131 cover.sortDecrElement(); 01132 01133 01134 // We have a violated cover inequality. 01135 // Construct a -minimal- violated cover 01136 // by testing and potentially tossing smallest 01137 // elements 01138 double oneLessCoverSum = coverSum - cover.getElements()[nCover-1]; 01139 while (oneLessCoverSum > b){ 01140 // move the excess cover member into the set of remainders 01141 remainder.insert(cover.getIndices()[nCover-1], 01142 cover.getElements()[nCover-1]); 01143 cover.truncate(nCover-1); 01144 nCover--; 01145 oneLessCoverSum -= cover.getElements()[nCover-1]; 01146 } 01147 01148 if (nCover<2){ 01149 #ifdef PRINT_DEBUG 01150 printf("nCover < 2...aborting\n"); 01151 #endif 01152 abort(); 01153 } 01154 01155 #ifdef PRINT_DEBUG /* debug */ 01156 printf("\ 01157 Lp relax of most violated minimal cover: row %i has cover of size %i.\n", 01158 row,nCover); 01159 //double sumCover = 0.0; 01160 for (i=0; i<cover.getNumElements(); i++){ 01161 printf("index %i, element %g, xstar value % g \n", 01162 cover.getIndices()[i],cover.getElements()[i], 01163 xstar[cover.getIndices()[i]]); 01164 //sumCover += cover.getElements()[i]; 01165 } 01166 printf("The b = %g, and the cover sum is %g\n\n", b, cover.sum()); 01167 #endif 01168 01169 #ifdef P0201 01170 double ppsum=0.0; 01171 for (i=0; i<nCover; i++){ 01172 ppsum += p0201[krow.getIndices()[i]]; 01173 } 01174 01175 if (ppsum > nCover-1){ 01176 printf("\ 01177 \nBad cover from lp relax of most violated cover..aborting\n"); 01178 abort(); 01179 } 01180 #endif 01181 01182 /* clean up */ 01183 delete [] ratio; 01184 return 1; 01185 } 01186 } |
|
Generate knapsack cover cuts for the model of the solver interface, si. Insert the generated cuts into OsiCut, cs. Implements CglCutGenerator. Definition at line 25 of file CglKnapsackCover.cpp. References deriveAKnapsack(), epsilon_, findExactMostViolatedMinCover(), findGreedyCover(), findJohnAndEllisCover(), findLPMostViolatedMinCover(), findPseudoJohnAndEllisCover(), CoinPackedVector::getIndices(), CoinPackedVector::getNumElements(), liftAndUncomplementAndAdd(), liftUpDownAndUncomplementAndAdd(), maxInKnapsack_, numRowsToCheck_, seqLiftAndUncomplementAndAdd(), and CoinPackedVector::setVector().
00027 { 00028 // Get basic problem information 00029 int nRows=si.getNumRows(); 00030 int nCols=si.getNumCols(); 00031 00032 // Create working space for "canonical" knapsack inequality 00033 // - krow will contain the coefficients and indices of the 00034 // (potentially complemented) variables in the knapsack inequality. 00035 // - b is the rhs of knapsack inequality. 00036 // - complement[i] is 1 if the index i in krow refers to the complement 00037 // of the variable, and 0 otherwise. 00038 CoinPackedVector krow; 00039 double b=0.0; 00040 int * complement= new int[nCols]; 00041 00042 // Create a local copy of the column solution (colsol), call it xstar, and 00043 // inititalize it. 00044 // Assumes the lp-relaxation has been solved, and the solver interface 00045 // has a meaningful colsol. 00046 double * xstar= new double[nCols]; 00047 00048 // To allow for vub knapsacks 00049 int * thisColumnIndex = new int [nCols]; 00050 double * thisElement = new double[nCols]; 00051 int * back = new int[nCols]; 00052 00053 const double *colsol = si.getColSolution(); 00054 int k; 00055 // For each row point to vub variable 00056 // -1 if no vub 00057 // -2 if can skip row for knapsacks 00058 00059 int * vub = new int [nRows]; 00060 00061 // For each column point to vub row 00062 int * vubRow = new int [nCols]; 00063 double * vubValue = new double [nRows]; 00064 00065 // Take out all fixed 00066 double * effectiveUpper = new double [nRows]; 00067 double * effectiveLower = new double [nRows]; 00068 const double * colUpper = si.getColUpper(); 00069 const double * colLower = si.getColLower(); 00070 for (k=0; k<nCols; k++){ 00071 back[k]=-1; 00072 xstar[k]=colsol[k]; 00073 complement[k]=0; 00074 vubRow[k]=-1; 00075 if (si.isBinary(k)) { 00076 if (si.isFreeBinary(k)) { 00077 vubRow[k]=-2; 00078 } else { 00079 vubRow[k]=-10; 00080 } 00081 } else if (colUpper[k]==colLower[k]) { 00082 vubRow[k]=-10; // fixed 00083 } 00084 } 00085 00086 int rowIndex; 00087 int numberVub=0; 00088 00089 const CoinPackedMatrix * matrixByRow = si.getMatrixByRow(); 00090 const double * elementByRow = matrixByRow->getElements(); 00091 const int * column = matrixByRow->getIndices(); 00092 const int * rowStart = matrixByRow->getVectorStarts(); 00093 const int * rowLength = matrixByRow->getVectorLengths(); 00094 const double * rowUpper = si.getRowUpper(); 00095 const double * rowLower = si.getRowLower(); 00096 00097 // Scan all rows looking for possibles 00098 00099 for (rowIndex=0;rowIndex<nRows;rowIndex++) { 00100 int start = rowStart[rowIndex]; 00101 int end = start + rowLength[rowIndex]; 00102 double upRhs = rowUpper[rowIndex]; 00103 double loRhs = rowLower[rowIndex]; 00104 int numberContinuous=0; 00105 int numberBinary=0; 00106 int iCont=-1; 00107 double sum = 0.0; 00108 double valueContinuous=0.0; 00109 double valueBinary=0.0; 00110 int iBinary=-1; 00111 int j; 00112 for (j=start;j<end;j++) { 00113 int iColumn=column[j]; 00114 if (colUpper[iColumn]>colLower[iColumn]) { 00115 sum += colsol[iColumn]*elementByRow[j]; 00116 if (vubRow[iColumn]==-2) { 00117 // binary 00118 numberBinary++; 00119 valueBinary=elementByRow[j]; 00120 iBinary=iColumn; 00121 } else if (vubRow[iColumn]==-1) { 00122 // only use if not at bound 00123 if (colsol[iColumn]<colUpper[iColumn]-1.0e-6&& 00124 colsol[iColumn]>colLower[iColumn]+1.0e-6) { 00125 // possible 00126 iCont=iColumn; 00127 numberContinuous++; 00128 valueContinuous=elementByRow[j]; 00129 } else { 00130 // ** needs more thought 00131 numberContinuous ++; 00132 iCont=-1; 00133 } 00134 } else { 00135 // ** needs more thought 00136 numberContinuous ++; 00137 iCont=-1; 00138 if (colsol[iColumn]<colUpper[iColumn]-1.0e-6&& 00139 colsol[iColumn]>colLower[iColumn]+1.0e-6) { 00140 // already assigned 00141 numberContinuous ++; 00142 iCont=-1; 00143 } 00144 } 00145 } else { 00146 // fixed 00147 upRhs -= colLower[iColumn]*elementByRow[j]; 00148 loRhs -= colLower[iColumn]*elementByRow[j]; 00149 } 00150 } 00151 // see if binding 00152 effectiveUpper[rowIndex] = upRhs; 00153 effectiveLower[rowIndex] = loRhs; 00154 vubValue[rowIndex] = valueContinuous; 00155 bool possible = false; 00156 if (fabs(sum-upRhs)<1.0e-5) { 00157 possible=true; 00158 } else { 00159 effectiveUpper[rowIndex]=DBL_MAX; 00160 } 00161 if (fabs(sum-loRhs)<1.0e-5) { 00162 possible=true; 00163 } else { 00164 effectiveLower[rowIndex]=-DBL_MAX; 00165 } 00166 if (!numberBinary||numberBinary+numberContinuous>maxInKnapsack_) 00167 possible=false; 00168 if (possible) { 00169 // binding with binary 00170 if(numberContinuous==1&&iCont>=0) { 00171 // vub 00172 if (numberBinary==1) 00173 #ifdef PRINT_DEBUG 00174 printf("vub (by row %d) %g <= 0-1 %g * %d + %g * %d <= %g\n", 00175 rowIndex,effectiveLower[rowIndex],valueBinary,iBinary, 00176 valueContinuous,iCont,effectiveUpper[rowIndex]); 00177 #endif 00178 vubRow[iCont]=rowIndex; 00179 vub[rowIndex]=iCont; 00180 numberVub++; 00181 } else if (numberBinary>1) { 00182 // could be knapsack 00183 vub[rowIndex]=-1; 00184 } else { 00185 // no point looking at this row 00186 vub[rowIndex]=-2; 00187 } 00188 } else { 00189 // no point looking at this row 00190 vub[rowIndex]=-2; 00191 } 00192 } 00193 // Main loop 00194 int numCheck = 0; 00195 int* toCheck = 0; 00196 if (!rowsToCheck_) { 00197 toCheck = new int[nRows]; 00198 CoinIotaN(toCheck, nRows, 0); 00199 numCheck = nRows; 00200 } else { 00201 numCheck = numRowsToCheck_; 00202 toCheck = rowsToCheck_; 00203 } 00204 00205 // Set up number of tries for each row 00206 int ntry; 00207 if (numberVub) 00208 ntry=4; 00209 else 00210 ntry=2; 00211 //ntry=2; // switch off 00212 for (int ii=0; ii < numCheck; ++ii){ 00213 rowIndex = toCheck[ii]; 00214 if (rowIndex < 0 || rowIndex >= nRows) 00215 continue; 00216 if (vub[rowIndex]==-2) 00217 continue; 00218 00219 #ifdef PRINT_DEBUG 00220 std::cout << "CGL: Processing row " << rowIndex << std::endl; 00221 #endif 00222 00223 // Get a tight row 00224 // (want to be able to 00225 // experiment by turning this on and off) 00226 // 00227 // const double * pi=si.rowprice(); 00228 // if (fabs(pi[row]) < epsilon_){ 00229 // continue; 00230 // } 00231 00232 00234 // Derive a "canonical" knapsack // 00235 // inequality (in binary variables) // 00236 // from the model row in mixed integer variables // 00238 #ifdef CGL_DEBUG 00239 assert(!krow.getNumElements()); 00240 #endif 00241 double effectiveRhs[4]; 00242 double rhs[4]; 00243 double sign[]={0.0,0.0,-1.0,1.0}; 00244 bool rowType[] = {false,true,false,true}; 00245 effectiveRhs[0] = effectiveLower[rowIndex]; 00246 rhs[0]=rowLower[rowIndex]; 00247 effectiveRhs[2] = effectiveRhs[0]; 00248 rhs[2]= effectiveRhs[0]; 00249 effectiveRhs[1] = effectiveUpper[rowIndex]; 00250 rhs[1]=rowUpper[rowIndex]; 00251 effectiveRhs[3] = effectiveRhs[1]; 00252 rhs[3]= effectiveRhs[1]; 00253 int itry; 00254 #ifdef CGL_DEBUG 00255 int kcuts[4]; 00256 memset(kcuts,0,4*sizeof(int)); 00257 #endif 00258 for (itry=0;itry<ntry;itry++) { 00259 #ifdef CGL_DEBUG 00260 int nlast=cs.sizeRowCuts(); 00261 #endif 00262 // see if to skip 00263 if (fabs(effectiveRhs[itry])>1.0e20) 00264 continue; 00265 int length = rowLength[rowIndex]; 00266 memcpy(thisColumnIndex,column+rowStart[rowIndex],length*sizeof(int)); 00267 memcpy(thisElement,elementByRow+rowStart[rowIndex], 00268 length*sizeof(double)); 00269 b=rhs[itry]; 00270 if (itry>1) { 00271 // see if we would be better off relaxing 00272 int i; 00273 // mark columns 00274 int length2=length; // for new length 00275 int numberReplaced=0; 00276 for (i=0;i<length;i++) { 00277 int iColumn = thisColumnIndex[i]; 00278 back[thisColumnIndex[i]]=i; 00279 if (vubRow[iColumn]==-10) { 00280 // fixed - take out 00281 thisElement[i]=0.0; 00282 } 00283 } 00284 for (i=0;i<length;i++) { 00285 int iColumn = thisColumnIndex[i]; 00286 int iRow = vubRow[iColumn]; 00287 if (iRow>=0&&vub[iRow]==iColumn&&iRow!=rowIndex) { 00288 double useRhs=0.0; 00289 double vubCoefficient = vubValue[iRow]; 00290 double thisCoefficient = thisElement[i]; 00291 int replace = 0; 00292 // break it out - may be able to do better 00293 if (sign[itry]*thisCoefficient>0.0) { 00294 // we want valid lower bound on continuous 00295 if (effectiveLower[iRow]>-1.0e20&&vubCoefficient>0.0) 00296 replace=-1; 00297 else if (effectiveUpper[iRow]<1.0e20&&vubCoefficient<0.0) 00298 replace=1; 00299 } else { 00300 // we want valid upper bound on continuous 00301 if (effectiveLower[iRow]>-1.0e20&&vubCoefficient<0.0) 00302 replace=-1; 00303 else if (effectiveUpper[iRow]<1.0e20&&vubCoefficient>0.0) 00304 replace=1; 00305 } 00306 if (replace) { 00307 numberReplaced++; 00308 if (replace<0) 00309 useRhs = effectiveLower[iRow]; 00310 else 00311 useRhs = effectiveUpper[iRow]; 00312 // now replace (just using vubRow==-2) 00313 // delete continuous 00314 thisElement[i]=0.0; 00315 double scale = thisCoefficient/vubCoefficient; 00316 // modify rhs 00317 b -= scale*useRhs; 00318 int start = rowStart[iRow]; 00319 int end = start+rowLength[iRow]; 00320 int j; 00321 for (j=start;j<end;j++) { 00322 int iColumn = column[j]; 00323 if (vubRow[iColumn]==-2) { 00324 double change = scale*elementByRow[j]; 00325 int iBack = back[iColumn]; 00326 if (iBack<0) { 00327 // element does not exist 00328 back[iColumn]=length2; 00329 thisElement[length2]=-change; 00330 thisColumnIndex[length2++]=iColumn; 00331 } else { 00332 // element does exist 00333 thisElement[iBack] -= change; 00334 } 00335 } 00336 } 00337 } 00338 } 00339 } 00340 if (numberReplaced) { 00341 length=0; 00342 for (i=0;i<length2;i++) { 00343 int iColumn = thisColumnIndex[i]; 00344 back[iColumn]=-1; // un mark 00345 if (thisElement[i]) { 00346 thisElement[length]=thisElement[i]; 00347 thisColumnIndex[length++]=iColumn; 00348 } 00349 } 00350 if (length>maxInKnapsack_) 00351 continue; // too long 00352 } else { 00353 for (i=0;i<length;i++) { 00354 int iColumn = thisColumnIndex[i]; 00355 back[iColumn]=-1; // un mark 00356 } 00357 continue; // no good 00358 } 00359 } 00360 if (!deriveAKnapsack(si, cs, krow, rowType[itry], b, complement, 00361 xstar, rowIndex, 00362 length,thisColumnIndex,thisElement)) { 00363 00364 // Reset local data and continue to the next iteration 00365 // of the rowIndex-loop 00366 for(k=0; k<krow.getNumElements(); k++) { 00367 if (complement[krow.getIndices()[k]]){ 00368 xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]]; 00369 complement[krow.getIndices()[k]]=0; 00370 } 00371 } 00372 krow.setVector(0,NULL,NULL); 00373 continue; 00374 } 00375 #ifdef PRINT_DEBUG 00376 { 00377 // Get the sense of the row 00378 int i; 00379 printf("rhs sense %c rhs %g\n",si.getRowSense()[rowIndex], 00380 si.getRightHandSide()[rowIndex]); 00381 const int * indices = si.getMatrixByRow()->getVector(rowIndex).getIndices(); 00382 const double * elements = si.getMatrixByRow()->getVector(rowIndex).getElements(); 00383 // for every variable in the constraint 00384 for (i=0; i<si.getMatrixByRow()->getVector(rowIndex).getNumElements(); i++){ 00385 printf("%d (s=%g) %g, ",indices[i],colsol[indices[i]],elements[i]); 00386 } 00387 printf("\n"); 00388 } 00389 #endif 00390 00392 // Look for a series of // 00393 // different types of minimal covers. // 00394 // If a minimal cover is found, // 00395 // lift the associated minimal cover inequality, // 00396 // uncomplement the vars // 00397 // and add it to the cut set. // 00398 // After the last type of cover is tried, // 00399 // restore xstar values // 00401 00403 // Try to generate a violated // 00404 // minimal cover greedily from fractional vars // 00406 00407 CoinPackedVector cover, remainder; 00408 00409 00410 if (findGreedyCover(rowIndex, krow, b, xstar, cover, remainder) == 1){ 00411 00412 // Lift cover inequality and add to cut set 00413 if (!liftAndUncomplementAndAdd(rowUpper[rowIndex], krow, b, 00414 complement, rowIndex, cover, 00415 remainder, cs)) { 00416 // Reset local data and continue to the next iteration 00417 // of the rowIndex-loop 00418 // I am not sure this is needed but I am just being careful 00419 for(k=0; k<krow.getNumElements(); k++) { 00420 if (complement[krow.getIndices()[k]]){ 00421 xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]]; 00422 complement[krow.getIndices()[k]]=0; 00423 } 00424 } 00425 krow.setVector(0,NULL,NULL); 00426 continue; 00427 } 00428 } 00429 00430 00432 // Try to generate a violated // 00433 // minimal cover using pseudo John and Ellis logic // 00435 00436 // Reset the cover and remainder 00437 cover.setVector(0,NULL,NULL); 00438 remainder.setVector(0,NULL,NULL); 00439 00440 if (findPseudoJohnAndEllisCover(rowIndex, krow, b, 00441 xstar, cover, remainder) == 1){ 00442 00443 // (Sequence Independent) Lift cover inequality and add to cut set 00444 if (!liftAndUncomplementAndAdd(rowUpper[rowIndex], krow, b, 00445 complement, rowIndex, cover, 00446 remainder, cs)) { 00447 // Reset local data and continue to the next iteration 00448 // of the rowIndex-loop 00449 // I am not sure this is needed but I am just being careful 00450 for(k=0; k<krow.getNumElements(); k++) { 00451 if (complement[krow.getIndices()[k]]){ 00452 xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]]; 00453 complement[krow.getIndices()[k]]=0; 00454 } 00455 } 00456 krow.setVector(0,NULL,NULL); 00457 continue; 00458 } 00459 00460 // Skip experiment for now... 00461 #if 0 00462 // experimenting here... 00463 // (Sequence Dependent) Lift cover inequality and add to cut set 00464 seqLiftAndUncomplementAndAdd(nCols, xstar, complement, rowIndex, 00465 krow.getNumElements(), b, cover, remainder, 00466 cs); 00467 #endif 00468 } 00469 00470 00471 00473 // Try to generate cuts using covers of unsat // 00474 // vars on reduced krows with John and Ellis logic // 00476 CoinPackedVector atOnes; 00477 CoinPackedVector fracCover; // different than cover 00478 00479 // reset the remainder 00480 remainder.setVector(0,NULL,NULL); 00481 00482 00483 if (findJohnAndEllisCover(rowIndex, krow, b, 00484 xstar, fracCover, atOnes, remainder) == 1){ 00485 00486 // experimenting here... 00487 // Sequence Dependent Lifting up on remainders and lifting down on the 00488 // atOnes 00489 liftUpDownAndUncomplementAndAdd(nCols, xstar, complement, rowIndex, 00490 krow.getNumElements(), b, fracCover, 00491 atOnes, remainder, cs); 00492 } 00493 00494 00496 // Try to generate a violated // 00497 // minimal cover by considering the // 00498 // most violated cover problem // 00500 00501 00502 // reset cover and remainder 00503 cover.setVector(0,NULL,NULL); 00504 remainder.setVector(0,NULL,NULL); 00505 00506 // if the size of the krow is "small", 00507 // use an exact algorithm to find the most violated (minimal) cover, 00508 // else, 00509 // use an lp-relaxation to find the most violated (minimal) cover. 00510 if(krow.getNumElements()<=15){ 00511 if (findExactMostViolatedMinCover(nCols, rowIndex, krow, b, 00512 xstar, cover, remainder) == 1){ 00513 00514 // Lift cover inequality and add to cut set 00515 if (!liftAndUncomplementAndAdd(rowUpper[rowIndex], krow, b, 00516 complement, rowIndex, cover, remainder, cs)) { 00517 // Reset local data and continue to the next iteration 00518 // of the rowIndex-loop 00519 // I am not sure this is needed but I am just being careful 00520 for(k=0; k<krow.getNumElements(); k++) { 00521 if (complement[krow.getIndices()[k]]){ 00522 xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]]; 00523 complement[krow.getIndices()[k]]=0; 00524 } 00525 } 00526 krow.setVector(0,NULL,NULL); 00527 continue; 00528 } 00529 } 00530 } 00531 else { 00532 if (findLPMostViolatedMinCover(nCols, rowIndex, krow, b, 00533 xstar, cover, remainder) == 1){ 00534 00535 // Lift cover inequality and add to cut set 00536 if (!liftAndUncomplementAndAdd(rowUpper[rowIndex], krow, b, 00537 complement, rowIndex, cover, remainder, cs)) { 00538 // Reset local data and continue to the next iteration 00539 // of the rowIndex-loop 00540 // I am not sure this is needed but I am just being careful 00541 for(k=0; k<krow.getNumElements(); k++) { 00542 if (complement[krow.getIndices()[k]]){ 00543 xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]]; 00544 complement[krow.getIndices()[k]]=0; 00545 } 00546 } 00547 krow.setVector(0,NULL,NULL); 00548 continue; 00549 } 00550 } 00551 } 00552 00553 00554 00555 // Reset xstar and complement to their initialized values for the next 00556 // go-around 00557 int k; 00558 if (fabs(b-rowUpper[rowIndex]) > epsilon_) { 00559 for(k=0; k<krow.getNumElements(); k++) { 00560 if (complement[krow.getIndices()[k]]){ 00561 xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]]; 00562 complement[krow.getIndices()[k]]=0; 00563 } 00564 } 00565 } 00566 krow.setVector(0,NULL,NULL); 00567 #ifdef CGL_DEBUG 00568 int nnow = cs.sizeRowCuts(); 00569 if (nnow>nlast) { 00570 const OsiRowCutDebugger * debugger = si.getRowCutDebugger(); 00571 if (debugger&&debugger->onOptimalPath(si)) { 00572 // check cuts okay 00573 int k; 00574 for (k=nlast;k<nnow;k++) { 00575 OsiRowCut rc=cs.rowCut(k); 00576 if(debugger->invalidCut(rc)) { 00577 printf("itry %d, rhs %g, length %d\n",itry,rhs[itry],length); 00578 int i; 00579 for (i=0;i<length;i++) { 00580 int iColumn = thisColumnIndex[i]; 00581 printf("column %d, coefficient %g, value %g, bounds %g %g\n",iColumn, 00582 thisElement[i],colsol[iColumn],colLower[iColumn], 00583 colUpper[iColumn]); 00584 } 00585 if (itry>1) { 00586 int length = rowLength[rowIndex]; 00587 memcpy(thisColumnIndex,column+rowStart[rowIndex], 00588 length*sizeof(int)); 00589 memcpy(thisElement,elementByRow+rowStart[rowIndex], 00590 length*sizeof(double)); 00591 printf("Original row had rhs %g and length %d\n", 00592 (itry==2 ? rowLower[rowIndex] :rowUpper[rowIndex]), 00593 length); 00594 for (i=0;i<length;i++) { 00595 int iColumn = thisColumnIndex[i]; 00596 printf("column %d, coefficient %g, value %g, bounds %g %g\n",iColumn, 00597 thisElement[i],colsol[iColumn],colLower[iColumn], 00598 colUpper[iColumn]); 00599 } 00600 } 00601 assert(!debugger->invalidCut(rc)); 00602 } 00603 } 00604 } 00605 if (itry>1&&nnow-nlast>kcuts[itry-2]) { 00606 printf("itry %d gave %d cuts as against %d for itry %d\n", 00607 itry,nnow-nlast,kcuts[itry-2],itry-2); 00608 } 00609 kcuts[itry]=nnow-nlast; 00610 nlast=nnow; 00611 } 00612 #endif 00613 } 00614 } 00615 // Clean up: free allocated memory 00616 if (toCheck != rowsToCheck_) 00617 delete[] toCheck; 00618 delete[] xstar; 00619 delete[] complement; 00620 delete [] thisColumnIndex; 00621 delete [] thisElement; 00622 delete [] back; 00623 delete [] vub; 00624 delete [] vubRow; 00625 delete [] vubValue; 00626 delete [] effectiveLower; 00627 delete [] effectiveUpper; 00628 } |
|
A method to set which rows should be tested for knapsack covers Definition at line 631 of file CglKnapsackCover.cpp. References numRowsToCheck_.
00632 { 00633 if (rowsToCheck_) 00634 delete[] rowsToCheck_; 00635 numRowsToCheck_ = num; 00636 if (num > 0) { 00637 rowsToCheck_ = new int[num]; 00638 CoinCopyN(ind, num, rowsToCheck_); 00639 } 00640 } |
|
A function that tests the methods in the CglKnapsackCover 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 20 of file CglKnapsackCoverTest.cpp.
00023 { 00024 int i; 00025 CoinRelFltEq eq(0.000001); 00026 00027 // Test default constructor 00028 { 00029 CglKnapsackCover kccGenerator; 00030 } 00031 00032 // Test copy & assignment 00033 { 00034 CglKnapsackCover rhs; 00035 { 00036 CglKnapsackCover kccGenerator; 00037 CglKnapsackCover cgC(kccGenerator); 00038 rhs=kccGenerator; 00039 } 00040 } 00041 00042 00043 // test exactSolveKnapsack 00044 { 00045 CglKnapsackCover kccg; 00046 const int n=7; 00047 double c=50; 00048 double p[n] = {70,20,39,37,7,5,10}; 00049 double w[n] = {31, 10, 20, 19, 4, 3, 6}; 00050 double z; 00051 int x[n]; 00052 int exactsol = kccg.exactSolveKnapsack(n, c, p, w, z, x); 00053 assert(exactsol==1); 00054 assert (z == 107); 00055 assert (x[0]==1); 00056 assert (x[1]==0); 00057 assert (x[2]==0); 00058 assert (x[3]==1); 00059 assert (x[4]==0); 00060 assert (x[5]==0); 00061 assert (x[6]==0); 00062 } 00063 00064 /* 00065 // Testcase /u/rlh/osl2/mps/scOneInt.mps 00066 // Model has 3 continous, 2 binary, and 1 general 00067 // integer variable. 00068 { 00069 OsiSolverInterface * siP = baseSiP->clone(); 00070 int * complement=NULL; 00071 double * xstar=NULL; 00072 00073 siP->readMps("../Mps/scOneInt","mps"); 00074 CglKnapsackCover kccg; 00075 int nCols=siP->getNumCols(); 00076 00077 // Test the siP methods for detecting 00078 // variable type 00079 int numCont=0, numBinary=0, numIntNonBinary=0, numInt=0; 00080 for (int thisCol=0; thisCol<nCols; thisCol++) { 00081 if ( siP->isContinuous(thisCol) ) numCont++; 00082 if ( siP->isBinary(thisCol) ) numBinary++; 00083 if ( siP->isIntegerNonBinary(thisCol) ) numIntNonBinary++; 00084 if ( siP->isInteger(thisCol) ) numInt++; 00085 } 00086 assert(numCont==3); 00087 assert(numBinary==2); 00088 assert(numIntNonBinary==1); 00089 assert(numInt==3); 00090 00091 00092 // Test initializeCutGenerator 00093 siP->initialSolve(); 00094 assert(xstar !=NULL); 00095 for (i=0; i<nCols; i++){ 00096 assert(complement[i]==0); 00097 } 00098 int nRows=siP->getNumRows(); 00099 for (i=0; i<nRows; i++){ 00100 int vectorsize = siP->getMatrixByRow()->vectorSize(i); 00101 assert(vectorsize==2); 00102 } 00103 00104 kccg.cleanUpCutGenerator(complement,xstar); 00105 delete siP; 00106 } 00107 */ 00108 00109 // Testcase /u/rlh/osl2/mps/tp3.mps 00110 // Models has 3 cols, 3 rows 00111 // Row 0 yields a knapsack, others do not. 00112 { 00113 // setup 00114 OsiSolverInterface * siP = baseSiP->clone(); 00115 std::string fn(mpsDir+"tp3"); 00116 siP->readMps(fn.c_str(),"mps"); 00117 // All integer variables should be binary. 00118 // Assert that this is true. 00119 for ( i = 0; i < siP->getNumCols(); i++ ) 00120 if ( siP->isInteger(i) ) 00121 assert(siP->getColUpper()[i]==1.0 && siP->isBinary(i)); 00122 OsiCuts cs; 00123 CoinPackedVector krow; 00124 double b=0; 00125 int nCols=siP->getNumCols(); 00126 int * complement=new int [nCols]; 00127 double * xstar=new double [nCols]; 00128 00129 CglKnapsackCover kccg; 00130 00131 // solve LP relaxation 00132 // a "must" before calling initialization 00133 siP->initialSolve(); 00134 assert( eq(siP->getObjValue(), 97.185) ); 00135 double mycs[] = {.627, .667558333333, .038}; 00136 siP->setColSolution(mycs); 00137 const double *colsol = siP->getColSolution(); 00138 int k; 00139 for (k=0; k<nCols; k++){ 00140 xstar[k]=colsol[k]; 00141 complement[k]=0; 00142 } 00143 00144 // test deriveAKnapsack 00145 int rind = ( siP->getRowSense()[0] == 'N' ) ? 1 : 0; 00146 int deriveaknap = kccg.deriveAKnapsack(*siP, cs, krow,b,complement,xstar,rind,siP->getMatrixByRow()->getVector(rind)); 00147 assert(deriveaknap ==1); 00148 assert(complement[0]==0); 00149 assert(complement[1]==1); 00150 assert(complement[2]==1); 00151 int inx[3] = {0,1,2}; 00152 double el[3] = {161, 120, 68}; 00153 CoinPackedVector r; 00154 r.setVector(3,inx,el); 00155 assert (krow == r); 00156 assert (b == 183.0); 00157 00158 00159 // test findGreedyCover 00160 CoinPackedVector cover,remainder; 00161 #if 0 00162 int findgreedy = kccg.findGreedyCover( 0, krow, b, xstar, cover, remainder ); 00163 assert( findgreedy == 1 ); 00164 int coveri = cover.getNumElements(); 00165 assert( cover.getNumElements() == 2); 00166 coveri = cover.getIndices()[0]; 00167 assert( cover.getIndices()[0] == 0); 00168 assert( cover.getIndices()[1] == 1); 00169 assert( cover.getElements()[0] == 161.0); 00170 assert( cover.getElements()[1] == 120.0); 00171 assert( remainder.getNumElements() == 1); 00172 assert( remainder.getIndices()[0] == 2); 00173 assert( remainder.getElements()[0] == 68.0); 00174 00175 // test liftCoverCut 00176 CoinPackedVector cut; 00177 double * rowupper = ekk_rowupper(model); 00178 double cutRhs = cover.getNumElements() - 1.0; 00179 kccg.liftCoverCut(b, krow.getNumElements(), 00180 cover, remainder, 00181 cut); 00182 assert ( cut.getNumElements() == 3 ); 00183 assert ( cut.getIndices()[0] == 0 ); 00184 assert ( cut.getIndices()[1] == 1 ); 00185 assert ( cut.getIndices()[2] == 2 ); 00186 assert( cut.getElements()[0] == 1 ); 00187 assert( cut.getElements()[1] == 1 ); 00188 assert( eq(cut.getElements()[2], 0.087719) ); 00189 00190 // test liftAndUncomplementAndAdd 00191 OsiCuts cuts; 00192 kccg.liftAndUncomplementAndAdd(*siP.getRowUpper()[0],krow,b,complement,0, 00193 cover,remainder,cuts); 00194 int sizerowcuts = cuts.sizeRowCuts(); 00195 assert ( sizerowcuts== 1 ); 00196 OsiRowCut testRowCut = cuts.rowCut(0); 00197 CoinPackedVector testRowPV = testRowCut.row(); 00198 OsiRowCut sampleRowCut; 00199 const int sampleSize = 3; 00200 int sampleCols[sampleSize]={0,1,2}; 00201 double sampleElems[sampleSize]={1.0,-1.0,-0.087719}; 00202 sampleRowCut.setRow(sampleSize,sampleCols,sampleElems); 00203 sampleRowCut.setLb(-DBL_MAX); 00204 sampleRowCut.setUb(-0.087719); 00205 bool equiv = testRowPV.equivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05) ); 00206 assert ( equiv ); 00207 #endif 00208 00209 // test find PseudoJohnAndEllisCover 00210 cover.setVector(0,NULL, NULL); 00211 remainder.setVector(0,NULL,NULL); 00212 00213 rind = ( siP->getRowSense()[0] == 'N' ) ? 1 : 0; 00214 int findPJE = kccg.findPseudoJohnAndEllisCover( rind, krow, 00215 b, xstar, cover, remainder ); 00216 assert( findPJE == 1 ); 00217 assert ( cover.getIndices()[0] == 0 ); 00218 assert ( cover.getIndices()[1] == 2 ); 00219 assert ( cover.getElements()[0] == 161 ); 00220 assert ( cover.getElements()[1] == 68 ); 00221 assert ( remainder.getIndices()[0] == 1 ); 00222 assert ( remainder.getElements()[0] == 120 ); 00223 OsiCuts cuts; 00224 kccg.liftAndUncomplementAndAdd((*siP).getRowUpper()[rind],krow,b, complement, rind, 00225 cover,remainder,cuts); 00226 assert (cuts.sizeRowCuts() == 1 ); 00227 00228 OsiRowCut testRowCut = cuts.rowCut(0); 00229 CoinPackedVector testRowPV = testRowCut.row(); 00230 00231 00232 const int sampleSize = 3; 00233 int sampleCols[sampleSize]={0,1,2}; 00234 double sampleElems[sampleSize]={1.0, -1.0, -1.0}; 00235 OsiRowCut sampleRowCut; 00236 sampleRowCut.setRow(sampleSize,sampleCols,sampleElems); 00237 sampleRowCut.setLb(-DBL_MAX); 00238 sampleRowCut.setUb(-1.0); 00239 00240 // test for 'close enough' 00241 assert( testRowPV.isEquivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05) ) ); 00242 // Reset complement & test next row 00243 for (i=0; i<nCols; i++){ 00244 complement[i]=0; 00245 } 00246 00247 rind++; 00248 deriveaknap = kccg.deriveAKnapsack(*siP,cuts,krow,b,complement,xstar,rind,siP->getMatrixByRow()->getVector(rind)); 00249 assert(deriveaknap==0); 00250 00251 // Reset complement & test next row 00252 for (i=0; i<nCols; i++){ 00253 complement[i]=0; 00254 } 00255 deriveaknap = kccg.deriveAKnapsack(*siP,cuts,krow,b,complement,xstar,2, 00256 siP->getMatrixByRow()->getVector(2)); 00257 assert(deriveaknap == 0); 00258 00259 // Clean up 00260 delete [] complement; 00261 delete [] xstar; 00262 00263 delete siP; 00264 } 00265 00266 #if 0 00267 // Testcase /u/rlh/osl2/mps/tp4.mps 00268 // Models has 6 cols, 1 knapsack row and 00269 // 3 rows explicily bounding variables 00270 // Row 0 yields a knapsack cover cut 00271 // using findGreedyCover which moves the 00272 // LP objective function value. 00273 { 00274 // Setup 00275 EKKContext * env=ekk_initializeContext(); 00276 EKKModel * model = ekk_newModel(env,""); 00277 OsiSolverInterface si(model); 00278 ekk_importModel(model, "tp4.mps"); 00279 CglKnapsackCover kccg; 00280 kccg.ekk_validateIntType(si); 00281 00282 // Solve the LP relaxation of the model and 00283 // print out ofv for sake of comparison 00284 ekk_allSlackBasis(model); 00285 ekk_crash(model,1); 00286 ekk_primalSimplex(model,1); 00287 double lpRelaxBefore=ekk_getRobjvalue(model); 00288 #ifdef CGL_DEBUG 00289 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00290 #endif 00291 00292 // Determine if lp sol is ip optimal 00293 // Note: no ekk_function to do this 00294 int nCols=ekk_getInumcols(model); 00295 double * optLpSol = ekk_colsol(model); 00296 int ipOpt = 1; 00297 i=0; 00298 while (i++<nCols && ipOpt){ 00299 if(optLpSol[i] < 1.0-1.0e-08 && optLpSol[i]> 1.0e-08) ipOpt = 0; 00300 } 00301 00302 if (ipOpt){ 00303 #ifdef CGL_DEBUG 00304 printf("Lp solution is within ip optimality tolerance\n"); 00305 #endif 00306 } 00307 else { 00308 OsiSolverInterface iModel(model); 00309 OsiCuts cuts; 00310 00311 // Test generateCuts method 00312 kccg.generateCuts(iModel,cuts); 00313 OsiSolverInterface::ApplyCutsReturnCode rc = iModel.applyCuts(cuts); 00314 00315 ekk_mergeBlocks(model,1); 00316 ekk_dualSimplex(model); 00317 double lpRelaxAfter=ekk_getRobjvalue(model); 00318 #ifdef CGL_DEBUG 00319 printf("\n\nFinal LP min=%f\n",lpRelaxAfter); 00320 #endif 00321 assert( lpRelaxBefore < lpRelaxAfter ); 00322 00323 // This may need to be updated as other 00324 // minimal cover finders are added 00325 assert( cuts.sizeRowCuts() == 1 ); 00326 OsiRowCut testRowCut = cuts.rowCut(0); 00327 CoinPackedVector testRowPV = testRowCut.row(); 00328 00329 OsiRowCut sampleRowCut; 00330 const int sampleSize = 6; 00331 int sampleCols[sampleSize]={0,1,2,3,4,5}; 00332 double sampleElems[sampleSize]={1.0,1.0,1.0,1.0,0.5, 2.0}; 00333 sampleRowCut.setRow(sampleSize,sampleCols,sampleElems); 00334 sampleRowCut.setLb(-DBL_MAX); 00335 sampleRowCut.setUb(3.0); 00336 bool equiv = testRowPV.equivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05) ); 00337 assert( testRowPV.equivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05) ) ); 00338 } 00339 00340 // Exit out of OSL 00341 ekk_deleteModel(model); 00342 ekk_endContext(env); 00343 00344 } 00345 #endif 00346 00347 00348 // Testcase /u/rlh/osl2/mps/tp5.mps 00349 // Models has 6 cols, 1 knapsack row and 00350 // 3 rows explicily bounding variables 00351 // Row 0 yields a knapsack cover cut 00352 // using findGreedyCover which moves the 00353 // LP objective function value. 00354 { 00355 // Setup 00356 OsiSolverInterface * siP = baseSiP->clone(); 00357 std::string fn(mpsDir+"tp5"); 00358 siP->readMps(fn.c_str(),"mps"); 00359 // All integer variables should be binary. 00360 // Assert that this is true. 00361 for ( i = 0; i < siP->getNumCols(); i++ ) 00362 if ( siP->isInteger(i) ) 00363 assert(siP->getColUpper()[i]==1.0 && siP->isBinary(i)); 00364 CglKnapsackCover kccg; 00365 00366 // Solve the LP relaxation of the model and 00367 // print out ofv for sake of comparison 00368 siP->initialSolve(); 00369 double lpRelaxBefore=siP->getObjValue(); 00370 assert( eq(lpRelaxBefore, -51.66666666667) ); 00371 double mycs[] = {.8999999999, .899999999999, .89999999999, 1.110223e-16, .5166666666667, 0}; 00372 siP->setColSolution(mycs); 00373 #ifdef CGL_DEBUG 00374 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00375 #endif 00376 00377 // Determine if lp sol is 0/1 optimal 00378 int nCols=siP->getNumCols(); 00379 const double * optLpSol = siP->getColSolution(); 00380 bool ipOpt = true; 00381 i=0; 00382 while (i++<nCols && ipOpt){ 00383 if(optLpSol[i] > kccg.epsilon_ && optLpSol[i] < kccg.onetol_) ipOpt = false; 00384 } 00385 00386 if (ipOpt){ 00387 #ifdef CGL_DEBUG 00388 printf("Lp solution is within ip optimality tolerance\n"); 00389 #endif 00390 } 00391 else { 00392 // set up 00393 OsiCuts cuts; 00394 CoinPackedVector krow; 00395 double b=0.0; 00396 int * complement=new int[nCols]; 00397 double * xstar=new double[nCols]; 00398 // initialize cut generator 00399 const double *colsol = siP->getColSolution(); 00400 for (i=0; i<nCols; i++){ 00401 xstar[i]=colsol[i]; 00402 complement[i]=0; 00403 } 00404 int row = ( siP->getRowSense()[0] == 'N' ) ? 1 : 0; 00405 // transform row into canonical knapsack form 00406 if (kccg.deriveAKnapsack(*siP, cuts, krow, b, complement, xstar, row,siP->getMatrixByRow()->getVector(row))){ 00407 CoinPackedVector cover, remainder; 00408 // apply greedy logic to detect violated minimal cover inequalities 00409 if (kccg.findGreedyCover(row, krow, b, xstar, cover, remainder) == 1){ 00410 // lift, uncomplements, and add cut to cut set 00411 kccg.liftAndUncomplementAndAdd((*siP).getRowUpper()[row],krow, b, complement, row, cover, remainder, cuts); 00412 } 00413 // reset optimal column solution (xstar) information in OSL 00414 const double * rowupper = siP->getRowUpper(); 00415 int k; 00416 if (fabs(b-rowupper[row]) > 1.0e-05) { 00417 for(k=0; k<krow.getNumElements(); k++) { 00418 if (complement[krow.getIndices()[k]]){ 00419 xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]]; 00420 complement[krow.getIndices()[k]]=0; 00421 } 00422 } 00423 } 00424 // clean up 00425 delete [] complement; 00426 delete [] xstar; 00427 } 00428 // apply the cuts 00429 OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts); 00430 00431 siP->resolve(); 00432 double lpRelaxAfter=siP->getObjValue(); 00433 assert( eq(lpRelaxAfter, -30.0) ); 00434 #ifdef CGL_DEBUG 00435 printf("\n\nFinal LP min=%f\n",lpRelaxAfter); 00436 #endif 00437 // test that expected cut was detected 00438 assert( lpRelaxBefore < lpRelaxAfter ); 00439 assert( cuts.sizeRowCuts() == 1 ); 00440 OsiRowCut testRowCut = cuts.rowCut(0); 00441 CoinPackedVector testRowPV = testRowCut.row(); 00442 OsiRowCut sampleRowCut; 00443 const int sampleSize = 6; 00444 int sampleCols[sampleSize]={0,1,2,3,4,5}; 00445 double sampleElems[sampleSize]={1.0,1.0,1.0,0.25,1.0,2.0}; 00446 sampleRowCut.setRow(sampleSize,sampleCols,sampleElems); 00447 sampleRowCut.setLb(-DBL_MAX); 00448 sampleRowCut.setUb(3.0); 00449 assert(testRowPV.isEquivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05))); 00450 } 00451 00452 delete siP; 00453 } 00454 00455 00456 // Testcase /u/rlh/osl2/mps/p0033 00457 // Miplib3 problem p0033 00458 // Test that no cuts chop off the optimal solution 00459 { 00460 // Setup 00461 OsiSolverInterface * siP = baseSiP->clone(); 00462 std::string fn(mpsDir+"p0033"); 00463 siP->readMps(fn.c_str(),"mps"); 00464 // All integer variables should be binary. 00465 // Assert that this is true. 00466 for ( i = 0; i < siP->getNumCols(); i++ ) 00467 if ( siP->isInteger(i) ) 00468 assert(siP->getColUpper()[i]==1.0 && siP->isBinary(i)); 00469 int nCols=siP->getNumCols(); 00470 CglKnapsackCover kccg; 00471 00472 // Solve the LP relaxation of the model and 00473 // print out ofv for sake of comparison 00474 siP->initialSolve(); 00475 double lpRelaxBefore=siP->getObjValue(); 00476 assert( eq(lpRelaxBefore, 2520.5717391304347) ); 00477 double mycs[] = {0, 1, 0, 0, -2.0837010502455788e-19, 1, 0, 0, 1, 00478 0.021739130434782594, 0.35652173913043478, 00479 -6.7220534694101275e-18, 5.3125906451789717e-18, 00480 1, 0, 1.9298798670241979e-17, 0, 0, 0, 00481 7.8875708048320448e-18, 0.5, 0, 00482 0.85999999999999999, 1, 1, 0.57999999999999996, 00483 1, 0, 1, 0, 0.25, 0, 0.67500000000000004}; 00484 siP->setColSolution(mycs); 00485 #ifdef CGL_DEBUG 00486 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00487 #endif 00488 00489 OsiCuts cuts; 00490 00491 // Test generateCuts method 00492 kccg.generateCuts(*siP,cuts); 00493 OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts); 00494 00495 siP->resolve(); 00496 double lpRelaxAfter=siP->getObjValue(); 00497 assert( eq(lpRelaxAfter, 2829.0597826086955) ); 00498 #ifdef CGL_DEBUG 00499 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00500 printf("\n\nFinal LP min=%f\n",lpRelaxAfter); 00501 #endif 00502 assert( lpRelaxBefore < lpRelaxAfter ); 00503 00504 // the CoinPackedVector p0033 is the optimal 00505 // IP solution to the miplib problem p0033 00506 int objIndices[14] = { 00507 0, 6, 7, 9, 13, 17, 18, 00508 22, 24, 25, 26, 27, 28, 29 }; 00509 CoinPackedVector p0033(14,objIndices,1.0); 00510 00511 // Sanity check 00512 const double * objective=siP->getObjCoefficients(); 00513 double ofv =0 ; 00514 int r; 00515 for (r=0; r<nCols; r++){ 00516 ofv=ofv + p0033[r]*objective[r]; 00517 } 00518 CoinRelFltEq eq; 00519 assert( eq(ofv,3089.0) ); 00520 00521 int nRowCuts = cuts.sizeRowCuts(); 00522 OsiRowCut rcut; 00523 CoinPackedVector rpv; 00524 for (i=0; i<nRowCuts; i++){ 00525 rcut = cuts.rowCut(i); 00526 rpv = rcut.row(); 00527 double p0033Sum = (rpv*p0033).sum(); 00528 assert (p0033Sum <= rcut.ub() ); 00529 } 00530 00531 delete siP; 00532 } 00533 00534 // if a debug file is there then look at it 00535 { 00536 FILE * fp = fopen("knapsack.debug","r"); 00537 if (fp) { 00538 int ncol,nel; 00539 double up; 00540 fscanf(fp,"%d %d %lg",&ncol,&nel,&up); 00541 printf("%d columns, %d elements, upper %g\n",ncol,nel,up); 00542 double * sol1 = new double[nel]; 00543 double * el1 = new double[nel]; 00544 int * col1 = new int[nel]; 00545 int * start = new int[ncol+1]; 00546 memset(start,0,ncol*sizeof(int)); 00547 int * row = new int[nel]; 00548 int i; 00549 for (i=0;i<nel;i++) { 00550 fscanf(fp,"%d %lg %lg",col1+i,el1+i,sol1+i); 00551 printf("[%d, e=%g, v=%g] ",col1[i],el1[i],sol1[i]); 00552 start[col1[i]]=1; 00553 row[i]=0; 00554 } 00555 printf("\n"); 00556 // Setup 00557 OsiSolverInterface * siP = baseSiP->clone(); 00558 00559 double lo=-1.0e30; 00560 double * upper = new double[ncol]; 00561 start[ncol]=nel; 00562 int last=0; 00563 for (i=0;i<ncol;i++) { 00564 upper[i]=1.0; 00565 int marked=start[i]; 00566 start[i]=last; 00567 if (marked) 00568 last++; 00569 } 00570 siP->loadProblem(ncol,1,start,row,el1,NULL,upper,NULL,&lo,&up); 00571 // use upper for solution 00572 memset(upper,0,ncol*sizeof(double)); 00573 for (i=0;i<nel;i++) { 00574 int icol=col1[i]; 00575 upper[icol]=sol1[i]; 00576 siP->setInteger(icol); 00577 } 00578 siP->setColSolution(upper); 00579 delete [] sol1; 00580 delete [] el1; 00581 delete [] col1; 00582 delete [] start; 00583 delete [] row; 00584 delete [] upper; 00585 CglKnapsackCover kccg; 00586 00587 OsiCuts cuts; 00588 00589 // Test generateCuts method 00590 kccg.generateCuts(*siP,cuts); 00591 // print out and compare to known cuts 00592 int numberCuts = cuts.sizeRowCuts(); 00593 if (numberCuts) { 00594 for (i=0;i<numberCuts;i++) { 00595 OsiRowCut * thisCut = cuts.rowCutPtr(i); 00596 int n=thisCut->row().getNumElements(); 00597 printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(), 00598 thisCut->ub()); 00599 int j; 00600 const int * index = thisCut->row().getIndices(); 00601 const double * element = thisCut->row().getElements(); 00602 for (j=0;j<n;j++) { 00603 printf(" (%d,%g)",index[j],element[j]); 00604 } 00605 printf("\n"); 00606 } 00607 } 00608 fclose(fp); 00609 } 00610 } 00611 00612 // Testcase /u/rlh/osl2/mps/p0201 00613 // Miplib3 problem p0282 00614 // Test that no cuts chop off the optimal ip solution 00615 { 00616 // Setup 00617 OsiSolverInterface * siP = baseSiP->clone(); 00618 std::string fn(mpsDir+"p0201"); 00619 siP->readMps(fn.c_str(),"mps"); 00620 // All integer variables should be binary. 00621 // Assert that this is true. 00622 for ( i = 0; i < siP->getNumCols(); i++ ) 00623 if ( siP->isInteger(i) ) 00624 assert(siP->getColUpper()[i]==1.0 && siP->isBinary(i)); 00625 00626 const int nCols=siP->getNumCols(); 00627 CglKnapsackCover kccg; 00628 00629 // Solve the LP relaxation of the model and 00630 // print out ofv for sake of comparisn 00631 siP->initialSolve(); 00632 double lpRelaxBefore=siP->getObjValue(); 00633 assert( eq(lpRelaxBefore, 6875.) ); 00634 double mycs[] = 00635 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 00636 0, 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5, 00637 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5, 0, 0, 00638 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5, 0, 0, 0, 0.5, 00639 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5, 0, 0, 0, 0.5, 0, 0, 00640 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 00641 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 00642 0, 0, 0, 0, 0, 0, 1, 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 00643 0, 0, 0, 0, 1, 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 00644 0, 0, 1, 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 00645 1}; 00646 siP->setColSolution(mycs); 00647 #ifdef CGL_DEBUG 00648 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00649 #endif 00650 00651 OsiCuts cuts; 00652 00653 // Test generateCuts method 00654 kccg.generateCuts(*siP,cuts); 00655 OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts); 00656 00657 siP->resolve(); 00658 double lpRelaxAfter=siP->getObjValue(); 00659 assert( eq(lpRelaxAfter, 7125) ); 00660 #ifdef CGL_DEBUG 00661 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00662 printf("\n\nFinal LP min=%f\n",lpRelaxAfter); 00663 #endif 00664 assert( lpRelaxBefore < lpRelaxAfter ); 00665 00666 // Optimal IP solution to p0201 00667 int objIndices[22] = { 8, 10, 21, 38, 39, 56, 00668 60, 74, 79, 92, 94, 110, 111, 128, 132, 146, 00669 151,164, 166, 182,183, 200 }; 00670 CoinPackedVector p0201(22,objIndices,1.0); 00671 00672 // Sanity check 00673 const double * objective=siP->getObjCoefficients(); 00674 double ofv =0 ; 00675 int r; 00676 for (r=0; r<nCols; r++){ 00677 ofv=ofv + p0201[r]*objective[r]; 00678 } 00679 CoinRelFltEq eq; 00680 assert( eq(ofv,7615.0) ); 00681 //printf("p0201 optimal ofv = %g\n",ofv); 00682 00683 int nRowCuts = cuts.sizeRowCuts(); 00684 OsiRowCut rcut; 00685 CoinPackedVector rpv; 00686 for (i=0; i<nRowCuts; i++){ 00687 rcut = cuts.rowCut(i); 00688 rpv = rcut.row(); 00689 double p0201Sum = (rpv*p0201).sum(); 00690 assert (p0201Sum <= rcut.ub() ); 00691 } 00692 00693 delete siP; 00694 } 00695 00696 00697 // see if I get the same covers that N&W get 00698 { 00699 OsiSolverInterface * siP=baseSiP->clone(); 00700 std::string fn(mpsDir+"nw460"); 00701 siP->readMps(fn.c_str(),"mps"); 00702 // All integer variables should be binary. 00703 // Assert that this is true. 00704 for ( i = 0; i < siP->getNumCols(); i++ ) 00705 if ( siP->isInteger(i) ) 00706 assert(siP->getColUpper()[i]==1.0 && siP->isBinary(i)); 00707 CglKnapsackCover kccg; 00708 00709 // Solve the LP relaxation of the model and 00710 // print out ofv for sake of comparison 00711 siP->initialSolve(); 00712 double lpRelaxBefore=siP->getObjValue(); 00713 assert( eq(lpRelaxBefore, -225.68951787852194) ); 00714 double mycs[] = {0.7099213482046447, 0, 0.34185802225477174, 1, 1, 0, 1, 1, 0}; 00715 siP->setColSolution(mycs); 00716 00717 OsiCuts cuts; 00718 00719 // Test generateCuts method 00720 kccg.generateCuts(*siP,cuts); 00721 OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts); 00722 00723 siP->resolve(); 00724 double lpRelaxAfter=siP->getObjValue(); 00725 assert( eq(lpRelaxAfter, -176) ); 00726 #ifdef CGL_DEBUG 00727 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00728 printf("\n\nFinal LP min=%f\n",lpRelaxAfter); 00729 #endif 00730 #ifdef MJS 00731 assert( lpRelaxBefore < lpRelaxAfter ); 00732 #endif 00733 00734 int nRowCuts = cuts.sizeRowCuts(); 00735 OsiRowCut rcut; 00736 CoinPackedVector rpv; 00737 for (i=0; i<nRowCuts; i++){ 00738 rcut = cuts.rowCut(i); 00739 rpv = rcut.row(); 00740 int j; 00741 printf("Row cut number %i has rhs = %g\n",i,rcut.ub()); 00742 for (j=0; j<rpv.getNumElements(); j++){ 00743 printf("index %i, element %g\n", rpv.getIndices()[j], rpv.getElements()[j]); 00744 } 00745 printf("\n"); 00746 } 00747 delete siP; 00748 } 00749 00750 // Debugging: try "exmip1.mps" 00751 { 00752 // Setup 00753 OsiSolverInterface * siP = baseSiP->clone(); 00754 std::string fn(mpsDir+"exmip1"); 00755 siP->readMps(fn.c_str(),"mps"); 00756 // All integer variables should be binary. 00757 // Assert that this is true. 00758 for ( i = 0; i < siP->getNumCols(); i++ ) 00759 if ( siP->isInteger(i) ) 00760 assert(siP->getColUpper()[i]==1.0 && siP->isBinary(i)); 00761 CglKnapsackCover kccg; 00762 00763 // Solve the LP relaxation of the model and 00764 // print out ofv for sake of comparison 00765 siP->initialSolve(); 00766 double lpRelaxBefore=siP->getObjValue(); 00767 assert( eq(lpRelaxBefore, 3.2368421052631575) ); 00768 double mycs[] = {2.5, 0, 0, 0.6428571428571429, 0.5, 4, 0, 0.26315789473684253}; 00769 siP->setColSolution(mycs); 00770 // Test generateCuts method 00771 OsiCuts cuts; 00772 kccg.generateCuts(*siP,cuts); 00773 OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts); 00774 00775 siP->resolve(); 00776 double lpRelaxAfter=siP->getObjValue(); 00777 assert( eq(lpRelaxAfter, 3.2368421052631575) ); 00778 #ifdef CGL_DEBUG 00779 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00780 printf("\n\nFinal LP min=%f\n",lpRelaxAfter); 00781 #endif 00782 assert( lpRelaxBefore <= lpRelaxAfter ); 00783 00784 delete siP; 00785 } 00786 00787 #ifdef CGL_DEBUG 00788 // See what findLPMostViolatedMinCover for knapsack with 2 elements does 00789 { 00790 int nCols = 2; 00791 int row = 1; 00792 CoinPackedVector krow; 00793 double e[2] = {5,10}; 00794 int ii[2] = {0,1}; 00795 krow.setVector(nCols,ii,e); 00796 double b=11; 00797 double xstar[2] = {.2,.9}; 00798 CoinPackedVector cover; 00799 CoinPackedVector remainder; 00800 CglKnapsackCover kccg; 00801 kccg.findLPMostViolatedMinCover(nCols, row, krow, b, xstar, cover, remainder); 00802 printf("num in cover = %i\n",cover.getNumElements()); 00803 int j; 00804 for (j=0; j<cover.getNumElements(); j++){ 00805 printf(" index %i element % g\n", cover.getIndices()[j], cover.getElements()[j]); 00806 } 00807 } 00808 #endif 00809 00810 #ifdef CGL_DEBUG 00811 // see what findLPMostViolatedMinCover does 00812 { 00813 int nCols = 5; 00814 int row = 1; 00815 CoinPackedVector krow; 00816 double e[5] = {1,1,1,1,10}; 00817 int ii[5] = {0,1,2,3,4}; 00818 krow.setVector(nCols,ii,e); 00819 double b=11; 00820 double xstar[5] = {.9,.9,1,1,.1}; 00821 CoinPackedVector cover; 00822 CoinPackedVector remainder; 00823 CglKnapsackCover kccg; 00824 kccg.findLPMostViolatedMinCover(nCols, row, krow, b, xstar, cover, remainder); 00825 printf("num in cover = %i\n",cover.getNumElements()); 00826 int j; 00827 for (j=0; j<cover.getNumElements(); j++){ 00828 printf(" index %i element % g\n", cover.getIndices()[j], cover.getElements()[j]); 00829 } 00830 } 00831 #endif 00832 00833 } |
|
which rows to look at. If specified, only these rows will be considered for generating knapsack covers. Otherwise all rows will be tried Definition at line 228 of file CglKnapsackCover.hpp. Referenced by CglKnapsackCover(), generateCuts(), operator=(), and setTestedRowIndices(). |