Changeset 1582 for trunk/Cbc/src/CbcHeuristicRENS.cpp

Timestamp:

Jan 7, 2011 12:16:00 PM (3 years ago)

Author:

forrest

Message:

add some sos heuristics

File:

• trunk/Cbc/src/CbcHeuristicRENS.cpp (modified) (9 diffs)

trunk/Cbc/src/CbcHeuristicRENS.cpp

@@ -90 +90 @@
         return 0;
     numberTries_++;
+    double saveFractionSmall=fractionSmall_;
     OsiSolverInterface * solver = model_->solver();
 
@@ -95 +96 @@
     const int * integerVariable = model_->integerVariable();
 
-    const double * currentSolution = solver->getColSolution();
     OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
-    newSolver->resolve();
+    const double * currentSolution = newSolver->getColSolution();
+    int type = rensType_&15;
+    if (type<11)
+      newSolver->resolve();
     double direction = newSolver->getObjSense();
     double cutoff=model_->getCutoff();
-    //newSolver->getDblParam(OsiDualObjectiveLimit, cutoff);
+    newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
     //cutoff *= direction;
     double gap = cutoff - newSolver->getObjValue() * direction ;
     double tolerance;
     newSolver->getDblParam(OsiDualTolerance, tolerance) ;
-    if (gap > 0.0 || !newSolver->isProvenOptimal()) {
+    if ((gap > 0.0 || !newSolver->isProvenOptimal())&&type<11) {
       gap += 100.0 * tolerance;
       int nFix = newSolver->reducedCostFix(gap);
@@ -115 +118 @@
           << CoinMessageEol;
       }
-    } else {
+    } else if (type<11) {
       return 0; // finished?
     }
     int numberColumns = solver->getNumCols();
     double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns);
-    int type = rensType_&15;
     double djTolerance = (type!=1) ? -1.0e30 : 1.0e-4;
     const double * colLower = newSolver->getColLower();
     const double * colUpper = newSolver->getColUpper();
+    double * contribution = NULL;
     int numberFixed = 0;
     if (type==3) {
@@ -143 +146 @@
         delete basis;
       }
-    } else if (type>=5&&type<=10) {
+    } else if (type>=5&&type<=11) {
       /* 5 fix sets at one
          6 fix on dj but leave unfixed SOS slacks
@@ -149 +152 @@
          8 fix all at zero but leave unfixed SOS slacks
          9 as 8 but only fix all at zero if just one in set nonzero
-         10 as 7 but pi other way
+         10 fix all "stable" ones
+         11 layered approach
       */
       // SOS type fixing
@@ -203 +207 @@
           }
         }
+        // Just leave one slack in each set
+        {
+          const double * objective = newSolver->getObjCoefficients();
+          int * best = new int [numberRows];
+          double * cheapest = new double[numberRows];
+          for (int i=0;i<numberRows;i++) {
+            best[i]=-1;
+            cheapest[i]=COIN_DBL_MAX;
+          }
+          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+            if (colUpper[iColumn]>colLower[iColumn]) {
+              if (columnLength[iColumn]==1) {
+                CoinBigIndex j = columnStart[iColumn];
+                int iRow = row[j];
+                if (bestDj[iRow]<1.0e30) {
+                  double obj = direction*objective[iColumn];
+                  if (obj<cheapest[iRow]) {
+                    cheapest[iRow]=obj;
+                    best[iRow]=iColumn;
+                  }
+                }
+              }
+            }
+          }
+          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+            if (colUpper[iColumn]>colLower[iColumn]) {
+              if (columnLength[iColumn]==1) {
+                CoinBigIndex j = columnStart[iColumn];
+                int iRow = row[j];
+                if (bestDj[iRow]<1.0e30) {
+                  if (best[iRow]!=-1&&iColumn!=best[iRow]) {
+                    newSolver->setColUpper(iColumn,0.0);
+                  }
+                }
+              }
+            }
+          }
+          delete [] best;
+          delete [] cheapest;
+        }
         int nSOS=0;
         double * sort = new double [numberRows];
         const double * pi = newSolver->getRowPrice();
