Changeset 1678 for trunk

Timestamp:

Jan 25, 2011 5:58:03 AM (9 years ago)

Author:

forrest

Message:

gub stuff

Location:

trunk/Clp/src

Files:

• ClpDynamicMatrix.cpp (modified) (4 diffs)
• ClpPackedMatrix.cpp (modified) (1 diff)
• ClpSimplex.cpp (modified) (2 diffs)
• ClpSimplexOther.cpp (modified) (4 diffs)
• ClpSimplexOther.hpp (modified) (1 diff)
• OsiClp/OsiClpSolverInterface.cpp (modified) (8 diffs)

trunk/Clp/src/ClpDynamicMatrix.cpp

@@ -177 +177 @@
      startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
      if (!numberGubColumns_) {
-          startColumn_ = new CoinBigIndex [1];
-          startColumn_[0] = 0;
+       //startColumn_ = new CoinBigIndex [1];
+       startColumn_[0] = 0;
      }
      CoinBigIndex numberElements = startColumn_[numberGubColumns_];
@@ -235 +235 @@
      numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
      matrix_ = originalMatrix;
+     //delete originalMatrixA;
      flags_ &= ~1;
      // resize model (matrix stays same)
+     // modify frequency
+     if (frequency>=50)
+       frequency = 50+(frequency-50)/2;
      int newRowSize = numberRows + CoinMin(numberSets_, frequency+numberRows) + 1;
      model->resize(newRowSize, numberNeeded);
@@ -903 +907 @@
 #endif
 #endif
