Changeset 1602 for trunk

Timestamp:

Sep 7, 2010 6:26:41 PM (9 years ago)

Author:

lou

Message:

Make simplex mode 1 (tableau access) work properly with maximisation. Note
that mode 2 (single pivot) is still minimisation only, though it shouldn't be
too hard to extend the fix.

Location:

trunk/Clp

Files:

• src/ClpLinearObjective.hpp (modified) (1 diff)
• src/OsiClp/OsiClpSolverInterface.cpp (modified) (18 diffs)
• src/OsiClp/OsiClpSolverInterface.hpp (modified) (11 diffs)
• test/OsiClpSolverInterfaceTest.cpp (modified) (1 diff)

trunk/Clp/src/ClpLinearObjective.hpp

@@ -20 +20 @@
      //@{
 
-     /** Returns gradient.  If Linear then solution may be NULL,
-         also returns an offset (to be added to current one)
-         If refresh is false then uses last solution
-         Uses model for scaling
-         includeLinear 0 - no, 1 as is, 2 as feasible
+     /** Returns objective coefficients.
+       
+       Offset is always set to 0.0. All other parameters unused.
      */
      virtual double * gradient(const ClpSimplex * model,

trunk/Clp/src/OsiClp/OsiClpSolverInterface.cpp

@@ -63 +63 @@
   double savedDualLimit=modelPtr_->dblParam_[ClpDualObjectiveLimit];
   if (fakeObjective_) {
+    // Clear (no objective, 0-1 and in B&B)
     modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()&(~128));
     // See if all with costs fixed
@@ -78 +79 @@
     }
     if (i==numberColumns) {
+      // Check (Clp fast dual)
       if ((specialOptions_&524288)==0) {
         // Set fake
@@ -84 +86 @@
         modelPtr_->dblParam_[ClpDualObjectiveLimit]=COIN_DBL_MAX;
       } else {
+        // Set (no objective, 0-1 and in B&B)
         modelPtr_->setMoreSpecialOptions(modelPtr_->moreSpecialOptions()|128);
       }
     }
   }
