Changeset 1366 for trunk

Timestamp:

May 16, 2009 4:23:18 AM (12 years ago)

Author:

forrest

Message:

try and improve quadratic and barrier

Location:

trunk/Clp/src

Files:

• ClpMain.cpp (modified) (2 diffs)
• ClpPredictorCorrector.cpp (modified) (8 diffs)
• ClpSimplexNonlinear.cpp (modified) (4 diffs)
• ClpSimplexOther.cpp (modified) (3 diffs)
• ClpSolve.cpp (modified) (7 diffs)

trunk/Clp/src/ClpMain.cpp

@@ -20 +20 @@
 #include "CoinSort.hpp"
 // History since 1.0 at end
-#define CLPVERSION "1.09.00"
+#define CLPVERSION "1.10.00"
 
 #include "CoinMpsIO.hpp"
@@ -2605 +2605 @@
           //printf("possible %d - but no integers\n",numberPossible);
         }
-        if (numberPossible>1&&(numberIntegers||true)) {
+        if (numberPossible>1&&(numberIntegers||false)) {
           // 
           printf("possible %d - %d integers\n",numberPossible,numberIntegers);

trunk/Clp/src/ClpPredictorCorrector.cpp

@@ -195 +195 @@
   int lastGoodIteration=0;
   double bestObjectiveGap=COIN_DBL_MAX;
+  double bestObjective=COIN_DBL_MAX;
   int saveIteration=-1;
   bool sloppyOptimal=false;
@@ -248 +249 @@
     historyInfeasibility_[LENGTH_HISTORY-1]=complementarityGap_;
     // switch off saved if changes
-    if (saveIteration+10<numberIterations_&&
-        complementarityGap_*2.0<historyInfeasibility_[0])
-      saveIteration=-1;
+    //if (saveIteration+10<numberIterations_&&
+    //complementarityGap_*2.0<historyInfeasibility_[0])
+    //saveIteration=-1;
     lastStep = CoinMin(actualPrimalStep_,actualDualStep_);
     double goodGapChange;
+    //#define KEEP_GOING_IF_FIXED 0
+#ifndef KEEP_GOING_IF_FIXED 
+#define KEEP_GOING_IF_FIXED 10000
+#endif
     if (!sloppyOptimal) {
       goodGapChange=0.93;
     } else {
       goodGapChange=0.7;
+      if (numberFixed>KEEP_GOING_IF_FIXED)
+        goodGapChange=0.99; // make more likely to carry on
     } 
     double gapO;
@@ -275 +282 @@
         <<CoinMessageEol;
       //start saving best
-      if (gapO<bestObjectiveGap) {
+      if (gapO<bestObjectiveGap) 
+        bestObjectiveGap=gapO;
+      if (primalObjective_<bestObjective) {
         saveIteration=numberIterations_;
-        bestObjectiveGap=gapO;
+        bestObjective=primalObjective_;
         if (!savePi) {
           savePi=new double[numberRows_];
@@ -333 +342 @@
         }
       } 
+      if (complementarityGap_>0.5*checkGap&&primalObjective_>
+          bestObjective+1.0e-9&&
+          (numberIterations_>saveIteration+5||numberIterations_>100)) {
+        handler_->message(CLP_BARRIER_EXIT2,messages_)
+          <<saveIteration
+          <<CoinMessageEol;
+        break;
+      }
     } 
     if ((gapO<1.0e-6||(gapO<1.0e-4&&complementarityGap_<0.1))&&!sloppyOptimal) {
@@ -365 +382 @@
         <<complementarityGap_<<"not decreasing"
         <<CoinMessageEol;
-      if (gapO>0.75*lastGood) {
+      if (gapO>0.75*lastGood&&numberFixed<KEEP_GOING_IF_FIXED) {
         break;
       } 
@@ -831 +848 @@
   delete [] saveSU;
   if (savePi) {
-    //std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
+    std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
     CoinMemcpyN(savePi,numberRows_,dual_);
     CoinMemcpyN(savePrimal,numberTotal,solution_);
@@ -3432 +3449 @@
       } else {
         numberKilled++;
+        if (solution_[iColumn]!=lower_[iColumn]&&
+            solution_[iColumn]!=upper_[iColumn]) {
+          printf("%d %g %g %g\n",iColumn,lower_[iColumn],
+                 solution_[iColumn],upper_[iColumn]);
+        }
         diagonal_[iColumn]=0.0;
         zVec_[iColumn]=0.0;
@@ -3565 +3587 @@
     solutionNorm_=rhsNorm_;
   } 
-  double scaledRHSError=maximumRHSError_/solutionNorm_;
+  double scaledRHSError=maximumRHSError_/(solutionNorm_+10.0);
   bool dualFeasible=true;
+#if KEEP_GOING_IF_FIXED > 5 
   if (maximumBoundInfeasibility_>primalTolerance()||
       scaledRHSError>primalTolerance())
     primalFeasible=false;
+#else
+  if (maximumBoundInfeasibility_>primalTolerance()||
+      scaledRHSError>CoinMax(CoinMin(100.0*primalTolerance(),1.0e-5),
+                             primalTolerance())
+    primalFeasible=false;
+#endif
   // relax dual test if obj big and gap smallish
   double gap=fabs(primalObjective_-dualObjective_);

trunk/Clp/src/ClpSimplexNonlinear.cpp

@@ -655 +655 @@
                          array1);
     if (result) {
+      if (result==2) {
+        // does not look good
+        double currentObj;
+        double thetaObj;
+        double predictedObj;
+        objective_->stepLength(this,solution_,solution_,0.0,
+                               currentObj,thetaObj,predictedObj);
+        if (currentObj==predictedObj) {
+#ifdef CLP_INVESTIGATE
+        printf("looks bad - no change in obj %g\n",currentObj);
+#endif
+          if (factorization_->pivots())
+            result=3;
+          else
+            problemStatus_=0;
+        }
+      }
       if (result==3) 
         break; // null vector not accurate
@@ -663 +680 @@
         double predictedObj;
         objective_->stepLength(this,solution_,solution_,0.0,
-                               currentObj,predictedObj,thetaObj);
+                               currentObj,thetaObj,predictedObj);
         printf("obj %g after interior move\n",currentObj);
       }
@@ -2062 +2079 @@
       printf("out of loop after %d passes\n",nPasses);
 #endif
+    if (nPasses>=1000)
+      returnCode=2;
     bool ordinaryDj=false;
     //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta) 
@@ -2271 +2290 @@
           } else {
             returnCode=2; // do single incoming
-            //returnCode=1;
+            returnCode=1;
           }
         }

trunk/Clp/src/ClpSimplexOther.cpp

@@ -65 +65 @@
     //columnArray_[0]->checkClear();
     int iSequence = which[i];
+    if (iSequence<0) {
+      costIncreased[i] = 0.0;
+      sequenceIncreased[i] = -1;
+      costDecreased[i] = 0.0;
+      sequenceDecreased[i] = -1;
+      continue;
+    }
     double costIncrease=COIN_DBL_MAX;
     double costDecrease=COIN_DBL_MAX;
@@ -196 +203 @@
     }
   }