-     return 0;
+     if (numberStaticRows_+numberActiveSets_<model->numberRows())
+       return 0;
+     else
+       return 1;
 }
 /*
@@ -2275 +2282 @@
           } else {
             // solo key
-            bool needKey;
+            bool needKey=false;
             if (numberActive) {
               if (whichKey<maximumGubColumns_) {

trunk/Clp/src/ClpPackedMatrix.cpp

@@ -4013 +4013 @@
      // Say no gaps
      flags_ &= ~2;
+     if (type_>=10)
+       return true; // gub
      if (check == 14 || check == 10) {
           if (matrix_->getNumElements() < columnStart[numberColumns]) {

trunk/Clp/src/ClpSimplex.cpp

@@ -1828 +1828 @@
 
      // Update hidden stuff e.g. effective RHS and gub
-     matrix_->updatePivot(this, oldIn, oldOut);
+     int invertNow=matrix_->updatePivot(this, oldIn, oldOut);
      objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_);
      if (handler_->logLevel() > 7) {
@@ -1985 +1985 @@
           }
           return 1;
-     } else if (factorization_->timeToRefactorize() && !dontInvert) {
+     } else if ((factorization_->timeToRefactorize() && !dontInvert)
+                ||invertNow) {
           //printf("ret after %d pivots\n",factorization_->pivots());
           return 1;

trunk/Clp/src/ClpSimplexOther.cpp

@@ -4453 +4453 @@
 // Sets basis from original
 void 
-ClpSimplexOther::setGubBasis(const ClpSimplex &original,const int * whichRows,
+ClpSimplexOther::setGubBasis(ClpSimplex &original,const int * whichRows,
                              const int * whichColumns)
 {
@@ -4470 +4470 @@
   //assert (firstOdd==numberNormal);
   double * solution = primalColumnSolution();
-  const double * originalSolution = original.primalColumnSolution();
+  double * originalSolution = original.primalColumnSolution();
   const double * upperSet = gubMatrix->upperSet();
+  // Column copy of GUB part
   int numberSets = gubMatrix->numberSets();
   const int * startSet = gubMatrix->startSets();
-  const CoinBigIndex * startColumn = gubMatrix->startColumn();
+  const CoinBigIndex * columnStart = gubMatrix->startColumn();
   const double * columnLower = gubMatrix->columnLower();
+#ifdef TRY_IMPROVE
+  const double * columnUpper = gubMatrix->columnUpper();
+  const double * lowerSet = gubMatrix->lowerSet();
+  const double * element = gubMatrix->element();
+  const int * row = gubMatrix->row();
+  bool allPositive=true;
+  double * rowActivity = new double[numberNonGub];
+  memset(rowActivity, 0, numberNonGub*sizeof(double));
+  {
+    // Non gub contribution
+    const double * element = matrix_->getElements();
+    const int * row = matrix_->getIndices();
+    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
+    const int * columnLength = matrix_->getVectorLengths();
+    for (int i=0;i<numberNormal;i++) {
+      int iColumn = whichColumns[i];
+      double value = originalSolution[iColumn];
+      if (value) {
+        for (CoinBigIndex j = columnStart[i];
+             j < columnStart[i] + columnLength[i]; j++) {
+          int iRow = row[j];
+          rowActivity[iRow] += value*element[j];
+        }
+      }
+    }
+  }
+  double * newSolution = new double [numberGubColumns];
+  int * slacks = new int [numberSets];
+  for (int i=0;i<numberSets;i++) {
+    double sum=0.0;
+    int iSlack=-1;
+    for (int j=startSet[i];j<startSet[i+1];j++) {
+      gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
+      int iColumn = whichColumns[j+numberNormal];
+      if (iColumn<numberColumns) {
+        columnIsGub[iColumn] = whichRows[numberNonGub+i];
+        double value = originalSolution[iColumn];
+        sum += value;
+        newSolution[j]=value;
+        for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
+          int iRow = row[k];
+          rowActivity[iRow] += value*element[k];
+          if (element[k] < 0.0)
+            allPositive=false;
+        }
+        if (columnStart[j]==columnStart[j+1])
+          iSlack=j;
+      } else {
+        newSolution[j]=0.0;
+        iSlack=j;
+        allPositive=false; // for now
+      }
+    }
+    slacks[i]=iSlack;
+    if (sum>upperSet[i]+1.0e-8) {
+      double gap = sum-upperSet[i];
+      if (iSlack>=0) {
+        double value=newSolution[iSlack];
+        if (value>0.0) {
+          double down = CoinMin(gap,value);
+          gap -= down;
+          sum -= down;
+          newSolution[iSlack] = value-down;
+        }
+      }
+      if (gap>1.0e-8) {
+        for (int j=startSet[i];j<startSet[i+1];j++) {
+          int iColumn = whichColumns[j+numberNormal];
+          if (newSolution[j]>0.0&&iColumn<numberColumns) {
+            double value = newSolution[j];
+            double down = CoinMin(gap,value);
+            gap -= down;
+            sum -= down;
+            newSolution[iSlack] = value-down;
+            for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
+              int iRow = row[k];
+              rowActivity[iRow] -= down*element[k];
+            }
+          }
+        }
+      }
+      assert (gap<1.0e-8);
+    } else if (sum<lowerSet[i]-1.0e-8) {
+      double gap = lowerSet[i]-sum;
+      if (iSlack>=0) {
+        double value=newSolution[iSlack];
+        if (value<columnUpper[iSlack]) {
+          double up = CoinMin(gap,columnUpper[iSlack]-value);
+          gap -= up;
+          sum += up;
+          newSolution[iSlack] = value+up;
+        }
+      }
+      if (gap>1.0e-8) {
+        for (int j=startSet[i];j<startSet[i+1];j++) {
+          int iColumn = whichColumns[j+numberNormal];
+          if (newSolution[j]<columnUpper[j]&&iColumn<numberColumns) {
+            double value = newSolution[j];
+            double up = CoinMin(gap,columnUpper[j]-value);
+            gap -= up;
+            sum += up;
+            newSolution[iSlack] = value+up;
+            for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
+              int iRow = row[k];
+              rowActivity[iRow] += up*element[k];
+            }
+          }
+        }
+      }
+      assert (gap<1.0e-8);
+    } 
+    if (fabs(sum-upperSet[i])>1.0e-7)
+      printf("Sum for set %d is %g - lower %g, upper %g\n",i,
+             sum,lowerSet[i],upperSet[i]);
+  }
+  if (allPositive) {
+    // See if we can improve solution
+    // first reduce if over
+    double * gaps = new double [numberNonGub];
+    double direction = optimizationDirection_; 
+    const double * cost = gubMatrix->cost();
+    bool over=false;
+    for (int i=0;i<numberNonGub;i++) {
+      double activity = rowActivity[i];
+      gaps[i]=0.0;
+      if (activity>rowUpper_[i]+1.0e-6) {
+        gaps[i]=activity-rowUpper_[i];
+        over=true;
+      }
+    }
+    double * weights = new double [numberGubColumns];
+    int * which = new int [numberGubColumns];
+    int * whichSet = new int [numberGubColumns];
+    if (over) {
+      int n=0;
+      for (int i=0;i<numberSets;i++) {
+        int iSlack = slacks[i];
+        if (iSlack<0||newSolution[iSlack]>upperSet[i]-1.0e-8)
+          continue;
+        double slackCost = cost[iSlack]*direction;
+        for (int j=startSet[i];j<startSet[i+1];j++) {
+          whichSet[j]=i;
+          double value = newSolution[j];
+          double thisCost = cost[j]*direction;
+          if (value>columnLower[j]&&j!=iSlack) {
+            if(thisCost<slackCost) {
+              double sum = 1.0e-30;
+              for (CoinBigIndex k = columnStart[j]; 
+                   k < columnStart[j+1] ; k++) {
+                int iRow = row[k];
+                sum += gaps[iRow]*element[k];
+              }
+              which[n]=j;
+              // big drop and small difference in cost better
+              weights[n++]=(slackCost-thisCost)/sum;
+            } else {
+              // slack better anyway
+              double move = value-columnLower[j];
+              newSolution[iSlack]=CoinMin(upperSet[i],
+                                          newSolution[iSlack]+move);
+              newSolution[j]=columnLower[j];
+              for (CoinBigIndex k = columnStart[j]; 
+                   k < columnStart[j+1] ; k++) {
+                int iRow = row[k];
+                rowActivity[iRow] -= move*element[k];
+              }
+            }
+          }
+        }
+      }
+      // sort
+      CoinSort_2(weights,weights+n,which);
+      for (int i=0;i<n;i++) {
+        int j= which[i];
+        int iSet = whichSet[j];
+        int iSlack = slacks[iSet];
+        assert (iSlack>=0);
+        double move = 0.0;
+        for (CoinBigIndex k = columnStart[j]; 
+             k < columnStart[j+1] ; k++) {
+          int iRow = row[k];
+          if(rowActivity[iRow]-rowUpper_[iRow]>move*element[k]) {
+            move = (rowActivity[iRow]-rowUpper_[iRow])/element[k];
+          }
+        }
+        move=CoinMin(move,newSolution[j]-columnLower[j]);
+        if (move) {
+          newSolution[j] -= move;
+          newSolution[iSlack] += move;
+          for (CoinBigIndex k = columnStart[j]; 
+               k < columnStart[j+1] ; k++) {
+            int iRow = row[k];
+            rowActivity[iRow] -= move*element[k];
+          }
+        }
+      }
+    }
+    delete [] whichSet;
+    delete [] which;
+    delete [] weights;
+    delete [] gaps;
+    // redo original status!
+    for (int i=0;i<numberSets;i++) {
+      int numberBasic=0;
+      int numberNewBasic=0;
+      int j1=-1;
+      int j2=-1;
+      for (int j=startSet[i];j<startSet[i+1];j++) {
+        if (newSolution[j]>columnLower[j]) {
+          numberNewBasic++;
+          j2=j;
+        }
+        int iOrig = whichColumns[j+numberNormal];
+        if (iOrig<numberColumns) {
+          if (original.getColumnStatus(iOrig)!=ClpSimplex::atLowerBound) {
+            numberBasic++;
+            j1=j;
+          }
+        } else {
+          int iSet = iOrig - numberColumns;
+          int iRow = whichRows[iSet+numberNonGub];
+          if (original.getRowStatus(iRow)==ClpSimplex::basic) {
+            numberBasic++;
+            j1=j;
+            abort();
+          }
+        }
+      }
+      if (numberBasic==1&&numberNewBasic==1&&
+          j1!=j2) {
+        int iOrig1=whichColumns[j1+numberNormal];
+        int iOrig2=whichColumns[j2+numberNormal];
+        ClpSimplex::Status status1 = original.getColumnStatus(iOrig1);
+        ClpSimplex::Status status2 = original.getColumnStatus(iOrig2);
+        originalSolution[iOrig1] = newSolution[j1];
+        originalSolution[iOrig2] = newSolution[j2];
+        original.setColumnStatus(iOrig1,status2);
+        original.setColumnStatus(iOrig2,status1);
+      }
+    }
+  }
+  delete [] newSolution;
+  delete [] slacks;
+  delete [] rowActivity;
+#else
   for (int i=0;i<numberSets;i++) {
     for (int j=startSet[i];j<startSet[i+1];j++) {
       gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
       int iColumn = whichColumns[j+numberNormal];
-      if (iColumn<numberColumns)
+      if (iColumn<numberColumns) {
         columnIsGub[iColumn] = whichRows[numberNonGub+i];
+      }
     }
   }
+#endif
   int * numberKey = new int [numberRows];
   memset(numberKey,0,numberRows*sizeof(int));
@@ -4530 +4778 @@
       int iSet = iOrig - numberColumns;
       int iRow = whichRows[iSet+numberNonGub];
-      if (original.getRowStatus(iRow)==ClpSimplex::basic) {
+      if (original.getRowStatus(iRow)==ClpSimplex::basic
+#ifdef TRY_IMPROVE
+          ||newSolution[i]>columnLower[i]+1.0e-8
+#endif
+          ) {
         assert(numberKey[iRow]);
         if (numberKey[iRow]==1)
@@ -4553 +4805 @@
       int chosen=-1;
       for (int j=startSet[i];j<startSet[i+1];j++) {
-        int length=startColumn[j+1]-startColumn[j];
+        int length=columnStart[j+1]-columnStart[j];
         int iOrig = whichColumns[j+numberNormal];
         double value;
         if (iOrig<numberColumns) {
-          value=originalSolution[iOrig]-columnLower[j];
+#ifdef TRY_IMPROVE
+          value=newSolution[j]-columnLower[j];
+#else
+          value = originalSolution[iOrig]-columnLower[j];
+#endif
           if (value>upper)
             gubMatrix->setStatus(i,ClpSimplex::atLowerBound);

trunk/Clp/src/ClpSimplexOther.hpp

@@ -197 +197 @@
                              int factorizationFrequency=50);
      /// Sets basis from original
-     void setGubBasis(const ClpSimplex &original,const int * whichRows,
+     void setGubBasis(ClpSimplex &original,const int * whichRows,
                       const int * whichColumns);
      /// Restores basis to original

trunk/Clp/src/OsiClp/OsiClpSolverInterface.cpp

@@ -57 +57 @@
   }
 #endif
-  if ((specialOptions_&2097152)!=0) {
-    resolveGub((9*modelPtr_->numberRows())/10);
-    return;
+  if ((specialOptions_&2097152)!=0||(specialOptions_&4194304)!=0) {
+    bool takeHint;
+    OsiHintStrength strength;
+    int algorithm = 0;
+    getHintParam(OsiDoDualInInitial,takeHint,strength);
+    if (strength!=OsiHintIgnore)
+      algorithm = takeHint ? -1 : 1;
+    if (algorithm>0||(specialOptions_&4194304)!=0) {
+      // Gub
+      resolveGub((9*modelPtr_->numberRows())/10);
+      return;
+    }
   }
   bool deleteSolver;
@@ -787 +796 @@
     return;
   }
-  if ((specialOptions_&2097152)!=0) {
-    resolveGub((9*modelPtr_->numberRows())/10);
-    return;
+  if ((specialOptions_&2097152)!=0||(specialOptions_&4194304)!=0) {
+    bool takeHint;
+    OsiHintStrength strength;
+    int algorithm = 0;
+    getHintParam(OsiDoDualInResolve,takeHint,strength);
+    if (strength!=OsiHintIgnore)
+      algorithm = takeHint ? -1 : 1;
+    if (algorithm>0||(specialOptions_&4194304)!=0) {
+      // Gub
+      resolveGub((9*modelPtr_->numberRows())/10);
+      return;
+    }
   }
   //void pclp(char *);
@@ -1286 +1304 @@
   OsiHintStrength strength;
   // Switch off printing if asked to
-  bool gotHint = (getHintParam(OsiDoReducePrint,takeHint,strength));
-  assert (gotHint);
+  getHintParam(OsiDoReducePrint,takeHint,strength);
   int saveMessageLevel=modelPtr_->logLevel();
   if (strength!=OsiHintIgnore&&takeHint) {
@@ -1296 +1313 @@
       modelPtr_->messageHandler()->setLogLevel(0);
   }
+  //modelPtr_->messageHandler()->setLogLevel(1);
   setBasis(basis_,modelPtr_);
   // find gub
@@ -1304 +1322 @@
   ClpSimplex * model2 = 
     static_cast<ClpSimplexOther *> (modelPtr_)->gubVersion(which,whichC,
-                                                           needed);
+                                                           needed,100);
   if (model2) {
     // move in solution
@@ -1312 +1330 @@
     ClpPrimalColumnSteepest steepest(5);
     model2->setPrimalColumnPivotAlgorithm(steepest);
-    double time1 = CoinCpuTime();
+    //double time1 = CoinCpuTime();
     model2->primal();
     //printf("Primal took %g seconds\n",CoinCpuTime()-time1);
@@ -1321 +1339 @@
     modelPtr_->primal(1);
     modelPtr_->setNumberIterations(totalIterations+modelPtr_->numberIterations());
+    delete model2;
   } else {
     modelPtr_->dual();
   }
+  delete [] which;
+  delete [] whichC;
   basis_ = getBasis(modelPtr_);
   modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
@@ -6333 +6354 @@
       }
       // set normal
-      specialOptions_ &= (2047|3*8192|15*65536|2097152);
+      specialOptions_ &= (2047|3*8192|15*65536|2097152|4194304);
       if (otherInformation!=NULL) {
         int * array = static_cast<int *> (otherInformation);