+  // Check (in branch and bound)
   if ((specialOptions_&1024)==0) {
     solver = new ClpSimplex(true);
@@ -2791 +2795 @@
 
 //------------------------------------------------------------------
-std::vector<double*> OsiClpSolverInterface::getDualRays(int maxNumRays,
+std::vector<double*> OsiClpSolverInterface::getDualRays(int /*maxNumRays*/,
                                                         bool fullRay) const
 {
@@ -2892 +2896 @@
   lastAlgorithm_=999;
   modelPtr_->whatsChanged_ &= (0xffff&~64);
-  CoinMemcpyN(array,modelPtr_->numberColumns(),
-                    modelPtr_->objective());
+  int n = modelPtr_->numberColumns() ;
+  if (fakeMinInSimplex_) {
+    std::transform(array,array+n,
+                   modelPtr_->objective(),std::negate<double>()) ;
+  } else {
+    CoinMemcpyN(array,n,modelPtr_->objective());
+  }
 }
 /* Set the lower bounds for all columns
@@ -3697 +3706 @@
 :
 OsiSolverInterface(),
-linearObjective_(NULL),
 rowsense_(NULL),
 rhs_(NULL),
@@ -3716 +3724 @@
 integerInformation_(NULL),
 whichRange_(NULL),
+fakeMinInSimplex_(false),
+linearObjective_(NULL),
 cleanupScaling_(0),
 specialOptions_(0x80000000),
@@ -3780 +3790 @@
 matrixByRowAtContinuous_(NULL),  
 integerInformation_(NULL),
-whichRange_(NULL)
+whichRange_(NULL),
+fakeMinInSimplex_(rhs.fakeMinInSimplex_)
 {
   //printf("in copy %x - > %x\n",&rhs,this);
@@ -3857 +3868 @@
 integerInformation_(NULL),
 whichRange_(NULL),
+fakeMinInSimplex_(false),
 cleanupScaling_(0),
 specialOptions_(0x80000000),
@@ -4962 +4974 @@
     //printf("obj a %g %g\n",modelPtr_->objectiveValue(),
     //     OsiSolverInterface::getObjValue());
-    return modelPtr_->objectiveValue();
+    if (fakeMinInSimplex_)
+      return -modelPtr_->objectiveValue() ;
+    else
+      return modelPtr_->objectiveValue();
   } else {
     return OsiSolverInterface::getObjValue();
@@ -4981 +4996 @@
   }
 #endif
-  modelPtr_->setObjectiveCoefficient(elementIndex,elementValue);
+  modelPtr_->setObjectiveCoefficient(elementIndex,
+    ((fakeMinInSimplex_)?-elementValue:elementValue));
 }
 
@@ -5284 +5300 @@
   }
 }
-/*Enables normal operation of subsequent functions.
-  This method is supposed to ensure that all typical things (like
-  reduced costs, etc.) are updated when individual pivots are executed
-  and can be queried by other methods 
+
+/*
+  Clp's copy-in/copy-out design paradigm is a challenge for the simplex modes.
+  Normal operation goes like this:
+    * startup() loads clp's work arrays, performing scaling for numerical
+      stability and compensating for max.
+    * clp solves the problem
+    * finish() unloads the work arrays into answer arrays, undoing scaling
+      and max compensation.
+  There are two solutions: undo scaling and max on demand, or make them
+  into noops. The various getBInv* methods undo scaling on demand (but
+  see special option 512) and do not need to worry about max. Other get
+  solution methods are not coded to do this, so the second approach is
+  used. For simplex modes, turn off scaling (necessary for both primal and
+  dual solutions) and temporarily convert max to min (necessary for dual
+  solution). This makes the unscaling in getBInv* superfluous, but don't
+  remove it. Arguably the better solution here would be to go through and
+  add unscaling and max compensation to the get solution methods. Look for
+  fakeMinInSimplex to see the places this propagates to.
+
+  TODO: setRowPrice never has worked properly, and I didn't try to fix it in
+        this go-round.
+
+  As of 100907, change applied to [enable|disable]Factorization (mode 1).
+  Limitation of [enable|disable]SimplexInterface (mode 2) noted in
+  documentation.  -- lh, 100907 --
+*/
+/*
+  Enables normal operation of subsequent functions.  This method is supposed
+  to ensure that all typical things (like reduced costs, etc.) are updated
+  when individual pivots are executed and can be queried by other methods
 */
 void 
@@ -5351 +5394 @@
   modelPtr_->setSolveType(1);
 }
+
+/*
+  Force scaling off. If the client thinks we're maximising, arrange it so
+  that clp sees minimisation while the client still sees maximisation. In
+  keeping with the spirit of the getBInv methods, special option 512 will
+  leave all work to the client.
+*/
 void 
 OsiClpSolverInterface::enableFactorization() const
 {
-  specialOptions_ &= ~0x80000000;
   saveData_.specialOptions_=specialOptions_;
-  int saveStatus = modelPtr_->problemStatus_;
+  // Try to preserve work regions, reuse factorization
   if ((specialOptions_&(1+8))!=1+8)
     setSpecialOptionsMutable((1+8)|specialOptions_);
+  // Are we allowed to make the output sensible to mere mortals?
+  if ((specialOptions_&512)==0) {
+    // Force scaling to off
+    saveData_.scalingFlag_ = modelPtr_->scalingFlag() ;
+    modelPtr_->scaling(0) ;
+    // Temporarily force to min but keep a copy of original objective.
+    if (getObjSense() < 0.0) {
+      fakeMinInSimplex_ = true ;
+      modelPtr_->setOptimizationDirection(1.0) ;
+      double *c = modelPtr_->objective() ;
+      int n = getNumCols() ;
+      linearObjective_ = new double[n] ;
+      CoinMemcpyN(c,n,linearObjective_) ;
+      std::transform(c,c+n,c,std::negate<double>()) ;
+    }
+  }
+  int saveStatus = modelPtr_->problemStatus_;
 #ifdef NDEBUG
   modelPtr_->startup(0);
@@ -5368 +5434 @@
 }
 
-//Undo whatever setting changes the above method had to make
+/*
+  Undo enableFactorization. Retrieve the special options and scaling and
+  remove the temporary objective used to fake minimisation in clp.
+*/
 void 
 OsiClpSolverInterface::disableFactorization() const
@@ -5378 +5447 @@
   int saveMessageLevel=modelPtr_->messageHandler()->logLevel();
   modelPtr_->messageHandler()->setLogLevel(0);
-  // Should re-do - for moment save arrays 
-  double * sol = CoinCopyOfArray(modelPtr_->columnActivity_,modelPtr_->numberColumns_);
-  double * dj = CoinCopyOfArray(modelPtr_->reducedCost_,modelPtr_->numberColumns_);
-  double * rsol = CoinCopyOfArray(modelPtr_->rowActivity_,modelPtr_->numberRows_);
-  double * dual = CoinCopyOfArray(modelPtr_->dual_,modelPtr_->numberRows_);
   modelPtr_->finish();
-  CoinMemcpyN(sol,modelPtr_->numberColumns_,modelPtr_->columnActivity_);
-  CoinMemcpyN(dj,modelPtr_->numberColumns_,modelPtr_->reducedCost_);
-  CoinMemcpyN(rsol,modelPtr_->numberRows_,modelPtr_->rowActivity_);
-  CoinMemcpyN(dual,modelPtr_->numberRows_,modelPtr_->dual_);
-  delete [] sol;
-  delete [] dj;
-  delete [] rsol;
-  delete [] dual;
   modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
-}
+  // Client asked for transforms on the way in, so back out.
+  if ((specialOptions_&512)==0) {
+    modelPtr_->scaling(saveData_.scalingFlag_) ;
+    if (fakeMinInSimplex_ == true) {
+      fakeMinInSimplex_ = false ;
+      modelPtr_->setOptimizationDirection(-1.0) ;
+      double *c = modelPtr_->objective() ;
+      int n = getNumCols() ;
+      std::transform(c,c+n,c,std::negate<double>()) ;
+      delete[] linearObjective_ ;
+    }
+  }
+}
+
+
+
 /* The following two methods may be replaced by the
    methods of OsiSolverInterface using OsiWarmStartBasis if:
@@ -5641 +5712 @@
 }
 
-//Get the reduced gradient for the cost vector c 
+/*
+  This method should not leave a permanent change in the solver. For
+  this reason, save a copy of the cost region and replace it after we've
+  calculated the duals and reduced costs.
+
+  On the good side, if we're maximising, we should negate the objective on
+  the way in and negate the duals on the way out. Since clp won't be doing
+  anything more with c, we can exploit (-1)(-1) = 1 and do nothing.
+*/
 void 
 OsiClpSolverInterface::getReducedGradient(
@@ -5648 +5727 @@
                                           const double * c) const
 {
-  assert (modelPtr_->solveType()==2);
+  //assert (modelPtr_->solveType()==2);
   // could do this faster with coding inside Clp
   // save current costs
   int numberColumns = modelPtr_->numberColumns();
   double * save = new double [numberColumns];
-  CoinMemcpyN(modelPtr_->costRegion(),numberColumns,save);
-  CoinMemcpyN(c,numberColumns,modelPtr_->costRegion());
+  double * obj = modelPtr_->costRegion();
+  CoinMemcpyN(obj,numberColumns,save);
+  // Compute new duals and reduced costs.
+  const double * columnScale = modelPtr_->columnScale();
+  if (!columnScale) {
+    CoinMemcpyN(c,numberColumns,obj) ;
+  } else {
+    // need to scale
+    for (int i=0;i<numberColumns;i++)
+      obj[i] = c[i]*columnScale[i];
+  }
   modelPtr_->computeDuals(NULL);
-  CoinMemcpyN(save,numberColumns,modelPtr_->costRegion());
+
+  // Restore previous cost vector
+  CoinMemcpyN(save,numberColumns,obj);
   delete [] save;
+
+  // Transfer results to parameters
   int numberRows = modelPtr_->numberRows();
-  CoinMemcpyN(modelPtr_->dualRowSolution(),numberRows,duals);
-  CoinMemcpyN(modelPtr_->djRegion(1),   numberColumns,columnReducedCosts);
+  const double * dualScaled = modelPtr_->dualRowSolution();
+  const double * djScaled = modelPtr_->djRegion(1);
+  if (!columnScale) {
+      CoinMemcpyN(dualScaled,numberRows,duals) ;
+      CoinMemcpyN(djScaled,numberColumns,columnReducedCosts) ;
+  } else {
+    // need to scale
+    const double * rowScale = modelPtr_->rowScale();
+    for (int i=0;i<numberRows;i++)
+      duals[i] = dualScaled[i]*rowScale[i];
+    for (int i=0;i<numberColumns;i++)
+      columnReducedCosts[i] = djScaled[i]/columnScale[i];
+  }
 }
 