+        if (type==11) {
+          contribution = new double [numberRows];
+          for (int i=0;i<numberRows;i++) {
+            if (bestDj[i]<1.0e30) 
+              contribution[i]=0.0;
+            else
+              contribution[i]=-1.0;
+          }
+        }
         for (int i=0;i<numberRows;i++) {
           if (bestDj[i]<1.0e30) {
@@ -217 +270 @@
         }
         if (10*nSOS>8*numberRows) {
-          std::sort(sort,sort+nSOS);
-          int last = static_cast<int>(nSOS*0.9*fractionSmall_);
-          double tolerance = sort[last];
-          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
-            if (colUpper[iColumn]>colLower[iColumn]) {
-              CoinBigIndex j;
-              if (currentSolution[iColumn]<=1.0e-6||
-                  currentSolution[iColumn]>=0.999999) {
-                if (fixSets) {
-                  for (j = columnStart[iColumn];
-                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
-                    int iRow = row[j];
-                    double useDj;
-                    if (type==5) 
-                      useDj = bestDj[iRow];
-                    else if (type==7)
-                      useDj= -fabs(pi[iRow]);
-                    else
-                      useDj= fabs(pi[iRow]);
-                    if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
-                      numberFixed++;
-                      if (currentSolution[iColumn]<=1.0e-6)
-                        newSolver->setColUpper(iColumn,0.0);
-                      else if (currentSolution[iColumn]>=0.999999) 
-                        newSolver->setColLower(iColumn,1.0);
+          if (type<10) {
+            std::sort(sort,sort+nSOS);
+            int last = static_cast<int>(nSOS*0.9*fractionSmall_);
+            double tolerance = sort[last];
+            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+              if (colUpper[iColumn]>colLower[iColumn]) {
+                CoinBigIndex j;
+                if (currentSolution[iColumn]<=1.0e-6||
+                    currentSolution[iColumn]>=0.999999) {
+                  if (fixSets) {
+                    for (j = columnStart[iColumn];
+                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
+                      int iRow = row[j];
+                      double useDj;
+                      if (type==5) 
+                        useDj = bestDj[iRow];
+                      else if (type==7)
+                        useDj= -fabs(pi[iRow]);
+                      else
+                        useDj= fabs(pi[iRow]);
+                      if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
+                        numberFixed++;
+                        if (currentSolution[iColumn]<=1.0e-6)
+                          newSolver->setColUpper(iColumn,0.0);
+                        else if (currentSolution[iColumn]>=0.999999) 
+                          newSolver->setColLower(iColumn,1.0);
+                      }
+                    }
+                  } else if (columnLength[iColumn]==1) {
+                    // leave more slacks
+                    int iRow = row[columnStart[iColumn]];
+                    if (bestDj[iRow]<1.0e30) {
+                      // fake dj
+                      dj[iColumn] *= 0.000001;
+                    }
+                  } else if (type==8||type==9) {
+                    if (currentSolution[iColumn]<=1.0e-6) {
+                      if (type==8) {
+                        dj[iColumn] *= 1.0e6;
+                      } else {
+                        bool fix=false;
+                        for (j = columnStart[iColumn];
+                             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
+                          int iRow = row[j];
+                          if (bestDj[iRow]<1.0e25) {
+                            fix=true;
+                            break;
+                          }
+                        }
+                        if (fix) {
+                          dj[iColumn] *= 1.0e6;
+                        }
+                      }
+                    } else {
+                      dj[iColumn] *= 0.000001;
                     }
                   }
-                } else if (columnLength[iColumn]==1) {
-                  // leave more slacks
-                  int iRow = row[columnStart[iColumn]];
-                  if (bestDj[iRow]<1.0e30) {
-                    // fake dj
-                    dj[iColumn] *= 0.000001;
-                  }
-                } else if (type==8||type==9) {
-                  if (currentSolution[iColumn]<=1.0e-6) {
-                    if (type==8) {
-                      dj[iColumn] *= 1.0e6;
-                    } else {
-                      bool fix=false;
-                      for (j = columnStart[iColumn];
-                           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
-                        int iRow = row[j];
-                        if (bestDj[iRow]<1.0e25) {
-                          fix=true;
-                          break;
-                        }
-                      }
-                      if (fix) {
-                        dj[iColumn] *= 1.0e6;
-                      }
-                    }
-                  } else {
-                    dj[iColumn] *= 0.000001;
-                  }
                 }
               }
             }
-          }
-          if (fixSets)
+            if (fixSets)
+              djTolerance = 1.0e30;
+          } else if (type==10) {
+            double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
+            char * mark = new char [numberColumns];
+            char * nonzero = new char [numberColumns];
+            double factor=CoinMax(1.000001,fractionSmall_);
+            fractionSmall_ = 0.5;
+            // loosen up
+            for (int i=0;i<numberRows;i++) {
+              if (bestDj[i]>=1.0e30) {
+                newSolver->setRowUpper(i,factor*saveUpper[i]);
+              }
+            }
+            newSolver->resolve();
+            const double * solution = newSolver->getColSolution();
+            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+              mark[iColumn]=0;
+              nonzero[iColumn]=0;
+              if (colUpper[iColumn]>colLower[iColumn]&&
+                  solution[iColumn]>0.9999)
+                mark[iColumn]=1;
+              else if (solution[iColumn]>0.00001)
+                nonzero[iColumn]=1;
+            }
+            // slightly small
+            for (int i=0;i<numberRows;i++) {
+              if (bestDj[i]>=1.0e30) {
+                newSolver->setRowUpper(i,saveUpper[i]*0.9999);
+              }
+            }
+            newSolver->resolve();
+            int nCheck=2;
+            if (newSolver->isProvenOptimal()) {
+              for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+                if (colUpper[iColumn]>colLower[iColumn]&&
+                    solution[iColumn]>0.9999)
+                  mark[iColumn]++;
+                else if (solution[iColumn]>0.00001)
+                  nonzero[iColumn]=1;
+              }
+            } else {
+              nCheck=1;
+            }
+            // correct values
+            for (int i=0;i<numberRows;i++) {
+              if (bestDj[i]>=1.0e30) {
+                newSolver->setRowUpper(i,saveUpper[i]);
+              }
+            }
+            newSolver->resolve();
+            int nFixed=0;
+            int nFixedToZero=0;
+            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+              if (colUpper[iColumn]>colLower[iColumn]) {
+                if (solution[iColumn]>0.9999&&mark[iColumn]==nCheck) {
+                  newSolver->setColLower(iColumn,1.0);
+                  nFixed++;
+                } else if (!mark[iColumn]&&!nonzero[iColumn]&&
+                           columnLength[iColumn]>1&&solution[iColumn]<0.00001) {
+                  newSolver->setColUpper(iColumn,0.0);
+                  nFixedToZero++;
+                }
+              }
+            }
+            char line[100];
+            sprintf(line,"Heuristic %s fixed %d to one (and %d to zero)",
+                    heuristicName(),
+                    nFixed,nFixedToZero);
+            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
+              << line
+              << CoinMessageEol;
+            delete [] mark;
+            delete []nonzero;
+            delete [] saveUpper;
+            numberFixed=numberColumns;
             djTolerance = 1.0e30;
+          }
         }
         delete basis;
         delete [] sort;
         delete [] bestDj;