-  //rowArray_[0]->clear();
+  rowArray_[0]->clear();
   //rowArray_[1]->clear();
   //columnArray_[1]->clear();
-  //columnArray_[0]->clear();
+  columnArray_[0]->clear();
   //rowArray_[3]->clear();
   if (!optimizationDirection_)
@@ -3130 +3137 @@
       int nRow=0;
       double obj=0.0;
-      CoinZeroN(rowActivity,nRow);
+      CoinZeroN(rowActivity,numberRows_);
       for (iColumn=0;iColumn<numJ;iColumn++) {
         int iValue = stack[iColumn];

trunk/Clp/src/ClpSolve.cpp

@@ -1875 +1875 @@
     bool scale=false;
     bool doKKT=false;
+    bool forceFixing=false;
     if (barrierOptions&16) {
       barrierOptions &= ~16;
@@ -1886 +1887 @@
       barrierOptions &= ~256;
       presolveInCrossover=true;
+    }
+    if (barrierOptions&512) {
+      barrierOptions &= ~512;
+      forceFixing=true;
     }
     if (barrierOptions&8) {
@@ -2053 +2058 @@
       CoinMemcpyN(barrier.rowUpper(),numberRows,saveUpper+numberColumns);
     }
-    if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&
+    if (((numberFixed*20>barrier.numberRows()&&numberFixed>5000)||forceFixing)&&
         presolveInCrossover) {
       // may as well do presolve
-      barrier.fixFixed();
+      if (!forceFixing) {
+        barrier.fixFixed();
+      } else {
+        // Fix
+        int n=barrier.numberColumns();
+        double * lower = barrier.columnLower();
+        double * upper = barrier.columnUpper();
+        double * solution = barrier.primalColumnSolution();
+        int nFix=0;
+        for (int i=0;i<n;i++) {
+          if (barrier.fixedOrFree(i)&&lower[i]<upper[i]) {
+            double value = solution[i];
+            if (value<lower[i]+1.0e-6&&value-lower[i]<upper[i]-value) {
+              solution[i]=lower[i];
+              upper[i]=lower[i];
+              nFix++;
+            } else if (value>upper[i]-1.0e-6&&value-lower[i]>upper[i]-value) {
+              solution[i]=upper[i];
+              lower[i]=upper[i];
+              nFix++;
+            }
+          }
+        }
+        printf("%d columns fixed\n",nFix);
+        int nr=barrier.numberRows();
+        lower = barrier.rowLower();
+        upper = barrier.rowUpper();
+        solution = barrier.primalRowSolution();
+        nFix=0;
+        for (int i=0;i<nr;i++) {
+          if (barrier.fixedOrFree(i+n)&&lower[i]<upper[i]) {
+            double value = solution[i];
+            if (value<lower[i]+1.0e-6&&value-lower[i]<upper[i]-value) {
+              solution[i]=lower[i];
+              upper[i]=lower[i];
+              nFix++;
+            } else if (value>upper[i]-1.0e-6&&value-lower[i]>upper[i]-value) {
+              solution[i]=upper[i];
+              lower[i]=upper[i];
+              nFix++;
+            }
+          }
+        }
+        printf("%d row slacks fixed\n",nFix);
+      }
       saveModel2=model2;
       extraPresolve=true;
     } else if (numberFixed) {
       // Set fixed to bounds (may have restored earlier solution)
-      barrier.fixFixed(false);
+      if (!forceFixing) {
+        barrier.fixFixed(false);
+      } else {
+        // Fix
+        int n=barrier.numberColumns();
+        double * lower = barrier.columnLower();
+        double * upper = barrier.columnUpper();
+        double * solution = barrier.primalColumnSolution();
+        int nFix=0;
+        for (int i=0;i<n;i++) {
+          if (barrier.fixedOrFree(i)&&lower[i]<upper[i]) {
+            double value = solution[i];
+            if (value<lower[i]+1.0e-8&&value-lower[i]<upper[i]-value) {
+              solution[i]=lower[i];
+              upper[i]=lower[i];
+              nFix++;
+            } else if (value>upper[i]-1.0e-8&&value-lower[i]>upper[i]-value) {
+              solution[i]=upper[i];
+              lower[i]=upper[i];
+              nFix++;
+            } else {
+              //printf("fixcol %d %g <= %g <= %g\n",
+              //     i,lower[i],solution[i],upper[i]);
+            }
+          }
+        }
+        printf("%d columns fixed\n",nFix);
+        int nr=barrier.numberRows();
+        lower = barrier.rowLower();
+        upper = barrier.rowUpper();
+        solution = barrier.primalRowSolution();
+        nFix=0;
+        for (int i=0;i<nr;i++) {
+          if (barrier.fixedOrFree(i+n)&&lower[i]<upper[i]) {
+            double value = solution[i];
+            if (value<lower[i]+1.0e-5&&value-lower[i]<upper[i]-value) {
+              solution[i]=lower[i];
+              upper[i]=lower[i];
+              nFix++;
+            } else if (value>upper[i]-1.0e-5&&value-lower[i]>upper[i]-value) {
+              solution[i]=upper[i];
+              lower[i]=upper[i];
+              nFix++;
+            } else {
+              printf("fixrow %d %g <= %g <= %g\n",
+                     i,lower[i],solution[i],upper[i]);
+            }
+          }
+        }
+        printf("%d row slacks fixed\n",nFix);
+      }
     }
 #ifdef BORROW    
@@ -2121 +2220 @@
         model2->createStatus();
         // solve
-        model2->setPerturbation(100);
+        if (!forceFixing)
+          model2->setPerturbation(100);
         if (model2->factorizationFrequency()==200) {
           // User did not touch preset
@@ -2128 +2228 @@
 #if 1
         // throw some into basis 
-        {
+        if(!forceFixing) {
           int numberRows = model2->numberRows();
           int numberColumns = model2->numberColumns();
@@ -2169 +2269 @@
           delete [] sort;
           delete [] dsort;
-        }
-        // model2->allSlackBasis();
-        if (gap<1.0e-3*static_cast<double> (numberRows+numberColumns)) {
-          if (saveUpper) {
-            int numberRows = model2->numberRows();
-            int numberColumns = model2->numberColumns();
-            CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
-            CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
-            delete [] saveLower;
-            CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
-            CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
-            delete [] saveUpper;
-            saveLower=NULL;
-            saveUpper=NULL;
-          }
-          int numberRows = model2->numberRows();
-          int numberColumns = model2->numberColumns();
-          // just primal values pass
-          double saveScale = model2->objectiveScale();
-          model2->setObjectiveScale(1.0e-3);
-          model2->primal(2);
-          model2->setObjectiveScale(saveScale);
-          // save primal solution and copy back dual
-          CoinMemcpyN(model2->primalRowSolution(),
-                      numberRows,rowPrimal);
-          CoinMemcpyN(rowDual,
-                      numberRows,model2->dualRowSolution());
-          CoinMemcpyN(model2->primalColumnSolution(),
-                      numberColumns,columnPrimal);
-          CoinMemcpyN(columnDual,
-                      numberColumns,model2->dualColumnSolution());
-          //model2->primal(1);
-          // clean up reduced costs and flag variables
-          {
-            double * dj = model2->dualColumnSolution();
-            double * cost = model2->objective();
-            double * saveCost = new double[numberColumns];
-            CoinMemcpyN(cost,numberColumns,saveCost);
-            double * saveLower = new double[numberColumns];
-            double * lower = model2->columnLower();
-            CoinMemcpyN(lower,numberColumns,saveLower);
-            double * saveUpper = new double[numberColumns];
-            double * upper = model2->columnUpper();
-            CoinMemcpyN(upper,numberColumns,saveUpper);
-            int i;
-            double tolerance = 10.0*dualTolerance_;
-            for ( i=0;i<numberColumns;i++) {
-              if (model2->getStatus(i)==basic) {
-                dj[i]=0.0;
-              } else if (model2->getStatus(i)==atLowerBound) {
-                if (optimizationDirection_*dj[i]<tolerance) {
-                  if (optimizationDirection_*dj[i]<0.0) {
-                    //if (dj[i]<-1.0e-3)
-                    //printf("bad dj at lb %d %g\n",i,dj[i]);
-                    cost[i] -= dj[i];
-                    dj[i]=0.0;
+          // model2->allSlackBasis();
+          if (gap<1.0e-3*static_cast<double> (numberRows+numberColumns)) {
+            if (saveUpper) {
+              int numberRows = model2->numberRows();
+              int numberColumns = model2->numberColumns();
+              CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
+              CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
+              CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
+              CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
+              delete [] saveLower;
+              delete [] saveUpper;
+              saveLower=NULL;
+              saveUpper=NULL;
+            }
+            int numberRows = model2->numberRows();
+            int numberColumns = model2->numberColumns();
+            // just primal values pass
+            double saveScale = model2->objectiveScale();
+            model2->setObjectiveScale(1.0e-3);
+            model2->primal(2);
+            model2->setObjectiveScale(saveScale);
+            // save primal solution and copy back dual
+            CoinMemcpyN(model2->primalRowSolution(),
+                        numberRows,rowPrimal);
+            CoinMemcpyN(rowDual,
+                        numberRows,model2->dualRowSolution());
+            CoinMemcpyN(model2->primalColumnSolution(),
+                        numberColumns,columnPrimal);
+            CoinMemcpyN(columnDual,
+                        numberColumns,model2->dualColumnSolution());
+            //model2->primal(1);
+            // clean up reduced costs and flag variables
+            {
+              double * dj = model2->dualColumnSolution();
+              double * cost = model2->objective();
+              double * saveCost = new double[numberColumns];
+              CoinMemcpyN(cost,numberColumns,saveCost);
+              double * saveLower = new double[numberColumns];
+              double * lower = model2->columnLower();
+              CoinMemcpyN(lower,numberColumns,saveLower);
+              double * saveUpper = new double[numberColumns];
+              double * upper = model2->columnUpper();
+              CoinMemcpyN(upper,numberColumns,saveUpper);
+              int i;
+              double tolerance = 10.0*dualTolerance_;
+              for ( i=0;i<numberColumns;i++) {
+                if (model2->getStatus(i)==basic) {
+                  dj[i]=0.0;
+                } else if (model2->getStatus(i)==atLowerBound) {
+                  if (optimizationDirection_*dj[i]<tolerance) {
+                    if (optimizationDirection_*dj[i]<0.0) {
+                      //if (dj[i]<-1.0e-3)
+                      //printf("bad dj at lb %d %g\n",i,dj[i]);
+                      cost[i] -= dj[i];
+                      dj[i]=0.0;
+                    }
+                  } else {
+                    upper[i]=lower[i];
                   }
-                } else {
-                  upper[i]=lower[i];
-                }
-              } else if (model2->getStatus(i)==atUpperBound) {
-                if (optimizationDirection_*dj[i]>tolerance) {
-                  if (optimizationDirection_*dj[i]>0.0) {
-                    //if (dj[i]>1.0e-3)
-                    //printf("bad dj at ub %d %g\n",i,dj[i]);
-                    cost[i] -= dj[i];
-                    dj[i]=0.0;
+                } else if (model2->getStatus(i)==atUpperBound) {
+                  if (optimizationDirection_*dj[i]>tolerance) {
+                    if (optimizationDirection_*dj[i]>0.0) {
+                      //if (dj[i]>1.0e-3)
+                      //printf("bad dj at ub %d %g\n",i,dj[i]);
+                      cost[i] -= dj[i];
+                      dj[i]=0.0;
+                    }
+                  } else {
+                    lower[i]=upper[i];
                   }
-                } else {
-                  lower[i]=upper[i];
                 }
               }
+              // just dual values pass
+              //model2->setLogLevel(63);
+              //model2->setFactorizationFrequency(1);
+              model2->dual(2);
+              CoinMemcpyN(saveCost,numberColumns,cost);
+              delete [] saveCost;
+              CoinMemcpyN(saveLower,numberColumns,lower);
+              delete [] saveLower;
+              CoinMemcpyN(saveUpper,numberColumns,upper);
+              delete [] saveUpper;
             }
-            // just dual values pass
-            //model2->setLogLevel(63);
-            //model2->setFactorizationFrequency(1);
-            model2->dual(2);
-            CoinMemcpyN(saveCost,numberColumns,cost);
-            delete [] saveCost;
-            CoinMemcpyN(saveLower,numberColumns,lower);
-            delete [] saveLower;
-            CoinMemcpyN(saveUpper,numberColumns,upper);
-            delete [] saveUpper;
           }
           // and finish
@@ -2295 +2395 @@
     }
     if (saveUpper) {
-      int numberRows = model2->numberRows();
-      int numberColumns = model2->numberColumns();
-      CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
-      CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
+      if (!forceFixing) {
+        int numberRows = model2->numberRows();
+        int numberColumns = model2->numberColumns();
+        CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
+        CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
+        CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
+        CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
+      }
       delete [] saveLower;
-      CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
-      CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
       delete [] saveUpper;
       saveLower=NULL;