@@ -6153 +6256 @@
   notOwned_=false;
   modelPtr_ = new ClpSimplex();
-  // This is also deleted by Clp --tkr 7/31/03
-  // delete linearObjective_;
   linearObjective_ = NULL;
   fillParamMaps();

trunk/Clp/src/OsiClp/OsiClpSolverInterface.hpp

@@ -63 +63 @@
   //@}
   
-  ///@name OsiSimplexInterface methods 
-  //@{
-  
-  /**Enables normal operation of subsequent functions.
-     This method is supposed to ensure that all typical things (like
-     reduced costs, etc.) are updated when individual pivots are executed
-     and can be queried by other methods 
-  */
-  virtual void enableSimplexInterface(bool doingPrimal);
-  
-  ///Undo whatever setting changes the above method had to make
-  virtual void disableSimplexInterface();
-  /** Returns 1 if can just do getBInv etc
-      2 if has all OsiSimplex methods
-      and 0 if it has none */
+  /*! @name OsiSimplexInterface methods
+      \brief Methods for the Osi Simplex API.
+
+    The current implementation should work for both minimisation and
+    maximisation in mode 1 (tableau access). In mode 2 (single pivot), only
+    minimisation is supported as of 100907.
+  */
+  //@{
+  /** \brief Simplex API capability.
+  
+    Returns
+     - 0 if no simplex API
+     - 1 if can just do getBInv etc
+     - 2 if has all OsiSimplex methods
+  */
   virtual int canDoSimplexInterface() const;
