Changeset 2151 for trunk/Cbc/src/CbcCutGenerator.cpp

Timestamp:

Mar 4, 2015 10:40:35 AM (6 years ago)

Author:

forrest

Message:

allow for automatically weakening some cuts

File:

• trunk/Cbc/src/CbcCutGenerator.cpp (modified) (1 diff)

trunk/Cbc/src/CbcCutGenerator.cpp

@@ -616 +616 @@
             const OsiRowCutDebugger * debugger = solver->getRowCutDebugger();
 #endif
+            //#define WEAKEN_CUTS 1
+#ifdef WEAKEN_CUTS
+            const double * lower = solver->getColLower();
+            const double * upper = solver->getColUpper();
+            const double * solution = solver->getColSolution();
+#endif
             for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) {
                 OsiRowCut * thisCut = cs.rowCutPtr(k) ;
+#ifdef WEAKEN_CUTS
+                // weaken cut if coefficients not integer
+                
+                double lb=thisCut->lb();
+                double ub=thisCut->ub();
+                if (lb<-1.0e100||ub>1.0e100) {
+                  // normal cut
+                  CoinPackedVector rpv = thisCut->row();
+                  const int n = rpv.getNumElements();
+                  const int * indices = rpv.getIndices();
+                  const double * elements = rpv.getElements();
+                  double bound=0.0;
+                  double sum=0.0;
+                  bool integral=true;
+                  int nInteger=0;
+                  for (int k=0; k<n; k++) {
+                    double value = fabs(elements[k]);
+                    int column=indices[k];
+                    sum += value;
+                    if (value!=floor(value+0.5))
+                      integral=false;
+                    if (solver->isInteger(column)) {
+                      nInteger++;
+                      double largerBound = CoinMax(fabs(lower[column]),
+                                                   fabs(upper[column]));
+                      double solutionBound=fabs(solution[column])+10.0;
+                      bound += CoinMin(largerBound,solutionBound);
+                    }
+                  }
+#if WEAKEN_CUTS ==1
+                  // leave if all 0-1
+                  if (nInteger==bound)
+                    integral=true;
+#endif
+                  if (!integral) {
+                    double weakenBy=1.0e-7*(bound+sum);
+#if WEAKEN_CUTS>2
+                    weakenBy *= 10.0;
+#endif              
+                    if (lb<-1.0e100)
+                      thisCut->setUb(ub+weakenBy);
+                    else
+                      thisCut->setLb(lb-weakenBy);
+                  }
+                }
+#endif
 #ifdef CGL_DEBUG
                 if (debugger && debugger->onOptimalPath(*solver)) {
+#if CGL_DEBUG>1
+                    const double * optimal = debugger->optimalSolution();
+                    CoinPackedVector rpv = thisCut->row();
+                    const int n = rpv.getNumElements();
+                    const int * indices = rpv.getIndices();
+                    const double * elements = rpv.getElements();
+                    
+                    double lb=thisCut->lb();
+                    double ub=thisCut->ub();
+                    double sum=0.0;
+                    
+                    for (int k=0; k<n; k++){
+                      int column=indices[k];
+                      sum += optimal[column]*elements[k];
+                    }
+                    // is it nearly violated
+                    if (sum >ub - 1.0e-8 ||sum < lb + 1.0e-8) { 
+                      double violation=CoinMax(sum-ub,lb-sum);
+                      std::cout<<generatorName_<<" cut with "<<n
+                               <<" coefficients, nearly cuts off known solutions by "<<violation
+                               <<", lo="<<lb<<", ub="<<ub<<std::endl;
+                      for (int k=0; k<n; k++){
+                        int column=indices[k];
+                        std::cout<<"( "<<column<<" , "<<elements[k]<<" ) ";
+                        if ((k%4)==3)
+                          std::cout <<std::endl;
+                      }
+                      std::cout <<std::endl;
+                      std::cout <<"Non zero solution values are"<<std::endl;
+                      int j=0;
+                      for (int k=0; k<n; k++){
+                        int column=indices[k];
+                        if (fabs(optimal[column])>1.0e-9) {
+                          std::cout<<"( "<<column<<" , "<<optimal[column]<<" ) ";
+                          if ((j%4)==3)
+                            std::cout <<std::endl;
+                          j++;
+                        }
+                      }
+                      std::cout <<std::endl;
+                    }
+#endif
                     assert(!debugger->invalidCut(*thisCut));
                     if(debugger->invalidCut(*thisCut))