+        if (10*nSOS<=8*numberRows) {
+          // give up
+          delete [] contribution;
+          delete newSolver;
+          return 0;
+        }
       }
     }
     // Do dj to get right number
-    if (type==4||type==6||type>7) {
+    if (type==4||type==6||(type>7&&type<10)) {
       double * sort = new double [numberColumns];
       for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
@@ -297 +433 @@
       djTolerance = CoinMax(sort[last],1.0e-5);
       delete [] sort;
+    } else if (type==11) {
+      // Do layered in a different way
+      int numberRows = solver->getNumRows();
+      // Column copy
+      const CoinPackedMatrix * matrix = newSolver->getMatrixByCol();
+      const double * element = matrix->getElements();
+      const int * row = matrix->getIndices();
+      const CoinBigIndex * columnStart = matrix->getVectorStarts();
+      const int * columnLength = matrix->getVectorLengths();
+      int * whichRow = new int[numberRows];
+      int * whichSet = new int [numberColumns];
+      int nSOS=0;
+      for (int i=0;i<numberRows;i++) {
+        whichRow[i]=0;
+        if (!contribution[i])
+          nSOS++;
+      }
+      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+        whichSet[iColumn]=-2;
+        if (colUpper[iColumn]>colLower[iColumn]) {
+          CoinBigIndex j;
+          double sum=0.0;
+          int iSOS=-1;
+          int n=0;
+          for (j = columnStart[iColumn];
+               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
+            int iRow = row[j];
+            if (contribution[iRow]>=0.0) {
+              iSOS=iRow;
+              n++;
+            } else {
+              sum += fabs(element[j]);
+            }
+          }
+          if (n>1)
+            printf("Too many SOS entries (%d) for column %d\n",
+                   n,iColumn);
+          if (sum) {
+            assert (iSOS>=0);
+            contribution[iSOS] += sum;
+            whichRow[iSOS]++;
+            whichSet[iColumn]=iSOS;
+          } else {
+            whichSet[iColumn]=iSOS+numberRows;
+          }
+        }
+      }
+      int * chunk = new int [numberRows];
+      for (int i=0;i<numberRows;i++) {
+        chunk[i]=-1;
+        if (whichRow[i]) {
+          contribution[i]= - contribution[i]/static_cast<double>(whichRow[i]);
+        } else {
+          contribution[i] = COIN_DBL_MAX;
+        }
+        whichRow[i]=i;
+      }
+      newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
+      double * saveLower = CoinCopyOfArray(colLower,numberColumns);
+      double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
+      CoinSort_2(contribution,contribution+numberRows,whichRow);
+      // Set do nothing solution
+      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+        if(whichSet[iColumn]>=numberRows)
+          newSolver->setColLower(iColumn,1.0);
+      }
+      newSolver->resolve();
+      int nChunk = (nSOS+9)/10;
+      int nPass=0;
+      int inChunk=0;
+      for (int i=0;i<nSOS;i++) {
+        chunk[whichRow[i]]=nPass;
+        inChunk++;
+        if (inChunk==nChunk) {
+          inChunk=0;
+          // last two together
+          if (i+nChunk<nSOS)
+            nPass++;
+        }
+      }
+      // adjust
+      nPass++;
+      for (int iPass=0;iPass<nPass;iPass++) {
+        // fix last chunk and unfix this chunk
+        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
+          int iSOS = whichSet[iColumn];
+          if (iSOS>=0) {
+            if (iSOS>=numberRows)
+              iSOS-=numberRows;
+            if (chunk[iSOS]==iPass-1&&betterSolution[iColumn]>0.9999) {
+              newSolver->setColLower(iColumn,1.0);
+            } else if (chunk[iSOS]==iPass) {
+              newSolver->setColLower(iColumn,saveLower[iColumn]);
+              newSolver->setColUpper(iColumn,saveUpper[iColumn]);
+            }
+          }
+        }
+        // solve
+        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
+                                         model_->getCutoff(), "CbcHeuristicRENS");
+        if (returnCode < 0) {
+            returnCode = 0; // returned on size
+            break;
+        } else if ((returnCode&1)==0) {
+          // no good
+          break;
+        }
+      }
+      if ((returnCode&2) != 0) {
+        // could add cut
+        returnCode &= ~2;
+      }
+      delete [] chunk;
+      delete [] saveLower;
+      delete [] saveUpper;
+      delete [] whichRow;
+      delete [] whichSet;
+      delete [] contribution;
+      delete newSolver;
+      return returnCode;
     }
     
@@ -490 +746 @@
         }
     }
-
+    //delete [] whichRow;
+    //delete [] contribution;
     delete newSolver;
+    fractionSmall_ = saveFractionSmall;
     return returnCode;
 }