-  /** Tells solver that calls to getBInv etc are about to take place.
-      Underlying code may need mutable as this may be called from 
-      CglCut:;generateCuts which is const.  If that is too horrific then
-      each solver e.g. BCP or CBC will have to do something outside
-      main loop.
+
+  /*! \brief Enables simplex mode 1 (tableau access)
+  
+    Tells solver that calls to getBInv etc are about to take place.
+    Underlying code may need mutable as this may be called from 
+    CglCut::generateCuts which is const.  If that is too horrific then
+    each solver e.g. BCP or CBC will have to do something outside
+    main loop.
   */
   virtual void enableFactorization() const;
-  /// and stop
+
+  /*! \brief Undo any setting changes made by #enableFactorization */
   virtual void disableFactorization() const;
   
-  /** Sets up solver for repeated use by Osi interface.
-      The normal usage does things like keeping factorization around so can be used.
-      Will also do things like keep scaling and row copy of matrix if
-      matrix does not change.
-      adventure:
-      0 - safe stuff as above
-      1 - will take more risks - if it does not work then bug which will be fixed
-      2 - don't bother doing most extreme termination checks e.g. don't bother
-      re-factorizing if less than 20 iterations.
-      3 - Actually safer than 1 (mainly just keeps factorization)
-      
-      printOut - -1 always skip round common messages instead of doing some work
-      0 skip if normal defaults
-      1 leaves
-  */
-  void setupForRepeatedUse(int senseOfAdventure=0, int printOut=0);
-  /// Synchronize model (really if no cuts in tree)
-  virtual void synchronizeModel();
-  // Sleazy methods to fool const requirements (no less safe as modelPtr_ mutable)
-  void setSpecialOptionsMutable(unsigned int value) const;
-
   /** Returns true if a basis is available
       AND problem is optimal.  This should be used to see if
@@ -145 +128 @@
   virtual int setBasisStatus(const int* cstat, const int* rstat);
   
+  ///Get the reduced gradient for the cost vector c 
+  virtual void getReducedGradient(double* columnReducedCosts, 
+                                  double * duals,
+                                  const double * c) const ;
+  
+  ///Get a row of the tableau (slack part in slack if not NULL)
+  virtual void getBInvARow(int row, double* z, double * slack=NULL) const;
+  
+  /** Get a row of the tableau (slack part in slack if not NULL)
+      If keepScaled is true then scale factors not applied after so
+      user has to use coding similar to what is in this method
+  */
+  virtual void getBInvARow(int row, CoinIndexedVector * z, CoinIndexedVector * slack=NULL,
+                           bool keepScaled=false) const;
+  
+  ///Get a row of the basis inverse
+  virtual void getBInvRow(int row, double* z) const;
+  
+  ///Get a column of the tableau
+  virtual void getBInvACol(int col, double* vec) const ;
+  
+  ///Get a column of the tableau
+  virtual void getBInvACol(int col, CoinIndexedVector * vec) const ;
+  
+  /** Update (i.e. ftran) the vector passed in.
+      Unscaling is applied after - can't be applied before
+  */
+  
+  virtual void getBInvACol(CoinIndexedVector * vec) const ;
+  
+  ///Get a column of the basis inverse
+  virtual void getBInvCol(int col, double* vec) const ;
+  
+  /** Get basic indices (order of indices corresponds to the
+      order of elements in a vector retured by getBInvACol() and
+      getBInvCol()).
+  */
+  virtual void getBasics(int* index) const;
+
+  /*! \brief Enables simplex mode 2 (individual pivot control)
+
+     This method is supposed to ensure that all typical things (like
+     reduced costs, etc.) are updated when individual pivots are executed
+     and can be queried by other methods.
+  */
+  virtual void enableSimplexInterface(bool doingPrimal);
+  
+  /*! \brief Undo setting changes made by #enableSimplexInterface */
+  virtual void disableSimplexInterface();
+
   /** Perform a pivot by substituting a colIn for colOut in the basis. 
       The status of the leaving variable is given in statOut. Where
@@ -179 +212 @@
                               double& t, CoinPackedVector* dx);
   
-  ///Get the reduced gradient for the cost vector c 
-  virtual void getReducedGradient(double* columnReducedCosts, 
-                                  double * duals,
-                                  const double * c) const ;
-  
-  ///Get a row of the tableau (slack part in slack if not NULL)
-  virtual void getBInvARow(int row, double* z, double * slack=NULL) const;
-  
-  /** Get a row of the tableau (slack part in slack if not NULL)
-      If keepScaled is true then scale factors not applied after so
-      user has to use coding similar to what is in this method
-  */
-  virtual void getBInvARow(int row, CoinIndexedVector * z, CoinIndexedVector * slack=NULL,
-                           bool keepScaled=false) const;
-  
-  ///Get a row of the basis inverse
-  virtual void getBInvRow(int row, double* z) const;
-  
-  ///Get a column of the tableau
-  virtual void getBInvACol(int col, double* vec) const ;
-  
-  ///Get a column of the tableau
-  virtual void getBInvACol(int col, CoinIndexedVector * vec) const ;
-  
-  /** Update (i.e. ftran) the vector passed in.
-      Unscaling is applied after - can't be applied before
-  */
-  
-  virtual void getBInvACol(CoinIndexedVector * vec) const ;
-  
-  ///Get a column of the basis inverse
-  virtual void getBInvCol(int col, double* vec) const ;
-  
-  /** Get basic indices (order of indices corresponds to the
-      order of elements in a vector retured by getBInvACol() and
-      getBInvCol()).
-  */
-  virtual void getBasics(int* index) const;
   
   //@}
