source: trunk/Cgl/src/CglKnapsackCover/CglKnapsackCover.cpp @ 508

Last change on this file since 508 was 508, checked in by jpfasano, 14 years ago

Fix invalid null pointer when running "cbc -miplib".
See https://projects.coin-or.org/Cbc/ticket/30

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 90.3 KB
Line 
1// Copyright (C) 2000, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cstdlib>
8#include <cstdio>
9#include <cmath>
10#include <cassert>
11#include <cfloat>
12#include <iostream>
13
14#include "CoinHelperFunctions.hpp"
15#include "CglKnapsackCover.hpp"
16#include "CoinPackedVector.hpp"
17#include "CoinSort.hpp"
18#include "CoinPackedMatrix.hpp"
19#include "OsiRowCutDebugger.hpp"
20//#define PRINT_DEBUG
21//#define CGL_DEBUG
22//-----------------------------------------------------------------------------
23// Generate knapsack cover cuts
24//-------------------------------------------------------------------
25void CglKnapsackCover::generateCuts(const OsiSolverInterface& si, OsiCuts& cs,
26                                    const CglTreeInfo info) const
27{
28  // Get basic problem information
29  int nRows=si.getNumRows(); 
30  int nCols=si.getNumCols(); 
31
32  // Create working space for "canonical" knapsack inequality
33  // - krow will contain the coefficients and indices of the
34  // (potentially complemented) variables in the knapsack inequality.
35  // - b is the rhs of knapsack inequality.
36  // - complement[i] is 1 if the index i in krow refers to the complement
37  // of the variable, and 0 otherwise.
38  CoinPackedVector krow; 
39  double b=0.0;
40  int * complement= new int[nCols];
41   
42  // Create a local copy of the column solution (colsol), call it xstar, and
43  // inititalize it.
44  // Assumes the lp-relaxation has been solved, and the solver interface
45  // has a meaningful colsol.
46  double * xstar= new double[nCols]; 
47
48  // To allow for vub knapsacks
49  int * thisColumnIndex = new int [nCols];
50  double * thisElement = new double[nCols];
51  int * back = new int[nCols];
52 
53  const double *colsol = si.getColSolution(); 
54  int k; 
55  // For each row point to vub variable
56  // -1 if no vub
57  // -2 if can skip row for knapsacks
58
59  int * vub = new int [nRows];
60
61  // Now vubValue are for positive coefficients and vlbValue for negative
62  // when L row
63  // For each column point to vub row
64  int * vubRow = new int [nCols];
65  double * vubValue = new double [nRows];
66
67  // For each column point to vlb row
68  int * vlbRow = new int [nCols];
69  double * vlbValue = new double [nRows];
70
71  // Take out all fixed
72  double * effectiveUpper = new double [nRows];
73  double * effectiveLower = new double [nRows];
74  const double * colUpper = si.getColUpper();
75  const double * colLower = si.getColLower();
76  for (k=0; k<nCols; k++){
77    back[k]=-1;
78    xstar[k]=colsol[k];
79    complement[k]=0;
80    vubRow[k]=-1;
81    vlbRow[k]=-1;
82    if (si.isBinary(k)) {
83      if (si.isFreeBinary(k)) {
84        vubRow[k]=-2;
85        vlbRow[k]=-2;
86      } else {
87        vubRow[k]=-10;
88        vlbRow[k]=-10;
89      }
90    } else if (colUpper[k]==colLower[k]) {
91      vubRow[k]=-10; // fixed
92      vlbRow[k]=-10; // fixed
93    }
94  }
95
96  int rowIndex;
97  int numberVub=0;
98
99  const CoinPackedMatrix * matrixByRow = si.getMatrixByRow();
100  const double * elementByRow = matrixByRow->getElements();
101  const int * column = matrixByRow->getIndices();
102  const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
103  const int * rowLength = matrixByRow->getVectorLengths();
104  const double * rowUpper = si.getRowUpper();
105  const double * rowLower = si.getRowLower();
106
107  // Scan all rows looking for possibles
108
109  for (rowIndex=0;rowIndex<nRows;rowIndex++) {
110    vub[rowIndex]=-1;
111    int start = rowStart[rowIndex];
112    int end = start + rowLength[rowIndex];
113    double upRhs = rowUpper[rowIndex]; 
114    double loRhs = rowLower[rowIndex]; 
115    double multiplier=0.0;
116    if (upRhs>1.0e20)
117      multiplier=-1.0;
118    else if (loRhs<-1.0e20)
119      multiplier=1.0;
120    int numberContinuous=0;
121    int numberBinary=0;
122    int iCont=-1;
123    double sum = 0.0;
124    double valueContinuous=0.0;
125    double valueBinary=0.0;
126    int iBinary=-1;
127    int j;
128    for (j=start;j<end;j++) {
129      int iColumn=column[j];
130      double value = elementByRow[j];
131      if (colUpper[iColumn]>colLower[iColumn]) {
132        sum += colsol[iColumn]*value;
133        if (vubRow[iColumn]==-2&&value*multiplier>0.0) {
134          // binary
135          numberBinary++;
136          valueBinary=value;
137          iBinary=iColumn;
138        } else if (vlbRow[iColumn]==-2&&value*multiplier<0.0) {
139          // binary
140          numberBinary++;
141          valueBinary=value;
142          iBinary=iColumn;
143        } else if (vubRow[iColumn]==-1) {
144          // only use if not at bound
145          //      if (colsol[iColumn]<colUpper[iColumn]-1.0e-6&&
146          //  colsol[iColumn]>colLower[iColumn]+1.0e-6) {
147            // possible
148            iCont=iColumn;
149            numberContinuous++;
150            valueContinuous=value;
151            //} else {
152            //// ** needs more thought
153            //numberContinuous ++;
154            //iCont=-1;
155            //}
156        } else {
157          // ** needs more thought
158          numberContinuous ++;
159          iCont=-1;
160          //if (colsol[iColumn]<colUpper[iColumn]-1.0e-6&&
161          //  colsol[iColumn]>colLower[iColumn]+1.0e-6) {
162          //// already assigned
163          //numberContinuous ++;
164          //iCont=-1;
165          //}
166        }
167      } else {
168        // fixed
169        upRhs -= colLower[iColumn]*value;
170        loRhs -= colLower[iColumn]*value;
171      }
172    }
173    // see if binding
174    effectiveUpper[rowIndex] = upRhs;
175    effectiveLower[rowIndex] = loRhs;
176    bool possible = false;
177    if (fabs(sum-upRhs)<1.0e-5) {
178      possible=true;
179    } else {
180      effectiveUpper[rowIndex]=DBL_MAX;
181    }
182    if (fabs(sum-loRhs)<1.0e-5) {
183      possible=true;
184    } else {
185      effectiveLower[rowIndex]=-DBL_MAX;
186    }
187    if (possible&&numberBinary&&numberBinary+numberContinuous<=maxInKnapsack_) {
188      // binding with binary
189      if(numberContinuous==1&&iCont>=0&&numberBinary==1) {
190        // vub
191#ifdef PRINT_DEBUG
192        printf("vub/vlb (by row %d) %g <= 0-1 %g * %d + %g * %d <= %g\n",
193               rowIndex,effectiveLower[rowIndex],valueBinary,iBinary,
194               valueContinuous,iCont,effectiveUpper[rowIndex]);
195#endif
196        if (multiplier*valueContinuous>0.0) {
197          vubValue[rowIndex] = valueContinuous;
198          vubRow[iCont]=rowIndex;
199        } else {
200          vlbValue[rowIndex] = valueContinuous;
201          vlbRow[iCont]=rowIndex;
202        }
203        vub[rowIndex]=iCont;
204        numberVub++;
205      } else if (numberBinary>1) {
206        // could be knapsack
207        vub[rowIndex]=-1;
208      } else {
209        // no point looking at this row
210        vub[rowIndex]=-2;
211      }
212    } else {
213      if (!possible||numberBinary+numberContinuous>maxInKnapsack_)
214        vub[rowIndex]=-2;      // no point looking at this row
215    }
216  }
217  // Main loop
218  int numCheck = 0;
219  int* toCheck = 0;
220  if (!rowsToCheck_) {
221     toCheck = new int[nRows];
222     CoinIotaN(toCheck, nRows, 0);
223     numCheck = nRows;
224  } else {
225     numCheck = numRowsToCheck_;
226     toCheck = rowsToCheck_;
227  }
228
229  // Set up number of tries for each row
230  int ntry;
231  if (numberVub) 
232    ntry=4;
233  else
234    ntry=2;
235  //ntry=2; // switch off
236  for (int ii=0; ii < numCheck; ++ii){
237     rowIndex = toCheck[ii];
238     if (rowIndex < 0 || rowIndex >= nRows)
239        continue;
240     if (vub[rowIndex]==-2)
241       continue;
242
243#ifdef PRINT_DEBUG
244    std::cout << "CGL: Processing row " << rowIndex << std::endl;
245#endif
246   
247    // Get a tight row
248    // (want to be able to
249    //  experiment by turning this on and off)
250    //
251    // const double * pi=si.rowprice();
252    // if (fabs(pi[row]) < epsilon_){
253    //  continue;
254    // }
255   
256   
257    //////////////////////////////////////////////////////
258    // Derive a "canonical"  knapsack                  //
259    // inequality (in binary variables)                 //
260    // from the model row in mixed integer variables    //
261    //////////////////////////////////////////////////////
262#ifdef CGL_DEBUG
263    assert(!krow.getNumElements());
264#endif
265    double effectiveRhs[4];
266    double rhs[4];
267    double sign[]={0.0,0.0,-1.0,1.0};
268    bool rowType[] = {false,true,false,true};
269    effectiveRhs[0] = effectiveLower[rowIndex]; 
270    rhs[0]=rowLower[rowIndex];
271    effectiveRhs[2] = effectiveRhs[0];
272    rhs[2]= effectiveRhs[0];
273    effectiveRhs[1] = effectiveUpper[rowIndex]; 
274    rhs[1]=rowUpper[rowIndex];
275    effectiveRhs[3] = effectiveRhs[1];
276    rhs[3]= effectiveRhs[1];
277    int itry;
278#ifdef CGL_DEBUG
279    int kcuts[4];
280    memset(kcuts,0,4*sizeof(int));
281#endif
282    for (itry=0;itry<ntry;itry++) {
283#ifdef CGL_DEBUG
284      int nlast=cs.sizeRowCuts();
285#endif
286      // see if to skip
287      if (fabs(effectiveRhs[itry])>1.0e20)
288        continue;
289      int length = rowLength[rowIndex];
290      memcpy(thisColumnIndex,column+rowStart[rowIndex],length*sizeof(int));
291      memcpy(thisElement,elementByRow+rowStart[rowIndex],
292             length*sizeof(double));
293      b=rhs[itry];
294      if (itry>1) {
295        // see if we would be better off relaxing
296        int i;
297        // mark columns
298        int length2=length; // for new length
299        int numberReplaced=0;
300        for (i=0;i<length;i++) {
301          int iColumn = thisColumnIndex[i];
302          back[thisColumnIndex[i]]=i;
303          if (vubRow[iColumn]==-10) {
304            // fixed - take out
305            thisElement[i]=0.0;
306          }
307        }
308        double dSign = sign[itry];
309        for (i=0;i<length;i++) {
310          int iColumn = thisColumnIndex[i];
311          int iRow=-1;
312          double vubCoefficient=0.0;
313          double thisCoefficient=thisElement[i];
314          int replace = 0;
315          if (vubRow[iColumn]>=0) {
316            iRow = vubRow[iColumn];
317            if (vub[iRow]==iColumn&&iRow!=rowIndex) {
318              vubCoefficient = vubValue[iRow];
319              // break it out - may be able to do better
320              if (dSign*thisCoefficient>0.0) {
321                // we want valid lower bound on continuous
322                if (effectiveLower[iRow]>-1.0e20&&vubCoefficient>0.0) 
323                  replace=-1;
324                else if (effectiveUpper[iRow]<1.0e20&&vubCoefficient<0.0) 
325                  replace=1;
326                // q assert (replace!=-1);
327                // q assert (replace!=1);
328              } else {
329                // we want valid upper bound on continuous
330                if (effectiveLower[iRow]>-1.0e20&&vubCoefficient<0.0) 
331                  replace=-1;
332                else if (effectiveUpper[iRow]<1.0e20&&vubCoefficient>0.0) 
333                  replace=1;
334                //assert (replace!=-1);
335              }
336            }
337          }
338          if (vlbRow[iColumn]>=0) {
339            iRow = vlbRow[iColumn];
340            if (vub[iRow]==iColumn&&iRow!=rowIndex) {
341              vubCoefficient = vlbValue[iRow];
342              // break it out - may be able to do better
343              if (dSign*thisCoefficient>0.0) {
344                // we want valid lower bound on continuous
345                if (effectiveLower[iRow]>-1.0e20&&vubCoefficient>0.0) 
346                  replace=-1;
347                else if (effectiveUpper[iRow]<1.0e20&&vubCoefficient<0.0) 
348                  replace=1;
349                //assert (replace!=1);
350              } else {
351                // we want valid upper bound on continuous
352                if (effectiveLower[iRow]>-1.0e20&&vubCoefficient<0.0) 
353                  replace=-1;
354                else if (effectiveUpper[iRow]<1.0e20&&vubCoefficient>0.0) 
355                  replace=1;
356                //q assert (replace!=-1);
357                //assert (replace!=1);
358              }
359            }
360          }
361          if (replace) {
362            double useRhs=0.0;
363            numberReplaced++;
364            if (replace<0)
365              useRhs = effectiveLower[iRow];
366            else
367              useRhs = effectiveUpper[iRow];
368            // now replace (just using vubRow==-2)
369            // delete continuous
370            thisElement[i]=0.0;
371            double scale = thisCoefficient/vubCoefficient;
372            // modify rhs
373            b -= scale*useRhs;
374            int start = rowStart[iRow];
375            int end = start+rowLength[iRow];
376            int j;
377            for (j=start;j<end;j++) {
378              int iColumn = column[j];
379              if (vubRow[iColumn]==-2) {
380                double change = scale*elementByRow[j];
381                int iBack = back[iColumn];
382                if (iBack<0) {
383                  // element does not exist
384                  back[iColumn]=length2;
385                  thisElement[length2]=-change;
386                  thisColumnIndex[length2++]=iColumn;
387                } else {
388                  // element does exist
389                  thisElement[iBack] -= change;
390                }
391              }
392            }
393          }
394        }
395        if (numberReplaced) {
396          length=0;
397          for (i=0;i<length2;i++) {
398            int iColumn = thisColumnIndex[i];
399            back[iColumn]=-1; // un mark
400            if (thisElement[i]) {
401              thisElement[length]=thisElement[i];
402              thisColumnIndex[length++]=iColumn;
403            }
404          }
405          if (length>maxInKnapsack_)
406            continue; // too long
407        } else {
408          for (i=0;i<length;i++) {
409            int iColumn = thisColumnIndex[i];
410            back[iColumn]=-1; // un mark
411          }
412          continue; // no good
413        }
414      }
415      if (!deriveAKnapsack(si, cs, krow, rowType[itry], b, complement, 
416                           xstar, rowIndex, 
417                           length,thisColumnIndex,thisElement)) {
418       
419        // Reset local data and continue to the next iteration
420        // of the rowIndex-loop
421        for(k=0; k<krow.getNumElements(); k++) {
422          if (complement[krow.getIndices()[k]]){
423            xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]];
424            complement[krow.getIndices()[k]]=0;       
425          }
426        }
427        krow.setVector(0,NULL,NULL);
428        continue;
429      }
430#ifdef PRINT_DEBUG
431      {
432        // Get the sense of the row
433        int i;
434        printf("rhs sense %c rhs %g\n",si.getRowSense()[rowIndex],
435               si.getRightHandSide()[rowIndex]);
436        const int * indices = si.getMatrixByRow()->getVector(rowIndex).getIndices();
437        const double * elements = si.getMatrixByRow()->getVector(rowIndex).getElements();
438        // for every variable in the constraint
439        for (i=0; i<si.getMatrixByRow()->getVector(rowIndex).getNumElements(); i++){
440          printf("%d (s=%g) %g, ",indices[i],colsol[indices[i]],elements[i]);
441        }
442        printf("\n");
443      }
444#endif
445     
446      //////////////////////////////////////////////////////
447      // Look for a series of                             //
448      // different types of minimal covers.               //
449      // If a minimal cover is found,                     //
450      // lift the associated minimal cover inequality,    //
451      // uncomplement the vars                            //
452      // and add it to the cut set.                       //
453      // After the last type of cover is tried,           //
454      // restore xstar values                             //
455      //////////////////////////////////////////////////////
456     
457      //////////////////////////////////////////////////////
458      // Try to generate a violated                       //
459      // minimal cover greedily from fractional vars      //
460      //////////////////////////////////////////////////////
461     
462      CoinPackedVector cover, remainder; 
463     
464     
465      if (findGreedyCover(rowIndex, krow, b, xstar, cover, remainder) == 1){
466       
467        // Lift cover inequality and add to cut set
468        if (!liftAndUncomplementAndAdd(rowUpper[rowIndex], krow, b,
469                                       complement, rowIndex, cover, 
470                                       remainder, cs)) {
471          // Reset local data and continue to the next iteration
472          // of the rowIndex-loop
473          // I am not sure this is needed but I am just being careful
474          for(k=0; k<krow.getNumElements(); k++) {
475            if (complement[krow.getIndices()[k]]){
476              xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]];
477              complement[krow.getIndices()[k]]=0;       
478            }
479          }
480          krow.setVector(0,NULL,NULL);
481          continue;
482        } 
483      }
484     
485     
486      //////////////////////////////////////////////////////
487      // Try to generate a violated                       //
488      // minimal cover using pseudo John and Ellis logic  //
489      //////////////////////////////////////////////////////
490     
491      // Reset the cover and remainder
492      cover.setVector(0,NULL,NULL);
493      remainder.setVector(0,NULL,NULL);
494     
495      if (findPseudoJohnAndEllisCover(rowIndex, krow, b,
496                                      xstar, cover, remainder) == 1){
497       
498        // (Sequence Independent) Lift cover inequality and add to cut set
499        if (!liftAndUncomplementAndAdd(rowUpper[rowIndex], krow, b,
500                                       complement, rowIndex, cover, 
501                                       remainder, cs)) {
502          // Reset local data and continue to the next iteration
503          // of the rowIndex-loop
504          // I am not sure this is needed but I am just being careful
505          for(k=0; k<krow.getNumElements(); k++) {
506            if (complement[krow.getIndices()[k]]){
507              xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]];
508              complement[krow.getIndices()[k]]=0;       
509            }
510          }
511          krow.setVector(0,NULL,NULL);
512          continue;
513        } 
514       
515        // Skip experiment for now...
516#if 0
517        // experimenting here...
518        // (Sequence Dependent) Lift cover inequality and add to cut set
519        seqLiftAndUncomplementAndAdd(nCols, xstar, complement, rowIndex,
520                                     krow.getNumElements(), b, cover, remainder,
521                                     cs);
522#endif
523      } 
524     
525     
526     
527      //////////////////////////////////////////////////////
528      // Try to generate cuts using covers of unsat       //
529      // vars on reduced krows with John and Ellis logic  //
530      //////////////////////////////////////////////////////
531      CoinPackedVector atOnes;
532      CoinPackedVector fracCover; // different than cover
533     
534      // reset the remainder
535      remainder.setVector(0,NULL,NULL);
536     
537      if (expensiveCuts_||krow.getNumElements()<=15||(!info.inTree&&krow.getNumElements()<=20)) {
538        if (findJohnAndEllisCover(rowIndex, krow, b,
539                                  xstar, fracCover, atOnes, remainder) == 1){
540         
541          // experimenting here...
542          // Sequence Dependent Lifting up on remainders and lifting down on the
543          // atOnes
544          liftUpDownAndUncomplementAndAdd(nCols, xstar, complement, rowIndex,
545                                          krow.getNumElements(), b, fracCover,
546                                          atOnes, remainder, cs);
547        }
548      }
549     
550      //////////////////////////////////////////////////////
551      // Try to generate a violated                       //
552      // minimal cover by considering the                 //
553      // most violated cover problem                      //
554      //////////////////////////////////////////////////////
555     
556     
557      // reset cover and remainder
558      cover.setVector(0,NULL,NULL);
559      remainder.setVector(0,NULL,NULL);
560     
561      // if the size of the krow is "small",
562      //    use an exact algorithm to find the most violated (minimal) cover,
563      // else,
564      //    use an lp-relaxation to find the most violated (minimal) cover.
565      if(krow.getNumElements()<=15||(!info.inTree&&krow.getNumElements()<=20)){
566        if (findExactMostViolatedMinCover(nCols, rowIndex, krow, b,
567                                          xstar, cover, remainder) == 1){
568         
569          // Lift cover inequality and add to cut set
570          if (!liftAndUncomplementAndAdd(rowUpper[rowIndex], krow, b,
571                                         complement, rowIndex, cover, remainder,
572                                         cs)) {
573            // Reset local data and continue to the next iteration
574            // of the rowIndex-loop
575            // I am not sure this is needed but I am just being careful
576            for(k=0; k<krow.getNumElements(); k++) {
577              if (complement[krow.getIndices()[k]]){
578                xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]];
579                complement[krow.getIndices()[k]]=0;       
580              }
581            }
582            krow.setVector(0,NULL,NULL);
583            continue;
584          } 
585        }
586      } 
587      else {
588        if (findLPMostViolatedMinCover(nCols, rowIndex, krow, b,
589                                       xstar, cover, remainder) == 1){
590         
591          // Lift cover inequality and add to cut set
592          if (!liftAndUncomplementAndAdd(rowUpper[rowIndex], krow, b,
593                                         complement, rowIndex, cover, remainder,
594                                         cs)) {
595            // Reset local data and continue to the next iteration
596            // of the rowIndex-loop
597            // I am not sure this is needed but I am just being careful
598            for(k=0; k<krow.getNumElements(); k++) {
599              if (complement[krow.getIndices()[k]]){
600                xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]];
601                complement[krow.getIndices()[k]]=0;       
602              }
603            }
604            krow.setVector(0,NULL,NULL);
605            continue;
606          } 
607        }
608      } 
609     
610     
611     
612      // Reset xstar and complement to their initialized values for the next
613      // go-around
614      int k;
615      if (fabs(b-rowUpper[rowIndex]) > epsilon_) {
616        for(k=0; k<krow.getNumElements(); k++) {
617          if (complement[krow.getIndices()[k]]){
618            xstar[krow.getIndices()[k]]= 1.0-xstar[krow.getIndices()[k]];
619            complement[krow.getIndices()[k]]=0;
620          }
621        }
622      }
623      krow.setVector(0,NULL,NULL);
624#ifdef CGL_DEBUG
625      int nnow = cs.sizeRowCuts();
626      if (nnow>nlast) {
627        const OsiRowCutDebugger * debugger = si.getRowCutDebugger();
628        if (debugger&&debugger->onOptimalPath(si)) {
629          // check cuts okay
630          int k;
631          for (k=nlast;k<nnow;k++) {
632            OsiRowCut rc=cs.rowCut(k);
633            if(debugger->invalidCut(rc)) {
634              printf("itry %d, rhs %g, length %d\n",itry,rhs[itry],length);
635              int i;
636              for (i=0;i<length;i++) {
637                int iColumn = thisColumnIndex[i];
638                printf("column %d, coefficient %g, value %g, bounds %g %g\n",iColumn,
639                       thisElement[i],colsol[iColumn],colLower[iColumn],
640                       colUpper[iColumn]);
641              }
642              if (itry>1) {
643                int length = rowLength[rowIndex];
644                memcpy(thisColumnIndex,column+rowStart[rowIndex],
645                       length*sizeof(int));
646                memcpy(thisElement,elementByRow+rowStart[rowIndex],
647                       length*sizeof(double));
648                printf("Original row had rhs %g and length %d\n",
649                       (itry==2 ? rowLower[rowIndex] :rowUpper[rowIndex]),
650                       length);
651                for (i=0;i<length;i++) {
652                  int iColumn = thisColumnIndex[i];
653                  printf("column %d, coefficient %g, value %g, bounds %g %g\n",iColumn,
654                         thisElement[i],colsol[iColumn],colLower[iColumn],
655                         colUpper[iColumn]);
656                }
657              }
658              assert(!debugger->invalidCut(rc));
659            }
660          }
661        }
662        if (itry>1&&nnow-nlast>kcuts[itry-2]) {
663          printf("itry %d gave %d cuts as against %d for itry %d\n",
664                 itry,nnow-nlast,kcuts[itry-2],itry-2);
665        }
666        kcuts[itry]=nnow-nlast;
667        nlast=nnow;
668      }
669#endif
670    }
671  }
672  // Clean up: free allocated memory
673  if (toCheck != rowsToCheck_)
674     delete[] toCheck;
675  delete[] xstar;
676  delete[] complement;
677  delete [] thisColumnIndex;
678  delete [] thisElement;
679  delete [] back;
680  delete [] vub;
681  delete [] vubRow;
682  delete [] vubValue;
683  delete [] vlbRow;
684  delete [] vlbValue;
685  delete [] effectiveLower;
686  delete [] effectiveUpper;
687}
688
689void
690CglKnapsackCover::setTestedRowIndices(int num, const int* ind)
691{
692   if (rowsToCheck_)
693      delete[] rowsToCheck_;
694   numRowsToCheck_ = num;
695   if (num > 0) {
696      rowsToCheck_ = new int[num];
697      CoinCopyN(ind, num, rowsToCheck_);
698   }
699}
700
701//-------------------------------------------------------------
702// Lift and uncomplement cut. Add cut to the cutset
703//-------------------------------------------------------------------
704int 
705CglKnapsackCover::liftAndUncomplementAndAdd(
706         double rowub,
707         CoinPackedVector & krow,
708         double & b,
709         int * complement,
710         int row,
711         CoinPackedVector & cover,
712         CoinPackedVector & remainder,
713         OsiCuts & cs ) const
714{
715  CoinPackedVector cut;
716  double cutRhs = cover.getNumElements() - 1.0;
717  int goodCut=1;
718 
719  if (remainder.getNumElements() > 0){
720    // Construct lifted cover cut
721    if (!liftCoverCut( 
722                      b, krow.getNumElements(), 
723                      cover, remainder,
724                      cut )) 
725      goodCut= 0; // no cut
726  }
727  // The cover consists of every variable in the knapsack.
728  // There is nothing to lift, so just add cut           
729  else {
730    cut.reserve(cover.getNumElements());
731    cut.setConstant(cover.getNumElements(),cover.getIndices(),1.0);
732  }
733 
734  // Uncomplement the complemented variables in the cut
735  int k;
736  if (fabs(b-rowub)> epsilon_) {
737    double * elements = cut.getElements();
738    int * indices = cut.getIndices();
739    for (k=0; k<cut.getNumElements(); k++){
740      if (complement[indices[k]]) {
741        // Negate the k'th element in packedVector cut
742        // and correspondingly adjust the rhs
743        elements[k] *= -1;
744        cutRhs += elements[k];
745      }
746    }
747  }
748  if (goodCut) {
749    // Create row cut. Effectiveness defaults to 0.
750    OsiRowCut rc;
751    rc.setRow(cut);
752#ifdef CGL_DEBUG
753    {
754      double * elements = cut.getElements();
755      int * indices = cut.getIndices();
756      int n=cut.getNumElements();
757      for (k=0; k<n; k++){
758        assert(indices[k]>=0);
759        assert(elements[k]);
760        assert (fabs(elements[k])>1.0e-12);
761      }
762    }
763#endif
764    rc.setLb(-DBL_MAX);
765    rc.setUb(cutRhs);
766    //  rc.setEffectiveness(0);
767    // Todo: put in a more useful measure such as  the violation.
768   
769    // Add row cut to the cut set 
770#ifdef PRINT_DEBUG
771    {
772      int k;
773      printf("cutrhs %g %d elements\n",cutRhs,cut.getNumElements());
774      double * elements = cut.getElements();
775      int * indices = cut.getIndices();
776      for (k=0; k<cut.getNumElements(); k++){
777        printf("%d %g\n",indices[k],elements[k]);
778      }
779    }
780#endif
781    cs.insert(rc);
782   
783    return 1;
784  } else {
785    return 0;
786  }
787}
788
789//-------------------------------------------------------------------
790// deriveAKnapsack - returns 1 if the method is able to
791//                  derive a canonical knapsack inequality
792//                  in binary variables of the form ax<=b
793//                  from the rowIndex-th row of the constraint matrix.
794//                  returns 0, otherwise.
795//  Precondition: complement must be 0'd out!!!
796//-------------------------------------------------------------------
797int 
798CglKnapsackCover::deriveAKnapsack(
799       const OsiSolverInterface & si, 
800       OsiCuts & cs,
801       CoinPackedVector & krow, 
802       bool treatAsLRow,
803       double & b,
804       int *  complement,
805       double *  xstar,
806       int rowIndex,
807       int numberElements,
808       const int * index,
809       const double * element) const
810{
811
812  // Fix to https://projects.coin-or.org/Cbc/ticket/30
813  {
814#ifndef COIN_DEVELOP
815    if (numberElements==0) return 0;
816#else
817    assert (numberElements>0);
818#endif
819  }
820
821  int i;
822
823  krow.clear();
824
825  // if the matrixRow represent a ge inequality, then
826  //     leMatrixRow == -matrixRow  // otherwise
827  //     leMatrixRow == matrixRow.
828
829  CoinPackedVector leMatrixRow(numberElements,index,element); 
830
831  double maxKrowElement = -DBL_MAX;
832  double minKrowElement = DBL_MAX;
833 
834
835  if (treatAsLRow) {
836    // treat as if L row
837  } else {
838    // treat as if G row
839    b=-b;
840    std::transform(leMatrixRow.getElements(),
841                   leMatrixRow.getElements() + leMatrixRow.getNumElements(),
842                   leMatrixRow.getElements(),
843                   std::negate<double>());
844  }
845 
846  // nBinUnsat is a counter for the number of unsatisfied
847  // (i.e. fractional) binary vars 
848  int nBinUnsat =0;
849  const double * colupper = si.getColUpper();
850  const double * collower = si.getColLower();
851 
852  // At this point, leMatrixRow and b represent a le inequality in general
853  // variables.
854  // To derive a canonical knapsack inequality in free binary variable,
855  // process out the continuous & non-binary integer & fixed binary variables.
856  // If the non-free-binary variables can be appropriately bounded,
857  // net them out of the constraint, otherwise abandon this row and return 0.
858
859  const int * indices = leMatrixRow.getIndices();
860  const double * elements = leMatrixRow.getElements();
861  // for every variable in the constraint
862  for (i=0; i<leMatrixRow.getNumElements(); i++){
863    // if the variable is not a free binary var
864    if ( !si.isFreeBinary(indices[i]) ) {
865      // and the coefficient is strictly negative
866      if(elements[i]<-epsilon_){
867        // and the variable has a finite upper bound
868        if (colupper[indices[i]] < si.getInfinity()){
869          // then replace the variable with its upper bound.
870          b=b-elements[i]*colupper[indices[i]];
871        } 
872        else {
873          return 0;
874        }
875      }
876      // if the coefficient is strictly positive
877      else if(elements[i]>epsilon_){
878        // and the variable has a finite lower bound
879        if (collower[indices[i]] > -si.getInfinity()){
880          // then replace the variable with its lower bound.
881          b=b-elements[i]*collower[indices[i]];
882        }
883        else {
884          return 0;
885        }
886      }
887      // note: if the coefficient is zero, the variable is not included in the
888      //       knapsack inequality.
889    }
890    // else the variable is a free binary var and is included in the knapsack
891    // inequality.
892    // note: the variable is included regardless of its solution value to the
893    // lp relaxation.
894    else{
895      krow.insert(indices[i], elements[i]);
896
897      // if the binary variable is unsatified (i.e. has fractional value),
898      // increment the counter.
899      if(xstar[indices[i]] > epsilon_ && xstar[indices[i]] < onetol_)
900        nBinUnsat++;
901
902      // keep track of the largest and smallest elements in the knapsack
903      // (the idea is if there is not a lot of variation in the knapsack
904      // coefficients, it is unlikely we will find a violated minimal
905      // cover from this knapsack so don't even bother trying)
906      if (fabs(elements[i]) > maxKrowElement) 
907        maxKrowElement = fabs(elements[i]);
908      if (fabs(elements[i]) < minKrowElement) 
909        minKrowElement = fabs(elements[i]);
910    }
911  }
912 
913  // If there's little variation in the knapsack coefficients, return 0.
914  // If there are no unsatisfied binary variables, return.
915  // If there's only one binary, return. 
916  // ToDo: but why return if 2 binary? ...there was some assumption in the
917  // findVioMinCover..(?)   
918  // Anyway probing will probably find something
919  if (krow.getNumElements() < 3 ||
920      nBinUnsat == 0 ||
921      maxKrowElement-minKrowElement < 1.0e-3*maxKrowElement ) {
922    return 0;
923  }
924
925  // However if we do decide to do when count is two - look carefully
926  if (krow.getNumElements()==2) {
927    const int * indices = krow.getIndices();
928    double * elements = krow.getElements();
929    double sum=0.0;
930    for(i=0; i<2; i++){
931      int iColumn = indices[i];
932      sum += elements[i]*xstar[iColumn];
933    }
934    if (sum<b-1.0e-4) {
935      return 0;
936    } else {
937
938#ifdef PRINT_DEBUG
939      printf("*** Doubleton Row is ");
940      for(i=0; i<2; i++){
941        int iColumn = indices[i];
942        sum += elements[i]*xstar[iColumn];
943        printf("%d (coeff = %g, value = %g} ",indices[i],
944               elements[i],xstar[iColumn]);
945      }
946      printf("<= %g - go for it\n",b);
947#endif
948
949    }
950  }
951
952
953  // At this point krow and b represent a le inequality in binary variables. 
954  // To obtain an le inequality with all positive coefficients, complement
955  // any variable with a negative coefficient by changing the sign of
956  // the coefficient, adjusting the rhs, and adjusting xstar, the column
957  // solution vector.
958  {
959     const int s = krow.getNumElements();
960     const int * indices = krow.getIndices();
961     double * elements = krow.getElements();
962     for(i=0; i<s; i++){
963         if (elements[i] < -epsilon_) {
964           complement[indices[i]]= 1;
965           elements[i] *= -1;
966           b+=elements[i]; 
967           xstar[indices[i]]=1.0-xstar[indices[i]];
968        }
969     }
970  }
971
972  // Quick feasibility check.
973  // If the problem is infeasible, add an infeasible col cut to cut set
974  // e.g. one that has lb > ub.
975  // TODO: test this scenario in BCP
976  if (b < 0 ){ 
977    OsiColCut cc;
978    int index = krow.getIndices()[0]; 
979    const double fakeLb = colupper[index] + 1.0;;  // yes, colupper.
980    const double fakeUb = collower[index];
981    assert( fakeUb < fakeLb );
982    cc.setLbs( 1, &index, &fakeLb);
983    cc.setUbs( 1, &index, &fakeLb);
984    cc.setEffectiveness(DBL_MAX);
985    cs.insert(cc);
986#ifdef PRINT_DEBUG
987    printf("Cgl: Problem is infeasible\n");
988#endif
989  }
990 
991  // At this point, krow and b represent a le inequality with postive
992  // coefficients.
993  // If any coefficient a_j > b, then x_j = 0, return 0
994  // If any complemented var has coef a_j > b, then x_j = 1, return 0
995  int fixed = 0;
996  CoinPackedVector fixedBnd; 
997  for(i=0; i<krow.getNumElements(); i++){
998    if (krow.getElements()[i]> b){
999      fixedBnd.insert(krow.getIndices()[i],complement[krow.getIndices()[i]]);
1000#ifdef PRINT_DEBUG   
1001      printf("Variable %i being fixed to %i due to row %i.\n",
1002             krow.getIndices()[i],complement[krow.getIndices()[i]],rowIndex); 
1003#endif
1004      fixed = 1;     
1005    }
1006  }
1007
1008  // After all possible variables are fixed by adding a column cut with
1009  // equivalent lower and upper bounds, return
1010  if (fixed) {
1011    OsiColCut cc;
1012    cc.setLbs(fixedBnd);
1013    cc.setUbs(fixedBnd);
1014    cc.setEffectiveness(DBL_MAX);
1015    return 0; 
1016  } 
1017
1018  return 1;
1019}
1020
1021
1022//-------------------------------------------------------------------
1023// deriveAKnapsack - returns 1 if the method is able to
1024//                  derive a cannonical knapsack inequality
1025//                  in binary variables of the form ax<=b
1026//                  from the rowIndex-th row of the constraint matrix.
1027//                  returns 0, otherwise.
1028//  Precondition: complement must be 0'd out!!!
1029//-------------------------------------------------------------------
1030int 
1031CglKnapsackCover::deriveAKnapsack(
1032       const OsiSolverInterface & si, 
1033       OsiCuts & cs,
1034       CoinPackedVector & krow, 
1035       double & b,
1036       int *  complement,
1037       double *  xstar,
1038       int rowIndex,
1039       const CoinPackedVectorBase & matrixRow ) const
1040{
1041  // Get the sense of the row
1042  const char  rowsense = si.getRowSense()[rowIndex];
1043 
1044  // Skip equality and unbounded rows
1045  if  (rowsense=='E' || rowsense=='N') {
1046    return 0; 
1047  }
1048 
1049  bool treatAsLRow =  (rowsense=='L');
1050  const int * indices = matrixRow.getIndices();
1051  const double * elements = matrixRow.getElements();
1052  int numberElements = matrixRow.getNumElements();
1053  return deriveAKnapsack( si, cs, krow, treatAsLRow, b, complement,
1054                          xstar, rowIndex, numberElements, indices,
1055                          elements);
1056}
1057
1058//--------------------------------------------------
1059// Find a violated minimal cover from
1060// a canonical form knapsack inequality by
1061// solving the lp relaxation of the
1062// -most- violated cover problem.
1063// Postprocess to ensure minimality.
1064// -----------------------------------------
1065int 
1066CglKnapsackCover::findLPMostViolatedMinCover(
1067      int nCols,
1068      int row,
1069      CoinPackedVector & krow,
1070      double & b,
1071      double * xstar, 
1072      CoinPackedVector & cover,
1073      CoinPackedVector & remainder) const
1074{
1075 
1076  // Assumes krow and b describe a knapsack inequality in canonical form
1077
1078  // Given a knapsack inequality sum a_jx_j <= b, and optimal lp solution
1079  // xstart, a violated minimal cover inequality exists if the following 0-1
1080  // programming problem has an optimal objective function value (oofv) < 1
1081  //     oofv = min sum (1-xstar_j)z_j
1082  //            s.t. sum a_jz_j > b
1083  //            z binary
1084 
1085  // The vector z is an incidence vector, defining the cover R with the
1086  // associated cover inequality:
1087  //    (sum j in R) x_j <= |R|-1
1088 
1089  // This problem is itself a (min version of the) knapsack problem
1090  // but with a unsightly strict inequalty.
1091 
1092  // To transform transform it into a max version,
1093  // complement the z's, z_j=1-y_j.
1094  // To compensate for the strict inequality, subtract epsilon from the rhs.
1095 
1096  //     oofv = (sum (1-xstar_j))-  max sum (1-xstar)y_j
1097  //                        s.t. sum a_jy_j <= (sum over j)a_j - b (- EPSILON)
1098  //                             y binary
1099 
1100  // If oofv < 1, then a violated min cover inequality has
1101  // incidence vector z with elements z_j=1-y_j and rhs= num of nonzero's in
1102  // z, i.e. the number of 0's in y.
1103
1104  // If the size of the knapsack is "small", we solve the problem exactly.
1105  // If the size of the knapsack is large, we solve the (simpler) lp relaxation
1106  // of the  knapsack problem and postprocess to ensure the construction of a
1107  // minimimal cover.
1108 
1109  // We also assume that testing/probing/fixing based on the knapsack structure
1110  // is done elsewhere. Only convenient-to-do sanity checks are done here.
1111  // (We do not assume that data is integer.)
1112
1113  double elementSum = krow.sum();
1114
1115  // Redundant/useless adjusted rows should have been trapped in the
1116  // transformation to the canonical form of knapsack inequality
1117  if (elementSum < b + epsilon_) {
1118    return -1; 
1119  }
1120 
1121  // Order krow in nonincreasing order of coefObj_j/a_j. 
1122  // (1-xstar_1)/a_1 >= (1-xstar_2)/a_2 >= ... >= (1-xstar_n)/a_n   
1123  // by defining this full-storage array "ratio" to be the external sort key.
1124  double * ratio= new double[nCols];
1125  memset(ratio, 0, (nCols*sizeof(double)));
1126
1127  int i;
1128  for (i=0; i<krow.getNumElements(); i++){
1129    if (fabs(krow.getElements()[i])> epsilon_ ){
1130      ratio[krow.getIndices()[i]]=
1131         (1.0-xstar[krow.getIndices()[i]]) / (krow.getElements()[i]);
1132    }
1133    else {
1134      ratio[krow.getIndices()[i]] = 0.0;
1135    }
1136  }
1137
1138  // ToDo: would be nice to have sortkey NOT be full-storage vector
1139  CoinDecrSolutionOrdered dso(ratio);
1140  krow.sort(dso);   
1141
1142  // Find the "critical" element index "r" in the knapsack lp solution
1143  int r = 0;
1144  double sum = krow.getElements()[0];
1145  while ( sum <= (elementSum - b - epsilon_ ) ){
1146    r++;
1147    sum += krow.getElements()[r];
1148  }   
1149 
1150  // Note: It is possible the r=0, and you get a violated minimal cover
1151  // if (r=0), then you've got a var with a really large coeff. compared
1152  //   to the rest of the row.
1153  // r=0 says trivially that the
1154  //   sum of ALL the binary vars in the row <= (cardinality of all the set -1)
1155  // Note: The cover may not be minimal if there are alternate optimals to the
1156  // maximization problem, so the cover must be post-processed to ensure
1157  // minimality.
1158 
1159  // "r" is the critical element
1160  // The lp relaxation column solution is:
1161  // y_j = 1 for  j=0,...,(r-1)
1162  // y_r = (elementSum - b - sum + krow.element()[r])/krow.element()[r]
1163  // y_j = 0 for j=r+1,...,krow.getNumElements()
1164 
1165  // The number of nonzeros in the lp column solution is r+1
1166 
1167  // if oofv to the lp knap >= 1, then no violated min cover is possible
1168  int nCover;
1169
1170  double lpoofv=0.0;
1171  for (i=r+1; i<krow.getNumElements(); i++){
1172    lpoofv += (1-xstar[krow.getIndices()[i]]);
1173  }
1174  double ipofv = lpoofv + (1-xstar[krow.getIndices()[r]]);
1175
1176  // Couldn't find an lp violated min cover inequality
1177  if (ipofv > 1.0 - epsilon_){   
1178    delete [] ratio;
1179    return -1;
1180  }
1181  else {
1182    // Partition knapsack into cover and noncover (i.e. remainder)
1183    // pieces
1184    nCover = krow.getNumElements() - r;
1185    double coverSum =0.0;
1186    cover.reserve(nCover);
1187    remainder.reserve(r);
1188   
1189    for (i=r; i<krow.getNumElements(); i++){
1190      cover.insert(krow.getIndices()[i],krow.getElements()[i]);
1191      coverSum += krow.getElements()[i];
1192    }
1193    for (i=0; i<r; i++){
1194      remainder.insert(krow.getIndices()[i],krow.getElements()[i]);
1195    }
1196   
1197    if (coverSum <= b){
1198#ifdef PRINT_DEBUG
1199      printf("The identified cover is NOT a cover\n");
1200      abort();
1201#endif
1202      delete [] ratio;
1203      return -1;
1204    }
1205   
1206    // Sort cover in terms of knapsack row coefficients   
1207    cover.sortDecrElement();
1208   
1209   
1210    // We have a violated cover inequality.
1211    // Construct a -minimal- violated cover
1212    // by testing and potentially tossing smallest
1213    // elements
1214    double oneLessCoverSum = coverSum - cover.getElements()[nCover-1];
1215    while (oneLessCoverSum > b+1.0e-12){
1216      // move the excess cover member into the set of remainders
1217      remainder.insert(cover.getIndices()[nCover-1],
1218                       cover.getElements()[nCover-1]);
1219      cover.truncate(nCover-1);
1220      nCover--;
1221      oneLessCoverSum -= cover.getElements()[nCover-1];
1222    }
1223
1224    if (nCover<2){
1225#ifdef PRINT_DEBUG
1226      printf("nCover < 2...aborting\n");
1227      abort();
1228#endif
1229      delete [] ratio;
1230      return -1;
1231    }
1232   
1233#ifdef PRINT_DEBUG   /* debug */
1234    printf("\
1235Lp relax of most violated minimal cover: row %i has cover of size %i.\n",
1236           row,nCover);
1237    //double sumCover = 0.0;
1238    for (i=0; i<cover.getNumElements(); i++){
1239      printf("index %i, element %g, xstar value % g \n",
1240             cover.getIndices()[i],cover.getElements()[i],
1241             xstar[cover.getIndices()[i]]);
1242      //sumCover += cover.getElements()[i];
1243    }
1244    printf("The b = %g, and the cover sum is %g\n\n", b, cover.sum());
1245#endif
1246
1247#ifdef P0201
1248      double ppsum=0.0;
1249      for (i=0; i<nCover; i++){
1250        ppsum += p0201[krow.getIndices()[i]];
1251      }
1252       
1253      if (ppsum > nCover-1){
1254          printf("\
1255\nBad cover from lp relax of most violated cover..aborting\n");
1256          abort();
1257      }
1258#endif
1259   
1260    /* clean up */
1261    delete [] ratio;
1262    return  1;
1263  }
1264}
1265
1266
1267//--------------------------------------------------
1268// Find a violated minimal cover from
1269// a canonical form knapsack inequality by
1270// solving the -most- violated cover problem
1271// and postprocess to ensure minimality
1272// -----------------------------------------
1273int 
1274CglKnapsackCover::findExactMostViolatedMinCover(
1275        int nCols,
1276        int row,
1277        CoinPackedVector & krow,
1278        double b, 
1279        double *  xstar, 
1280        CoinPackedVector & cover,
1281        CoinPackedVector & remainder) const 
1282{
1283 
1284  // assumes the row is in canonical knapsack form
1285 
1286  // A violated min.cover inequality exists if the
1287  // opt obj func value (oofv) < 1:
1288  //     oofv = min sum (1-xstar_j)z_j
1289  //            s.t. sum a_jz_j > b
1290  //            x binary
1291
1292  //     The vector z is the incidence vector
1293  //     defines the set R and the cover inequality.
1294  //      (sum j in R) x_j <= |R|-1
1295
1296  //    This is the min version of the knapsack problem.
1297  //    (note that strict inequalty...bleck)
1298
1299  //    To obtain the max version, complement the z's, z_j=1-y_j and
1300  //    adjust the constraint.
1301
1302  //    oofv = (sum (1-xstar_j))-  max sum (1-xstar)y_j
1303  //                       s.t. sum a_jy_j <= (sum over j)a_j - b (- EPSILON)]
1304  //                   y binary
1305
1306  // If oofv < 1, violated min cover inequality has
1307  //    incidence vector z=1-y and rhs= num of nonzero's in z, i.e.
1308  //    the number 0 in y.
1309
1310  //    We solve the  0-1 knapsack problem by explicit ennumeration
1311
1312  double elementSum = krow.sum();
1313
1314  // Redundant/useless adjusted rows should have been trapped in
1315  // transformation to canonical form of knapsack inequality
1316  if (elementSum < b + epsilon_) {
1317#ifdef PRINT_DEBUG
1318    printf("Redundant/useless adjusted row\n");
1319#endif
1320    return -1; 
1321  }
1322
1323  // Order krow in nonincreasing order of coefObj_j/a_j. 
1324  // (1-xstar_1)/a_1 >= (1-xstar_2)/a_2 >= ... >= (1-xstar_n)/a_n   
1325  // by defining this full-storage array "ratio" to be the external sort key. 
1326  double * ratio= new double[nCols];
1327  memset(ratio, 0, (nCols*sizeof(double)));
1328
1329  int i;
1330  {
1331     const int * indices = krow.getIndices();
1332     const double * elements = krow.getElements();
1333     for (i=0; i<krow.getNumElements(); i++){
1334        if (fabs(elements[i])> epsilon_ ){
1335           ratio[indices[i]]= (1.0-xstar[indices[i]]) / elements[i];
1336        }
1337        else {
1338           ratio[indices[i]] = 0.0;
1339        }
1340     }
1341  }
1342
1343  // ToDo: would be nice to have sortkey NOT be full-storage vector
1344  CoinDecrSolutionOrdered dso(ratio);
1345  krow.sort(dso);   
1346
1347#ifdef CGL_DEBUG
1348  // sanity check
1349  for ( i=1; i<krow.getNumElements(); i++ ) {
1350    double ratioim1 =  ratio[krow.getIndices()[i-1]];
1351    double ratioi= ratio[krow.getIndices()[i]];
1352    assert( ratioim1 >= ratioi );
1353  } 
1354#endif 
1355 
1356  // Recall:
1357  // oofv = (sum (1-xstar_j))-  max sum (1-xstar)y_j
1358  //           s.t. sum a_jy_j <= (sum over j)a_j - b (- epsilon_)]
1359  //           y binary
1360
1361  double objConst = 0.0;
1362  double exactOptVal = -1.0;
1363  int * exactOptSol = new int[krow.getNumElements()];
1364  double * p = new double[krow.getNumElements()];
1365  double * w = new double[krow.getNumElements()];
1366  int kk;
1367  for (kk=0; kk<krow.getNumElements(); kk++){
1368    p[kk]=1.0-xstar[krow.getIndices()[kk]];
1369    w[kk]=krow.getElements()[kk];
1370    objConst+=p[kk];
1371  }
1372 
1373  // vectors are indexed in ratioSortIndex order
1374  exactSolveKnapsack(krow.getNumElements(), (elementSum-b-epsilon_), p, w,
1375                     exactOptVal, exactOptSol);
1376 
1377  if(objConst-exactOptVal < 1){
1378    cover.reserve(krow.getNumElements());
1379    remainder.reserve(krow.getNumElements());
1380
1381    // Partition the krow into the cover and the remainder.
1382    // The cover is complement of solution.
1383    double coverElementSum = 0;
1384    for(kk=0;kk<krow.getNumElements();kk++){
1385      if(exactOptSol[kk]==0){
1386        cover.insert(krow.getIndices()[kk],krow.getElements()[kk]);
1387        coverElementSum += krow.getElements()[kk];
1388      }
1389      else {
1390        remainder.insert(krow.getIndices()[kk],krow.getElements()[kk]);
1391      }
1392    }
1393
1394    cover.sortDecrElement();
1395
1396    // We have a violated cover inequality.
1397    // Construct a -minimal- violated cover
1398    // by testing and potentially tossing smallest
1399    // elements
1400    double oneLessCoverElementSum =
1401       coverElementSum - cover.getElements()[cover.getNumElements()-1];
1402    while (oneLessCoverElementSum > b){
1403      // move the excess cover member into the set of remainders
1404      remainder.insert(cover.getIndices()[cover.getNumElements()-1],
1405                       cover.getElements()[cover.getNumElements()-1]);
1406      cover.truncate(cover.getNumElements()-1);
1407      oneLessCoverElementSum -= cover.getElements()[cover.getNumElements()-1];
1408    }
1409
1410#ifdef PRINT_DEBUG
1411    printf("Exact Most Violated Cover: row %i has cover of size %i.\n",
1412           row,cover.getNumElements());
1413    //double sumCover = 0.0;
1414    for (i=0; i<cover.getNumElements(); i++){
1415      printf("index %i, element %g, xstar value % g \n",
1416             cover.getIndices()[i], cover.getElements()[i],
1417             xstar[cover.getIndices()[i]]);
1418      //sumCover += cover.getElements()[i];
1419    }
1420    printf("The b = %g, and the cover sum is %g\n\n", b, cover.sum() );
1421#endif
1422   
1423    // local clean up
1424    delete [] exactOptSol;
1425    delete [] p;
1426    delete [] w;
1427    delete [] ratio;
1428   
1429    return  1; // found an exact one
1430  }
1431
1432  // local clean up
1433  delete [] exactOptSol;
1434  delete [] p;
1435  delete [] w;
1436  delete [] ratio;
1437 
1438  return 0; // didn't find an exact one
1439
1440} 
1441
1442//-------------------------------------------------------------------
1443// Find Pseudo John-and-Ellis Cover
1444//
1445// only generates -violated- minimal covers
1446//-------------------------------------------------------------------
1447int
1448CglKnapsackCover::findPseudoJohnAndEllisCover(
1449     int row,
1450     CoinPackedVector & krow,
1451     double & b,
1452     double * xstar, 
1453     CoinPackedVector & cover, 
1454     CoinPackedVector & remainder) const
1455
1456{
1457  // semi-mimic of John&Ellis's approach without taking advantage of SOS info
1458  // RLH: They find a minimal cover on unsatisfied variables, but it is
1459  // not guaranteed to be violated by currently solution
1460
1461  // going for functional now, will make efficient when working
1462 
1463  // look at unstatisfied binary vars with nonzero row coefficients only
1464  // get row in canonical form (here row is in canonical form)
1465  // if complement var, complement soln val too. (already done)
1466  // (*) sort in increasing value of soln val
1467  // track who is the biggest coef and it's index.
1468  // if biggest > adjRhs, skip row. Bad knapsack.
1469  // margin = adjRhs
1470  // idea: if (possibly compl) soln >= .5 round up, else round down
1471  // they do more, but that's the essence
1472  // go through the elements {
1473  // if round down, skip
1474  // if round up, add to element to cover. adjust margin
1475  // if current element = biggest, then get next biggest
1476  // if biggest > marg, you've got a cover. stop looking               
1477  // else try next element in the loop
1478  // }       
1479  // (*)RLH: I'm going to sort in decreasing order of soln val
1480  // b/c var's whose soln < .5 in can. form get rounded down
1481  // and skipped. If you can get a min cover of the vars
1482  // whose soln is >= .5, I believe this gives the same as J&E.
1483  // But if not, maybe I can get something more.
1484       
1485  // (**)By checking largest value left, they ensure a minimal cover
1486  // on the unsatisfied variables
1487     
1488  // if you have a cover
1489  // sort the elements to be lifted in order of their reduced costs.
1490  // lift in this order.
1491  // ...I don't understand their lifting, so for now use sequence-indep lifting
1492
1493  // J&E employ lifting up and down.
1494  // Here I'm including the vars at one in the cover.
1495  // Adding these vars back in may cause the minimality of the cover to lost.
1496  // So, post-processing to establish minimality is required.
1497 
1498  cover.reserve(krow.getNumElements());
1499  remainder.reserve(krow.getNumElements());
1500
1501  double unsatRhs = b;
1502
1503  // working info on unsatisfied vars
1504  CoinPackedVector unsat;
1505  unsat.reserve(krow.getNumElements());
1506
1507  // working info on vars with value one
1508  CoinPackedVector atOne;
1509  atOne.reserve(krow.getNumElements());
1510
1511  // partition the (binary) variables in the canonical knapsack
1512  // into those at zero, those at fractions, and those at one.
1513  // Note: no consideration given to whether variables are free
1514  // or fixed at binary values.
1515  // Note: continuous and integer vars have already been netted out
1516  // to derive the canonical knapsack form
1517  int i;
1518  for (i=0; i<krow.getNumElements(); i++){
1519
1520    if (xstar[krow.getIndices()[i]] > onetol_){
1521      atOne.insert(krow.getIndices()[i],krow.getElements()[i]); 
1522      unsatRhs -= krow.getElements()[i];
1523    }   
1524    else if (xstar[krow.getIndices()[i]] >= epsilon_){
1525      unsat.insert(krow.getIndices()[i],krow.getElements()[i]) ;
1526    }
1527    else { 
1528      remainder.insert(krow.getIndices()[i],krow.getElements()[i]);
1529    }
1530  }
1531
1532  // sort the indices of the unsat var in order of decreasing solution value
1533  CoinDecrSolutionOrdered decrSol(xstar);
1534  unsat.sort(decrSol);
1535 
1536#ifdef CGL_DEBUG
1537  // sanity check
1538  for (i=1; i<unsat.getNumElements(); i++){
1539    double xstarim1= xstar[unsat.getIndices()[i-1]];
1540    double xstari= xstar[unsat.getIndices()[i]];
1541    assert( xstarim1 >= xstari );
1542  }
1543#endif
1544
1545  // get the largest coefficient among the unsatisfied variables   
1546  double bigCoef= 0.0;
1547  // double temp;
1548  int bigIndex = 0;
1549  for (i=0; i<unsat.getNumElements(); i++){
1550    if (unsat.getElements()[i]>bigCoef){
1551      bigCoef = unsat.getElements()[i];
1552      bigIndex = i;
1553    }
1554  }
1555 
1556  // initialize
1557  i=0;
1558  double margin = unsatRhs;
1559  int gotCover=0;
1560  int j;
1561 
1562  // look in order through the unsatisfied vars which along with the
1563  // the max element defines a cover
1564  while (i<unsat.getNumElements() && !gotCover){
1565    margin -= unsat.getElements()[i];
1566   
1567    // get the biggest row coef downstream in the given order   
1568    if (i == bigIndex){
1569      bigCoef = 0.0;
1570      bigIndex = 0;
1571      for (j=i+1; j<unsat.getNumElements(); j++){
1572        double temp = unsat.getElements()[j];
1573        if (temp > bigCoef ){
1574          bigCoef = temp;
1575          bigIndex = j;
1576        }
1577      }
1578    }
1579   
1580    if (bigCoef > margin+epsilon2_) gotCover = 1;
1581    i++;         
1582  }
1583
1584  // J&E approach; get first single one element that fills the margin   
1585  if(gotCover){
1586    j=i;
1587    if (j<unsat.getNumElements()){ // need this "if" incase nUnsat=1   
1588      while (unsat.getElements()[j]< margin) {
1589        j++;
1590      }
1591      // switch members so that first nUnsat define the cover
1592      unsat.swap(i,j);
1593      i++;
1594    }
1595   
1596    // check that detected cover is violated
1597    // (would we want to save incase it's violated later?)
1598    int nCover = i;
1599    double coverElementSum = 0.0;
1600    double coverXstarSum = 0.0;
1601    int k;
1602    for (k=0; k<nCover; k++){
1603      coverElementSum += unsat.getElements()[k];
1604      coverXstarSum +=  xstar[unsat.getIndices()[k]];
1605    }
1606   
1607    // Split the unsatisfied elements into those in the cover and those
1608    // not in the cover. The elements not in the cover are considered
1609    // remainders. Variables atOne belong to the cover
1610
1611    // Test if the detected cover is violated
1612    if (coverXstarSum > (nCover-1) && coverElementSum > unsatRhs+epsilon2_){
1613      for (i=nCover; i<unsat.getNumElements(); i++) {
1614        remainder.insert(unsat.getIndices()[i],unsat.getElements()[i]);
1615      }
1616      unsat.truncate(nCover);
1617      cover = unsat;
1618      cover.append(atOne);
1619
1620      for (k=nCover; k<cover.getNumElements(); k++){
1621        coverElementSum+=cover.getElements()[k];
1622        coverXstarSum+=xstar[cover.getIndices()[k]];
1623      }
1624     
1625      // Sanity check
1626      int size = cover.getNumElements() + remainder.getNumElements(); 
1627      int krowsize = krow.getNumElements();
1628      assert( size == krowsize );
1629     
1630      // Sort cover in terms of knapsack row coefficients   
1631      cover.sortDecrElement();
1632
1633    // New!
1634    // We have a violated cover inequality.
1635    // Construct a -minimal- violated cover
1636    // by testing and potentially tossing smallest
1637    // elements
1638    double oneLessCoverElementSum =
1639       coverElementSum - cover.getElements()[cover.getNumElements()-1];
1640    while (oneLessCoverElementSum > b){
1641      // move the excess cover member into the set of remainders
1642      remainder.insert(cover.getIndices()[cover.getNumElements()-1],
1643                       cover.getElements()[cover.getNumElements()-1]);
1644      cover.truncate(cover.getNumElements()-1);
1645      oneLessCoverElementSum -= cover.getElements()[cover.getNumElements()-1];
1646    }
1647 
1648#ifdef PRINT_DEBUG
1649    if (coverXstarSum > (nCover-1) && coverElementSum > b){
1650      printf("John and Ellis: row %i has cover of size %i.\n",
1651             row,cover.getNumElements());
1652      //double sumCover = 0.0;
1653      for (i=0; i<cover.getNumElements(); i++){
1654        printf("index %i, element %g, xstar value % g \n",
1655               cover.getIndices()[i], cover.getElements()[i],
1656               xstar[cover.getIndices()[i]]);
1657      }
1658      printf("The b = %g, and the cover element sum is %g\n\n",b,cover.sum());
1659    }
1660#endif
1661
1662    }
1663    // if minimal cover is not violated, turn gotCover off
1664    else {
1665//    printf("heuristically found minimal cover is NOT violated by current lp solution");
1666      gotCover = 0;
1667    }         
1668  }
1669
1670  // If no minimal cover was found, pack it in   
1671  if (!gotCover || cover.getNumElements() < 2) {
1672    return -1;
1673  }
1674 
1675  return  1;
1676}
1677
1678
1679
1680
1681
1682//-------------------------------------------------------------------
1683// Find a "approx" John-and-Ellis Cover
1684// (i.e. that this approximates John & Ellis code)
1685//  Test to see if we generate the same covers, and lifted cuts
1686//
1687// generates minimal covers, not necessarily violated ones.
1688//-------------------------------------------------------------------
1689int
1690CglKnapsackCover::findJohnAndEllisCover(
1691     int row,
1692     CoinPackedVector & krow,
1693     double & b,
1694     double * xstar, 
1695     CoinPackedVector & fracCover, 
1696     CoinPackedVector & atOne,
1697     CoinPackedVector & remainder) const
1698
1699{
1700  // John Forrest and Ellis Johnson's approach as I see it.
1701  // RLH: They find a minimal cover on unsatisfied variables,
1702  // which may not be violated by the solution to the lp relaxation
1703
1704  // "functional before efficient" is my creed.
1705 
1706  // look at unstatisfied binary vars with nonzero row coefficients only
1707  // get row in canonical form (here krow is in canonical form)
1708  // if complement var, complement soln val too. (already done)
1709  // (*) sort in increasing value of soln val
1710  // track who is the biggest coef and it's index.
1711  // if biggest > adjRhs, skip row. Bad knapsack.
1712  // margin = adjRhs
1713  // idea: if (possibly compl) soln >= .5 round up, else round down
1714  // they do more, but that's the essence
1715  // go through the elements {
1716  // if round down, skip
1717  // if round up, add to element to cover. adjust margin
1718  // if current element = biggest, then get next biggest
1719  // if biggest > marg, you've got a cover. stop looking               
1720  // else try next element in the loop
1721  // }       
1722  // (*)RLH: I'm going to sort in decreasing order of soln val
1723  // b/c var's whose soln < .5 in can. form get rounded down
1724  // and skipped. If you can get a min cover of the vars
1725  // whose soln is >= .5, I believe this gives the same as J&E.
1726  // But if not, maybe I can get something more.
1727       
1728  // (**)By checking largest value left, they ensure a minimal cover
1729  // on the unsatisfied variables
1730     
1731  // if you have a cover
1732  // sort the elements to be lifted in order of their reduced costs.
1733  // lift in this order.
1734  // They lift down on variables at one, in a sequence-dependent manner.
1735  // Partion the variables into three sets: those in the cover, those
1736  // not in the cover at value one, and those remaining.
1737 
1738  fracCover.reserve(krow.getNumElements());
1739  remainder.reserve(krow.getNumElements());
1740  atOne.reserve(krow.getNumElements());
1741
1742  double unsatRhs = b;
1743
1744  // working info on unsatisfied vars
1745  CoinPackedVector unsat;
1746  unsat.reserve(krow.getNumElements());
1747
1748  // partition the (binary) variables in the canonical knapsack
1749  // into those at zero, those at fractions, and those at one.
1750  //
1751  // essentially, temporarily fix to one the free vars with lp soln value of
1752  // one by calculating the "unsatRhs". Call the result the "reduced krow".
1753  //
1754  // Note: continuous and integer vars, and variables fixed at
1755  // binary values have already been netted out
1756  // in deriving the canonical knapsack form
1757  int i;
1758  for (i=0; i<krow.getNumElements(); i++){
1759
1760    if (xstar[krow.getIndices()[i]] > onetol_){
1761      atOne.insert(krow.getIndices()[i],krow.getElements()[i]); 
1762      unsatRhs -= krow.getElements()[i];
1763    }   
1764    else if (xstar[krow.getIndices()[i]] >= epsilon_){
1765      unsat.insert(krow.getIndices()[i],krow.getElements()[i]) ;
1766    }
1767    else { 
1768      remainder.insert(krow.getIndices()[i],krow.getElements()[i]);
1769    }
1770  }
1771
1772  // sort the indices of the unsat var in order of decreasing solution value
1773  CoinDecrSolutionOrdered decrSol(xstar);
1774  unsat.sort(decrSol);
1775 
1776#ifdef CGL_DEBUG
1777  // sanity check
1778  for (i=1; i<unsat.getNumElements(); i++){
1779    double xstarim1 = xstar[unsat.getIndices()[i-1]];
1780    double xstari=  xstar[unsat.getIndices()[i]];
1781    assert( xstarim1 >= xstari );
1782  }
1783#endif
1784
1785  // get the largest coefficient among the unsatisfied variables   
1786  double bigCoef= 0.0;
1787  // double temp;
1788  int bigIndex = 0;
1789  for (i=0; i<unsat.getNumElements(); i++){
1790    if (unsat.getElements()[i]>bigCoef){
1791      bigCoef = unsat.getElements()[i];
1792      bigIndex = i;
1793    }
1794  }
1795 
1796  // initialize
1797  i=0;
1798  double margin = unsatRhs;
1799  int gotCover=0;
1800  int j;
1801 
1802  // look in order through the unsatisfied vars which along with the
1803  // the max element defines a cover
1804  while (i<unsat.getNumElements() && !gotCover){
1805    margin -= unsat.getElements()[i];
1806   
1807    // get the biggest row coef downstream in the given order   
1808    if (i == bigIndex){
1809      bigCoef = 0.0;
1810      bigIndex = 0;
1811      for (j=i+1; j<unsat.getNumElements(); j++){
1812        double temp = unsat.getElements()[j];
1813        if (temp > bigCoef ){
1814          bigCoef = temp;
1815          bigIndex = j;
1816        }
1817      }
1818    }
1819   
1820    if (bigCoef > margin+epsilon2_) gotCover = 1;
1821    i++;         
1822  }
1823
1824  // J&E approach; get first single one element that fills the margin   
1825  if(gotCover){ 
1826    j=i;
1827    if (j<unsat.getNumElements()){ // need this "if" incase nUnsat=1   
1828      while (unsat.getElements()[j]< margin) {
1829        j++;
1830      }
1831      // switch members so that first nUnsat define the cover
1832      unsat.swap(i,j);
1833      i++;
1834    }
1835   
1836    // DEBUG: verify we have a cover over the reduced krow
1837    // (may not be violated)
1838    int nCover = i;
1839    double coverElementSum = 0.0;
1840    int k;
1841    for (k=0; k<nCover; k++){
1842      coverElementSum += unsat.getElements()[k];
1843    }
1844   
1845    // Split the unsatisfied elements into those in the "reduced krow" cover
1846    // and those not in the cover. The elements not in the cover are
1847    // considered remainders. Variables atOne are reported as atOne.
1848
1849
1850    // Test if the detected cover is violated
1851    if (coverElementSum > unsatRhs+epsilon2_){
1852      for (i=nCover; i<unsat.getNumElements(); i++) {
1853        remainder.insert(unsat.getIndices()[i],unsat.getElements()[i]);
1854      }
1855      unsat.truncate(nCover);
1856      fracCover = unsat;
1857      // cover.append(atOne);
1858
1859      // Sanity check
1860      int size = (fracCover.getNumElements() + remainder.getNumElements() +
1861                  atOne.getNumElements());
1862      int krowsize =  krow.getNumElements();
1863      assert( size == krowsize );
1864     
1865      // Sort cover in terms of knapsack row coefficients   
1866      fracCover.sortDecrElement();
1867
1868    // We have a not-necessarily-violated "reduced krow" cover inequality.
1869    // Minimal on the "reduced krow"
1870
1871#if 0
1872     
1873      double oneLessCoverElementSum =
1874         coverElementSum-fracCover.getElements()[fracCover.getNumElements()-1];
1875      while (oneLessCoverElementSum > b){
1876        // move the excess cover member into the set of remainders
1877        remainder.insert(fracCover.getIndices()[fracCover.getNumElements()-1],
1878                         fracCover.getElements()[fracCover.getNumElements()-1]);
1879        fracCover.truncate(fracCover.getNumElements()-1);
1880        oneLessCoverElementSum -=
1881           fracCover.getElements()[fracCover.getNumElements()-1];
1882      }
1883#endif
1884 
1885#ifdef PRINT_DEBUG
1886      printf("More Exactly John and Ellis:");
1887      printf(" row %i has -reduced--fractional- cover of size %i.\n",
1888             row,fracCover.getNumElements());
1889      double sumFracCover = 0.0;
1890      for (i=0; i<fracCover.getNumElements(); i++){
1891        printf("index %i, element %g, xstar value % g \n",
1892               fracCover.getIndices()[i],fracCover.getElements()[i],
1893               xstar[fracCover.getIndices()[i]]);
1894        sumFracCover += fracCover.getElements()[i];
1895      }
1896      double sumAtOne = 0.0;
1897      printf("There are %i variables at one:\n",atOne.getNumElements());
1898      for (i=0; i<atOne.getNumElements(); i++){
1899        printf("index %i, element %g, xstar value % g \n",
1900               atOne.getIndices()[i],atOne.getElements()[i],
1901               xstar[atOne.getIndices()[i]]);
1902        sumAtOne += atOne.getElements()[i];
1903      }
1904      printf("The b = %g, sumAtOne = %g, unsatRhs = b-sumAtOne = %g, ",
1905             b, sumAtOne, unsatRhs);
1906      printf("and the (fractional) cover element sum is %g\n\n", sumFracCover);
1907#endif
1908
1909    }
1910    // if minimal cover is not violated, turn gotCover off
1911    else {
1912//    printf("heuristically found minimal cover is NOT violated by current lp solution");
1913      gotCover = 0;
1914    }         
1915  }
1916
1917  // If no minimal cover was found, pack it in   
1918  //  if (!gotCover || cover.getNumElements() < 2) {
1919  if (!gotCover) {
1920    return -1;
1921  }
1922 
1923  return  1;
1924}
1925
1926//-------------------------------------------------------------------
1927// findGreedyCover: attempts to find a violated minimal
1928//                  cover using a greedy approach
1929//
1930// If a cover is found, it the cover and the remainder are
1931// sorted in nonincreasing order of the coefficients.
1932//-------------------------------------------------------------------
1933int
1934CglKnapsackCover::findGreedyCover(
1935     int row,
1936     CoinPackedVector & krow,
1937     double & b,
1938     double * xstar,
1939     CoinPackedVector & cover,
1940     CoinPackedVector & remainder
1941     ) const
1942  // the row argument is a hold over from debugging
1943  // ToDo: move the print cover statement out to the mainprogram
1944  // and remove the row argument
1945
1946{ 
1947  int i;
1948  int gotCover =0;
1949 
1950  cover.reserve(krow.getNumElements());
1951  remainder.reserve(krow.getNumElements());
1952 
1953  // sort knapsack in non-increasing size of row Coefficients
1954  krow.sortDecrElement(); 
1955 
1956  // greedily pack them in
1957  // looking only at unsatisfied vars, i.e. 0<xstar[.]<1
1958  double greedyElementSum = 0.0;
1959  double greedyXstarSum = 0.0;
1960 
1961  for (i=0;i<krow.getNumElements();i++){
1962    // if xstar fractional && no cover yet, consider it for the cover
1963    if (xstar[krow.getIndices()[i]] >= epsilon_ &&
1964        xstar[krow.getIndices()[i]] <= onetol_ &&
1965        !gotCover){
1966      greedyElementSum += krow.getElements()[i];
1967      greedyXstarSum += xstar[krow.getIndices()[i]];
1968      cover.insert(krow.getIndices()[i],krow.getElements()[i]);
1969      if (greedyElementSum > b+epsilon2_){
1970        gotCover = 1;
1971      }
1972    }
1973    else{
1974      remainder.insert(krow.getIndices()[i],krow.getElements()[i]);
1975    }
1976  }
1977 
1978  // sanity check
1979  int size =  remainder.getNumElements()+cover.getNumElements();
1980  int krowsize = krow.getNumElements();
1981  assert( size==krowsize );
1982 
1983  // if no violated minimal cover was found, pack it in
1984  if ( (greedyXstarSum<=(cover.getNumElements()-1)+epsilon2_) ||
1985       (!gotCover) ||
1986       (cover.getNumElements() < 2)){
1987    return -1;
1988  }
1989 
1990#ifdef PRINT_DEBUG
1991  printf("Greedy cover: row %i has cover of size %i\n",
1992         row,cover.getNumElements());
1993  for (i=0; i<cover.getNumElements(); i++){
1994    printf("index %i, element %g, xstar value % g \n", 
1995           cover.getIndices()[i], cover.getElements()[i],
1996           xstar[cover.getIndices()[i]]);
1997  }
1998  printf("The b = %g, and the cover sum is %g\n\n", b, greedyElementSum);
1999#endif
2000
2001  return  1;
2002}
2003
2004//-------------------------------------------------------------
2005// Lift Up, Down, and Uncomplement. Add the resutling cut to the cutset
2006//
2007// In the solution to the lp relaxtion,
2008// the binary variable's solution value is either 0, 1 or fractional.
2009//
2010// Input:
2011// The variables in fracCover form a cover, when the vars atOne take value one.
2012// A cover for the krow would consist of the union of the fracCover and atOne vars
2013// (which may not be violated, and may need to be processed to acheieve minimal-ness).
2014//
2015// Rather than take the "union" cover and lift up the remainder variables,
2016// we do something a little bit more interesting with the vars at one.
2017//
2018// The ip theory says that the lifted minimal cover cut can be strengthen by
2019// "lifting down" the vars atOne.
2020// -- this is what I believe John&Ellis were doing in OSL's knapsack cover cuts
2021//    with a lifting heuristic.
2022//
2023//-------------------------------------------------------------------
2024void 
2025CglKnapsackCover::liftUpDownAndUncomplementAndAdd(
2026         int nCols,
2027         double * xstar, 
2028         int * complement,
2029         int row,
2030         int nRowElem,
2031         double & b,
2032
2033         // the following 3 packed vectors partition the krow:
2034
2035         // vars have frac soln values in lp relaxation
2036         // and form cover with the vars atOne
2037         CoinPackedVector & fracCover, 
2038         // vars have soln value of 1 in lp relaxation
2039         CoinPackedVector & atOne,
2040         // and together with fracCover form minimal (?) cover.
2041         CoinPackedVector & remainder,
2042         OsiCuts & cs ) const
2043{
2044  CoinPackedVector cut;
2045
2046  // reserve storage for the cut
2047  cut.reserve(nRowElem);
2048 
2049  // the cut coefficent for the members of the cover is 1.0
2050  cut.setConstant(fracCover.getNumElements(),fracCover.getIndices(),1.0);
2051 
2052  // Preserve the cutRhs which is |C|-1, where |C| is the size of the
2053  // fracCover.
2054  double cutRhs=fracCover.getNumElements()-1;
2055
2056  // local variables
2057  // unsatRhs is the rhs for the reduced krow
2058  double unsatRhs = 0, sumAtOne = 0;
2059  int i;
2060  for (i=0; i<atOne.getNumElements(); i++){
2061    sumAtOne += atOne.getElements()[i];
2062  }
2063  unsatRhs=b-sumAtOne;
2064
2065#ifdef PRINT_DEBUG
2066  int firstFrac = fracCover.getIndices()[0];
2067  if (unsatRhs<=0.0&&fabs(xstar[firstFrac])>epsilon2_) {
2068    printf("At one %d\n",atOne.getNumElements());
2069    for (i=0; i<atOne.getNumElements(); i++){
2070      int iColumn = atOne.getIndices()[i];
2071      printf("%d %g %g\n",atOne.getIndices()[i],atOne.getElements()[i],
2072             xstar[iColumn]);
2073    }
2074    printf("frac %d\n",fracCover.getNumElements());
2075    for (i=0; i<fracCover.getNumElements(); i++){
2076      int iColumn = fracCover.getIndices()[i];
2077      printf("%d %g %g\n",fracCover.getIndices()[i],fracCover.getElements()[i],
2078             xstar[iColumn]);
2079    }
2080  }
2081#endif
2082
2083  //assert ( unsatRhs > 0 );
2084
2085  // If there is something to lift, then calculate the lifted coefficients
2086  if (unsatRhs>0.0&&(remainder.getNumElements()+atOne.getNumElements())> 0){
2087   
2088    // What order to lift?
2089    // Take the remainder vars in decreasing order of their
2090    // xstar solution value. Sort remainder in order of decreasing
2091    // xstar value.
2092    // Lift them "up"
2093    // (The lift "down" the variables atOne.
2094    CoinDecrSolutionOrdered dso1(xstar);
2095    remainder.sort(dso1);   
2096   
2097    // a is the part of krow corresponding to vars which have been lifted
2098    // alpha are the lifted coefficients with explicit storage of lifted zero
2099    // coefficients the a.getIndices() and alpha.getIndices() are identical
2100    CoinPackedVector a(fracCover);
2101    CoinPackedVector alpha;
2102    int i;
2103    for (i=0; i<fracCover.getNumElements(); i++){
2104      alpha.insert(fracCover.getIndices()[i],1.0);
2105    }
2106    // needed as an argument for exactSolveKnapsack
2107    int * x = new int[nRowElem];
2108    double psi_j=0.0;
2109   
2110    // Order alpha and a in nonincreasing order of alpha_j/a_j. 
2111    // alpha_1/a_1 >= alpha_2/a_2 >= ... >= alpha_n/a_n   
2112    // by defining this full-storage array "ratio" to be the external sort key.
2113    // right now external sort key must be full-storage.
2114
2115    double * ratio= new double[nCols];
2116    memset(ratio, 0, (nCols*sizeof(double)));
2117    double alphasize = alpha.getNumElements();
2118    double asize = a.getNumElements();
2119    assert( alphasize == asize );
2120   
2121    for (i=0; i<a.getNumElements(); i++){
2122      if (fabs(a.getElements()[i])> epsilon_ ){
2123        ratio[a.getIndices()[i]]=alpha.getElements()[i]/a.getElements()[i];
2124      }
2125      else {
2126        ratio[a.getIndices()[i]] = 0.0;
2127      }
2128    }
2129
2130    CoinDecrSolutionOrdered dso2(ratio);
2131    a.sort(dso2);   
2132    alpha.sort(dso2);
2133   
2134#ifdef CGL_DEBUG
2135    // sanity check
2136    for ( i=1; i<a.getNumElements(); i++ ) {
2137      int alphai=  alpha.getIndices()[i];
2138      int ai = a.getIndices()[i];
2139      assert( alphai == ai);
2140    } 
2141#endif 
2142   
2143    // Loop through the remainder variables to be lifted "up", and lift.
2144    int j;
2145    for (j=0; j<remainder.getNumElements(); j++){
2146      // calculate the lifted coefficient of x_j = cutRhs-psi_j
2147      // where
2148      // psi_j =  max of the current lefthand side of the cut
2149      //          s.t. the reduced knapsack corresponding to vars that have
2150      //          been lifted <= unsatRhs-a_j
2151     
2152      // Note: For exact solve, must be sorted in
2153      // alpha_1/a_1>=alpha_2/a_2>=...>=alpha_n/a_n order
2154      // check if lifted var can take value 1
2155      if (unsatRhs - remainder.getElements()[j] < epsilon_){
2156        psi_j = cutRhs;
2157      }
2158      else {     
2159        exactSolveKnapsack(alpha.getNumElements(),
2160                 unsatRhs-remainder.getElements()[j],
2161                 alpha.getElements(),a.getElements(),psi_j,x);
2162      }
2163     
2164      // assert the new coefficient is nonegative?
2165      alpha.insert(remainder.getIndices()[j],cutRhs-psi_j);
2166      a.insert(remainder.getIndices()[j],remainder.getElements()[j]);
2167     
2168      // if the lifted coefficient is non-zero
2169      // (i.e. psi_j != cutRhs), add it to the cut
2170      if (fabs(cutRhs-psi_j)>epsilon_)
2171         cut.insert(remainder.getIndices()[j],cutRhs-psi_j);
2172     
2173      ratio[remainder.getIndices()[j]]=
2174         (cutRhs-psi_j)/remainder.getElements()[j];
2175      CoinDecrSolutionOrdered dso(ratio);
2176      a.sort(dso);   
2177      alpha.sort(dso);   
2178    }
2179
2180    // Loop throught the variables atOne and lift "down"
2181    for (j=0; j<atOne.getNumElements(); j++){
2182      // calculate the lifted coefficient of x_j = psi_j-cutRhs (now cutRhs
2183      // gets updated) where
2184      // psi_j =  max of the current lefthand side of the cut
2185      //          s.t. the reduced knapsack corresponding to vars that have
2186      //          been lifted <= unsatRhs+a_j
2187     
2188      // Note: For exact solve, must be sorted in
2189      // alpha_1/a_1>=alpha_2/a_2>=...>=alpha_n/a_n order
2190      exactSolveKnapsack(alpha.getNumElements(),
2191                         unsatRhs+atOne.getElements()[j],
2192                         alpha.getElements(),a.getElements(),psi_j,x);
2193      alpha.insert(atOne.getIndices()[j],psi_j-cutRhs);
2194      a.insert(atOne.getIndices()[j],atOne.getElements()[j]);
2195      // if the lifted coefficient is non-zero (i.e. psi_j != cutRhs), add it
2196      // to the cut
2197      if (fabs(psi_j-cutRhs)>epsilon_)
2198         cut.insert(atOne.getIndices()[j],psi_j-cutRhs);
2199
2200#ifdef CGL_DEBUG
2201      assert ( fabs(atOne.getElements()[j])> epsilon_ );
2202#else
2203      if ( fabs(atOne.getElements()[j])<= epsilon_ ) {
2204        // exit gracefully
2205        cutRhs = DBL_MAX;
2206        break;
2207      }
2208#endif
2209      ratio[atOne.getIndices()[j]]=(psi_j-cutRhs)/atOne.getElements()[j];
2210
2211      // update cutRhs and unsatRhs
2212      cutRhs = psi_j ;     
2213      unsatRhs += atOne.getElements()[j];
2214
2215      CoinDecrSolutionOrdered dso(ratio);
2216      a.sort(dso);   
2217      alpha.sort(dso);   
2218    }
2219    delete [] x;
2220    delete [] ratio;
2221  }
2222
2223  // If the cut is violated, add it to the pool
2224  // if (sum over cut.getIndices())
2225  //    cut.element()*xstar > cover.getNumElements()-1, un-complement
2226  // and add it to the pool.
2227  double sum=0;
2228  for (i=0; i<cut.getNumElements(); i++){
2229    sum+= cut.getElements()[i]*xstar[cut.getIndices()[i]];
2230  }
2231  if (sum > cutRhs+epsilon2_){
2232#ifdef PRINT_DEBUG
2233    printf("Sequentially lifted UpDown cover cut of ");
2234    printf("size %i derived from fracCover of size %i.\n",
2235           cut.getNumElements(), fracCover.getNumElements());
2236    for (i=0; i<cut.getNumElements(); i++){
2237      printf("index %i, element %g, xstar value % g \n",
2238             cut.getIndices()[i],cut.getElements()[i],
2239             xstar[cut.getIndices()[i]]);
2240    }
2241    printf("The cutRhs = %g, and the alpha_j*xstar_j sum is %g\n\n",
2242           cutRhs, sum);
2243#endif
2244   
2245    // de-complement
2246    int k;
2247    const int s = cut.getNumElements();
2248    const int * indices = cut.getIndices();
2249    double * elements = cut.getElements();
2250    for (k=0; k<s; k++){
2251      if (complement[indices[k]]) {
2252        // Negate the k'th element in packedVector cut
2253        // and correspondingly adjust the rhs
2254        elements[k] *= -1;
2255        cutRhs += elements[k];
2256      }
2257    }
2258   
2259   
2260    // Create row cut
2261    OsiRowCut rc;
2262    rc.setRow(cut);
2263#ifdef CGL_DEBUG
2264    {
2265      double * elements = cut.getElements();
2266      int * indices = cut.getIndices();
2267      int n=cut.getNumElements();
2268      for (k=0; k<n; k++){
2269        assert(indices[k]>=0);
2270        assert(elements[k]);
2271        assert (fabs(elements[k])>1.0e-12);
2272      }
2273    }
2274#endif
2275    rc.setLb(-DBL_MAX);
2276    rc.setUb(cutRhs);
2277    // ToDo: what's the default effectiveness?
2278    //  rc.setEffectiveness(1.0);
2279    // Add row cut to the cut set 
2280#ifdef PRINT_DEBUG
2281    {
2282      int k;
2283      printf("cutrhs %g %d elements\n",cutRhs,cut.getNumElements());
2284      double * elements = cut.getElements();
2285      int * indices = cut.getIndices();
2286      for (k=0; k<cut.getNumElements(); k++){
2287        printf("%d %g\n",indices[k],elements[k]);
2288      }
2289    }
2290#endif
2291    cs.insert(rc);
2292  }
2293}
2294
2295//-------------------------------------------------------------------
2296// seqLiftCoverCut:  Given a canonical knapsack inequality and a
2297//                cover, performs sequence-dependent lifting.
2298//                Reference: Nemhauser & Wolsey
2299//
2300// NW suggest a lifting heuristic order that requires an argmax operation.
2301// What's the strength vs performance tradeoff of using argmax over
2302// a heuristic that depends soley on an a-prioi ordering based on
2303// the optimal solution to the lp relaxation? ToDo:  Do both, and report.
2304//
2305//-------------------------------------------------------------------
2306void
2307CglKnapsackCover::seqLiftAndUncomplementAndAdd(
2308      int nCols,
2309      double * xstar, 
2310      int * complement,
2311      int row,                       // row index number: used for debugging
2312                                     //     and to index into row bounds
2313      int nRowElem,                  // number of elements in the row, aka row
2314                                     //     size, row length.
2315      double & b,                    // rhs of the canonical knapsack
2316                                     //     inequality derived from row
2317      CoinPackedVector & cover,       // need not be violated
2318      CoinPackedVector & remainder,
2319      OsiCuts & cs) const
2320{
2321  CoinPackedVector cut;
2322 
2323  // reserve storage for the cut
2324  cut.reserve(nRowElem);
2325 
2326  // the cut coefficent for the members of the cover is 1.0
2327  cut.setConstant(cover.getNumElements(),cover.getIndices(),1.0);
2328 
2329  // so preserve the cutRhs which is |C|-1, where |C| is the size of the cover.
2330  double cutRhs=cover.getNumElements()-1;
2331 
2332  // If there is something to lift, then calcualte the lifted coefficients
2333  if (remainder.getNumElements()> 0){
2334   
2335    // What order to lift?
2336    // Take the to-be-lifted vars in decreasing order of their
2337    // xstar solution value. Sort remainder in order of decreasing
2338    // xstar value.
2339    CoinDecrSolutionOrdered dso1(xstar);
2340    remainder.sort(dso1);   
2341   
2342    // a is the part of krow corresponding to vars which have been lifted
2343    // alpha are the lifted coefficients with explicit storage of lifted zero
2344    // coefficients the a.getIndices() and alpha.getIndices() are identical
2345    CoinPackedVector a(cover);
2346    CoinPackedVector alpha;
2347    int i;
2348    for (i=0; i<cover.getNumElements(); i++){
2349      alpha.insert(cover.getIndices()[i],1.0);
2350    }
2351    // needed as an argument for exactSolveKnapsack
2352    int * x = new int[nRowElem]; 
2353    double psi_j=0.0;
2354   
2355    // Order alpha and a in nonincreasing order of alpha_j/a_j. 
2356    // alpha_1/a_1 >= alpha_2/a_2 >= ... >= alpha_n/a_n   
2357    // by defining this full-storage array "ratio" to be the external sort key.
2358   
2359    double * ratio= new double[nCols];
2360    memset(ratio, 0, (nCols*sizeof(double)));
2361    int alphasize =  alpha.getNumElements();
2362    int asize = a.getNumElements();
2363    assert( alphasize == asize );
2364   
2365    for (i=0; i<a.getNumElements(); i++){
2366      if (fabs(a.getElements()[i])> epsilon_ ){
2367        ratio[a.getIndices()[i]]=alpha.getElements()[i]/a.getElements()[i];
2368      }
2369      else {
2370        ratio[a.getIndices()[i]] = 0.0;
2371      }
2372    }
2373   
2374    // ToDo: would be nice to have sortkey NOT be full-storage vector
2375    CoinDecrSolutionOrdered dso2(ratio);
2376    // RLH:  JP, Is there a more efficient way?
2377    // The sort is identical for a and alpha, but I'm having to sort twice
2378    // here, and at every iteration in the loop below.
2379    a.sort(dso2);   
2380    alpha.sort(dso2);
2381   
2382#ifdef CGL_DEBUG
2383    // sanity check
2384    for ( i=1; i<a.getNumElements(); i++ ) {
2385      int alphai= alpha.getIndices()[i];
2386      int ai = a.getIndices()[i];
2387      assert( alphai == ai);
2388    } 
2389#endif 
2390   
2391    // Loop through the variables to be lifted, and lift.
2392    int j;
2393    for (j=0; j<remainder.getNumElements(); j++){
2394      // calculate the lifted coefficient of x_j = cutRhs-psi_j
2395      // where psi_j =  max of the current lefthand side of the cut
2396      // s.t. the knapsack corresponding to vars that have been lifted <= b-a_j
2397     
2398      // Note: For exact solve, must be sorted in
2399      // alpha_1/a_1>=alpha_2/a_2>=...>=alpha_n/a_n order
2400      exactSolveKnapsack(alpha.getNumElements(),
2401                         b-remainder.getElements()[j],
2402                         alpha.getElements(),a.getElements(),psi_j,x);
2403      alpha.insert(remainder.getIndices()[j],cutRhs-psi_j);
2404      a.insert(remainder.getIndices()[j],remainder.getElements()[j]);
2405      // if the lifted coefficient is non-zero (i.e. psi_j != cutRhs), add it
2406      // to the cut
2407      if (fabs(cutRhs-psi_j)>epsilon_)
2408         cut.insert(remainder.getIndices()[j],cutRhs-psi_j);
2409     
2410      ratio[remainder.getIndices()[j]]=
2411         (cutRhs-psi_j)/remainder.getElements()[j];
2412      CoinDecrSolutionOrdered dso(ratio);
2413      a.sort(dso);   
2414      alpha.sort(dso);   
2415    }
2416    delete [] x;
2417    delete [] ratio;
2418  }
2419
2420  // If the cut is violated, add it to the pool
2421  // if (sum over cut.getIndices())
2422  //    cut.element()*xstar > cover.getNumElements()-1, un-complement
2423  // and add it to the pool.
2424  double sum=0;
2425  int i;
2426  for (i=0; i<cut.getNumElements(); i++){
2427    sum+= cut.getElements()[i]*xstar[cut.getIndices()[i]];
2428  }
2429  if (sum > cutRhs+epsilon2_){
2430#ifdef PRINT_DEBUG
2431    printf("Sequentially lifted cover cut of size %i derived from cover of size %i.\n",cut.getNumElements(), cover.getNumElements());
2432    for (i=0; i<cut.getNumElements(); i++){
2433      printf("index %i, element %g, xstar value % g \n", cut.getIndices()[i],cut.getElements()[i], xstar[cut.getIndices()[i]]);
2434    }
2435    printf("The cutRhs = %g, and the alpha_j*xstar_j sum is %g\n\n", cutRhs, sum);
2436#endif
2437   
2438    int k;
2439    const int s = cut.getNumElements();
2440    const int * indices = cut.getIndices();
2441    double * elements = cut.getElements();
2442    for (k=0; k<s; k++){
2443      if (complement[indices[k]]) {
2444        // Negate the k'th element in packedVector cut
2445        // and correspondingly adjust the rhs
2446        elements[k] *= -1;
2447        cutRhs += elements[k];
2448      }
2449    }
2450   
2451    // Create a row cut and add it to the cut set
2452    OsiRowCut rc;
2453    rc.setRow(cut);
2454#ifdef CGL_DEBUG
2455    {
2456      double * elements = cut.getElements();
2457      int * indices = cut.getIndices();
2458      int n=cut.getNumElements();
2459      for (k=0; k<n; k++){
2460        assert(indices[k]>=0);
2461        assert(elements[k]);
2462        assert (fabs(elements[k])>1.0e-12);
2463      }
2464    }
2465#endif
2466    rc.setLb(-DBL_MAX);
2467    rc.setUb(cutRhs);
2468    // ToDo: what's a meaningful effectivity?
2469    //  rc.setEffectiveness(1.0);
2470#ifdef PRINT_DEBUG
2471    {
2472      int k;
2473      printf("cutrhs %g\n",cutRhs);
2474      double * elements = cut.getElements();
2475      int * indices = cut.getIndices();
2476      for (k=0; k<cut.getNumElements(); k++){
2477        printf("%d %g\n",indices[k],elements[k]);
2478      }
2479    }
2480#endif
2481    cs.insert(rc);
2482  }
2483}
2484
2485//-------------------------------------------------------------------
2486// liftCoverCut:  Given a canonical knapsack inequality and a
2487//                cover, constructs a lift cover cut via
2488//                sequence-independent lifting.
2489//-------------------------------------------------------------------
2490int
2491CglKnapsackCover::liftCoverCut(
2492      double & b,
2493      int nRowElem,
2494      CoinPackedVector & cover,
2495      CoinPackedVector & remainder,
2496      CoinPackedVector & cut) const
2497{
2498  int i;
2499  int  goodCut=1;
2500  // Given knapsack ax <=b, and a cover (e.g. cover corresponds to {0,...,nCover-1})
2501  // coverIndices are assumed in nondecr order of coverElements   
2502  // a_0>=a_1>=...>=a_(nCover-1)   
2503
2504  // TODO: right now if the lifted coefficient is zero,
2505  // then it's still in the cut.
2506  // Should not carry explicit zero coefficients
2507
2508  // Calculate the sum of the knapsack coefficients of the cover variables
2509  double sum = cover.sum();
2510
2511  // Define lambda to be the "cover excess".
2512  // By definition, lambda > 0. If this is not the case, something's screwy. Exit gracefully.
2513  double lambda = sum-b;
2514  if (lambda < epsilon_) {
2515#ifdef PRINT_DEBUG
2516    printf("lambda < epsilon....aborting. \n");
2517    std::cout << "lambda " << lambda << " epsilon " << epsilon_ << std::endl;
2518    //abort();
2519    goodCut=0;
2520#else
2521    //std::cout << "lambda " << lambda << " exiting" << std::endl;
2522    goodCut=0;
2523#endif
2524  }
2525
2526  // mu is vector of partial sums:
2527  //   mu[i] = sum(j=0 to i) a_j where the cover is C={0,1,..,r}
2528  //   mu[0] = 0, mu[1]=a_0, mu[2]=a_0+a_1, etc.
2529  //   and C is assumed to be sorted in nondecreasing knapsack coefficient order.
2530  double * mu= new double[cover.getNumElements()+1];
2531  double * muMinusLambda= new double[cover.getNumElements()+1];
2532  memset(mu, 0, (cover.getNumElements()+1)*sizeof(double));
2533  memset(muMinusLambda, 0, (cover.getNumElements()+1)*sizeof(double));
2534 
2535  // mu[0] = 0, mu[1]= knapsack coef of cover element 0, etc.
2536  muMinusLambda[0]= -lambda;
2537  for(i=1; i<(cover.getNumElements()+1); i++){
2538    mu[i]=mu[i-1]+ cover.getElements()[i-1];
2539    muMinusLambda[i]=mu[i]-lambda;
2540  }
2541
2542  cut.reserve(nRowElem);
2543
2544  // the cut coefficent for the members of the cover is 1.0
2545  cut.setConstant(cover.getNumElements(),cover.getIndices(),1.0);
2546 
2547  // if f(z) is superadditive
2548  int h;
2549  if (muMinusLambda[1] >= cover.getElements()[1]-epsilon_){
2550    for (h=0; h<remainder.getNumElements(); h++){
2551      if (remainder.getElements()[h] <= muMinusLambda[1]){
2552        // cutCoef[nCut] is 0, so don't bother storing
2553      }   
2554      else{ 
2555        // Todo: searching is inefficient. sort not in cover...
2556        // change so that I sort remainder before the call to lift.
2557        int found=0;
2558        i=2;
2559        while (!found && i<(cover.getNumElements()+1)){
2560          if (remainder.getElements()[h] <= muMinusLambda[i]+epsilon_){
2561            bool e = cut.isExistingIndex(remainder.getIndices()[h]);
2562            assert( !e );
2563            cut.insert( remainder.getIndices()[h], i-1.0 );
2564            found=1;
2565          }
2566          i++;
2567        }
2568        if (!found) {
2569#ifdef CGL_DEBUG
2570          printf("Error: Unable to fix lifted coefficient\n");
2571          abort();
2572#else
2573          goodCut=0;
2574#endif
2575        }
2576      } // end else
2577    }// end for each j not in C
2578  } // end if f superadditive
2579
2580  // else use superadditive function g
2581  else {
2582    double * rho= new double[cover.getNumElements()];
2583    rho[0]=lambda;
2584    for (i=1; i<cover.getNumElements(); i++){
2585      rho[i]=CoinMax(0.0, cover.getElements()[i]- muMinusLambda[1]);
2586    }
2587   
2588    int h;
2589    for (h=0; h<remainder.getNumElements(); h++){
2590     
2591      int found=0; // Todo: searcing is inefficient: sort...
2592      i=0;
2593      while(!found && i<cover.getNumElements()){
2594        if (remainder.getElements()[h] <= muMinusLambda[i+1]+epsilon_){
2595          bool notE = !cut.isExistingIndex(remainder.getIndices()[h]);
2596          assert( notE );
2597          if (i)
2598            cut.insert( remainder.getIndices()[h], (double)i );
2599          found=1;
2600        }
2601        else if (remainder.getElements()[h] < muMinusLambda[i+1]+rho[i+1]){
2602          bool notE = !cut.isExistingIndex(remainder.getIndices()[h]); 
2603          assert( notE );
2604          double cutCoef = i+1 
2605              - (muMinusLambda[i+1]+rho[i+1]-remainder.getElements()[h])/rho[1];   
2606          if (fabs(cutCoef)>epsilon_)
2607            cut.insert( remainder.getIndices()[h], cutCoef );
2608          found=1;
2609        }
2610        i++;
2611      } // endwhile
2612    } // end for j not in C
2613    delete [] rho;
2614  } // end else use g
2615
2616  delete [] muMinusLambda;
2617  delete [] mu;
2618
2619  return goodCut;
2620}
2621
2622//-------------------------------------------------------------------
2623// A goto-less implementation of the Horowitz-Sahni exact solution
2624// procedure for solving knapsack problem.
2625//
2626// Reference: Martello and Toth, Knapsack Problems, Wiley, 1990, p30-31.
2627//
2628// ToDo: Implement a dynamic programming appraoch for case
2629//       of knapsacks with integral coefficients
2630//-------------------------------------------------------------------
2631int
2632CglKnapsackCover::exactSolveKnapsack(
2633       int n, 
2634       double c, 
2635       double const *pp, 
2636       double const *ww, 
2637       double & z, 
2638       int * x) const
2639{
2640  // The knapsack problem is to find:
2641
2642  // max {sum(j=1,n) p_j*x_j st. sum (j=1,n)w_j*x_j <= c, x binary}
2643
2644  // Notation:
2645  //     xhat : current solution vector
2646  //     zhat : current solution value = sum (j=1,n) p_j*xhat_j
2647  //     chat : current residual capacity = c - sum (j=1,n) w_j*xhat_j
2648  //     x    : best solution so far, n-vector.
2649  //     z    : value of the best solution so far =  sum (j=1,n) p_j*x_j
2650     
2651
2652  // Input: n, the number of variables;
2653  //        c, the rhs;
2654  //        p, n-vector of objective func. coefficients;
2655  //        w, n-vector of the row coeff.
2656
2657  // Output: z, the optimal objective function value;
2658  //         x, the optimal (binary) solution n-vector;
2659
2660  // Assumes items are sorted  p_1/w_1 >= p_2/w_2 >= ... >= p_n/w_n
2661 
2662  memset(x, 0, (n)*sizeof(int));
2663  int * xhat = new int[n+1];
2664  memset(xhat, 0, (n+1)*sizeof(int));
2665  int j;
2666
2667  // set up: adding the extra element and
2668  // accounting for the FORTRAN vs C difference in indexing arrays.
2669  double * p = new double[n+2];
2670  double * w = new double[n+2];
2671  int ii;
2672  for (ii=1; ii<n+1; ii++){
2673    p[ii]=pp[ii-1];
2674    w[ii]=ww[ii-1];
2675  }
2676
2677  // 1. initialize
2678  double zhat = 0.0;
2679  z = 0.0;
2680  double chat = c+epsilon2_;
2681  p[n+1] = 0.0;
2682  w[n+1]= DBL_MAX;
2683  j=1;
2684
2685  while (1){
2686    // 2. compute upper bound u
2687    // "find r = min {i: sum k=j,i w_k>chat};"
2688    ii=j;
2689    double wSemiSum = w[j];
2690    double pSemiSum = p[j];
2691    while (wSemiSum <= chat && ii<n+2){
2692      ii++;
2693      wSemiSum+=w[ii];
2694      pSemiSum+=p[ii];
2695    }
2696    if (ii==n+2){
2697      printf("Exceeded iterator limit. Aborting...\n");
2698      abort();
2699    }
2700    // r = ii at this point
2701    wSemiSum -= w[ii];
2702    pSemiSum -= p[ii];
2703    double u = pSemiSum + floor((chat - wSemiSum)*p[ii]/w[ii]);
2704   
2705    // "if (z >= zhat + u) goto 5: backtrack;"
2706    if (!(z >= zhat + u)) {
2707      do {
2708        // 3. perform a forward step
2709        while (w[j] <= chat){
2710          chat = chat - w[j];
2711          zhat = zhat + p[j];
2712          xhat[j] = 1;
2713          j+=1;
2714        }
2715        if (j<=n) {
2716          xhat[j]= 0;
2717          j+=1;
2718        }
2719      } while(j==n); 
2720
2721      // "if (j<n) goto 2: compute_ub;"
2722      if (j<n)
2723        continue;
2724     
2725      // 4. up date the best solution so far
2726      if (zhat > z) {
2727        z=zhat;
2728        int k;
2729        for (k=0; k<n; k++){
2730          x[k]=xhat[k+1];
2731        }
2732      }
2733      j=n;
2734      if (xhat[n] == 1){
2735        chat = chat+ w[n];
2736        zhat = zhat-p[n];
2737        xhat[n]=0;
2738      }
2739    }
2740    // 5. backtrack
2741    // "find i=max{k<j:xhat[k]=1};"
2742    int i=j-1; 
2743    while (!(xhat[i]==1) && i>0){
2744      i--;
2745    }
2746   
2747    // "if (no such i exists) return;"
2748    if (i==0){
2749      delete [] p;
2750      delete [] w;
2751      delete [] xhat;
2752      return 1;
2753    }
2754   
2755    chat = chat + w[i];
2756    zhat=zhat -p[i];
2757    xhat[i]=0;
2758    j=i+1;
2759    // "goto 2: compute_ub;"
2760  }
2761}
2762
2763//-------------------------------------------------------------------
2764// Default Constructor
2765//-------------------------------------------------------------------
2766CglKnapsackCover::CglKnapsackCover ()
2767:
2768CglCutGenerator(),
2769epsilon_(1.0e-08),
2770epsilon2_(1.0e-5),
2771onetol_(1-epsilon_),
2772maxInKnapsack_(50),
2773numRowsToCheck_(-1),
2774rowsToCheck_(0),
2775expensiveCuts_(false)
2776{
2777  // nothing to do here
2778}
2779
2780//-------------------------------------------------------------------
2781// Copy constructor
2782//-------------------------------------------------------------------
2783CglKnapsackCover::CglKnapsackCover (const CglKnapsackCover & source) :
2784   CglCutGenerator(source),
2785   epsilon_(source.epsilon_),
2786   epsilon2_(source.epsilon2_),
2787   onetol_(source.onetol_),
2788   maxInKnapsack_(source.maxInKnapsack_),
2789   numRowsToCheck_(source.numRowsToCheck_),
2790   rowsToCheck_(0),
2791   expensiveCuts_(source.expensiveCuts_)
2792{
2793   if (numRowsToCheck_ > 0) {
2794      rowsToCheck_ = new int[numRowsToCheck_];
2795      CoinCopyN(source.rowsToCheck_, numRowsToCheck_, rowsToCheck_);
2796   }
2797  // Nothing to do here
2798}
2799
2800//-------------------------------------------------------------------
2801// Clone
2802//-------------------------------------------------------------------
2803CglCutGenerator *
2804CglKnapsackCover::clone() const
2805{
2806  return new CglKnapsackCover(*this);
2807}
2808
2809//-------------------------------------------------------------------
2810// Destructor
2811//-------------------------------------------------------------------
2812CglKnapsackCover::~CglKnapsackCover ()
2813{
2814   delete[] rowsToCheck_;
2815  // Nothing to do here
2816}
2817
2818//----------------------------------------------------------------
2819// Assignment operator
2820//-------------------------------------------------------------------
2821CglKnapsackCover &
2822CglKnapsackCover::operator=(const CglKnapsackCover& rhs)
2823{
2824   if (this != &rhs) {
2825      CglCutGenerator::operator=(rhs);
2826      epsilon_=rhs.epsilon_;
2827      epsilon2_=rhs.epsilon2_;
2828      onetol_=rhs.onetol_;
2829      maxInKnapsack_=rhs.maxInKnapsack_;
2830      delete[] rowsToCheck_;
2831      numRowsToCheck_ = rhs.numRowsToCheck_;
2832      if (numRowsToCheck_ > 0) {
2833         rowsToCheck_ = new int[numRowsToCheck_];
2834         CoinCopyN(rhs.rowsToCheck_, numRowsToCheck_, rowsToCheck_);
2835      } else {
2836         rowsToCheck_ = 0;
2837      }
2838      expensiveCuts_ = rhs.expensiveCuts_;
2839   }
2840   return *this;
2841}
2842// Create C++ lines to get to current state
2843std::string
2844CglKnapsackCover::generateCpp( FILE * fp) 
2845{
2846  CglKnapsackCover other;
2847  fprintf(fp,"0#include \"CglKnapsackCover.hpp\"\n");
2848  fprintf(fp,"3  CglKnapsackCover knapsackCover;\n");
2849  if (maxInKnapsack_!=other.maxInKnapsack_)
2850    fprintf(fp,"3  knapsackCover.setMaxInKnapsack(%d);\n",maxInKnapsack_);
2851  else
2852    fprintf(fp,"4  knapsackCover.setMaxInKnapsack(%d);\n",maxInKnapsack_);
2853  if (expensiveCuts_ != other.expensiveCuts_) {
2854    if (expensiveCuts_) 
2855      fprintf(fp,"3  knapsackCover.switchOnExpensive();\n");
2856    else
2857      fprintf(fp,"3  knapsackCover.switchOffExpensive();\n");
2858  } else {
2859    if (expensiveCuts_) 
2860      fprintf(fp,"4  knapsackCover.switchOnExpensive();\n");
2861    else
2862      fprintf(fp,"4  knapsackCover.switchOffExpensive();\n");
2863  }
2864  if (getAggressiveness()!=other.getAggressiveness())
2865    fprintf(fp,"3  knapsackCover.setAggressiveness(%d);\n",getAggressiveness());
2866  else
2867    fprintf(fp,"4  knapsackCover.setAggressiveness(%d);\n",getAggressiveness());
2868  return "knapsackCover";
2869}
Note: See TracBrowser for help on using the repository browser.