@@ -310 +305 @@
   
   //---------------------------------------------------------------------------
-  /**@name Hotstart related methods (primarily used in strong branching). <br>
+  /**@name Hotstart related methods (primarily used in strong branching).
      The user can create a hotstart (a snapshot) of the optimization process
      then reoptimize over and over again always starting from there.<br>
@@ -423 +418 @@
   /// Get pointer to array[getNumCols()] of objective function coefficients
   virtual const double * getObjCoefficients() const 
-  { return modelPtr_->objective(); }
+  { if (fakeMinInSimplex_)
+      return linearObjective_ ;
+    else
+      return modelPtr_->objective(); }
   
   /// Get objective function sense (1 for min (default), -1 for max)
   virtual double getObjSense() const 
-  { return modelPtr_->optimizationDirection(); }
+  { return ((fakeMinInSimplex_)?-modelPtr_->optimizationDirection():
+                                 modelPtr_->optimizationDirection()); }
   
   /// Return true if column is continuous
@@ -1077 +1076 @@
   /// Set fake objective
   void setFakeObjective(double * fakeObjective);
+  /*! \brief Set up solver for repeated use by Osi interface.
+
+    The normal usage does things like keeping factorization around so can be
+    used.  Will also do things like keep scaling and row copy of matrix if
+    matrix does not change.
+
+    \p senseOfAdventure:
+    - 0 - safe stuff as above
+    - 1 - will take more risks - if it does not work then bug which will be
+          fixed
+    - 2 - don't bother doing most extreme termination checks e.g. don't bother
+          re-factorizing if less than 20 iterations.
+    - 3 - Actually safer than 1 (mainly just keeps factorization)
+      
+    \p printOut
+    - -1 always skip round common messages instead of doing some work
+    -  0 skip if normal defaults
+    -  1 leaves
+  */
+  void setupForRepeatedUse(int senseOfAdventure=0, int printOut=0);
+  /// Synchronize model (really if no cuts in tree)
+  virtual void synchronizeModel();
+  /*! \brief Set special options in underlying clp solver.
+
+    Safe as const because #modelPtr_ is mutable.
+  */
+  void setSpecialOptionsMutable(unsigned int value) const;
+
   //@}
   
@@ -1199 +1226 @@
   /// Sort of lexicographic resolve
   void lexSolve();
+  //@}
+  
 protected:
-  //@}
-  
   /**@name Protected member data */
   //@{
   /// Clp model represented by this class instance
   mutable ClpSimplex * modelPtr_;
-  /// Linear objective - just points to ClpModel
-  double * linearObjective_;
+  //@}
   /**@name Cached information derived from the OSL model */
   //@{
@@ -1254 +1280 @@
   int itlimOrig_;
   
-  /// Last algorithm used , 1 = primal, 2 = dual
+  /*! \brief Last algorithm used
+  
+    Coded as
+    -    0 invalid
+    -    1 primal
+    -    2 dual
+    - -911 disaster in the algorithm that was attempted
+    -  999 current solution no longer optimal due to change in problem or
+           basis
+  */
   mutable int lastAlgorithm_;
   
@@ -1279 +1314 @@
   //std::map<OsiStrParam, ClpStrParam> strParamMap_;
   
+  /*! \brief Faking min to get proper dual solution signs in simplex API */
+  mutable bool fakeMinInSimplex_ ;
+  /*! \brief Linear objective
+  
+    Normally a pointer to the linear coefficient array in the clp objective.
+    An independent copy when #fakeMinInSimplex_ is true, because we need
+    something permanent to point to when #getObjCoefficients is called.
+  */
+  mutable double *linearObjective_;
+
   /// To save data in OsiSimplex stuff
   mutable ClpDataSave saveData_;
@@ -1308 +1353 @@
       128 Model will only change in column bounds
       256 Clean up model before hot start
-      512 Give user direct access to Clp regions in getBInvARow etc
+      512 Give user direct access to Clp regions in getBInvARow etc (i.e.,
+          do not unscale, and do not return result in getBInv parameters;
+          you have to know where to look for the answer)
       1024 Don't "borrow" model in initialSolve
       2048 Don't crunch
@@ -1421 +1468 @@
   /// Are we in trouble
   bool inTrouble_;
-      
   //@}
 };

trunk/Clp/test/OsiClpSolverInterfaceTest.cpp

@@ -1390 +1390 @@
       const double * rowScale = clp->rowScale();
       const double * columnScale = clp->columnScale();
+#     if 0
       if (!iPass)
         assert (!rowScale);
       else
         assert (rowScale); // only true for this example
+#     endif
       /* has to be exactly correct as in OsiClpsolverInterface.cpp
          (also redo each pass as may change