Some testing of the LP interface: bugs got fixed.
authorathos
Thu, 07 Jul 2005 15:00:04 +0000
changeset 15420219ee65ffcc
parent 1541 305ce06287c9
child 1543 a88ccf686a61
Some testing of the LP interface: bugs got fixed.
lemon/lp_base.h
lemon/lp_cplex.cc
test/lp_test.cc
     1.1 --- a/lemon/lp_base.h	Thu Jul 07 09:04:39 2005 +0000
     1.2 +++ b/lemon/lp_base.h	Thu Jul 07 15:00:04 2005 +0000
     1.3 @@ -138,19 +138,19 @@
     1.4        INFINITE = 4
     1.5      };
     1.6  
     1.7 -      ///\e The type of the investigated LP problem
     1.8 -      enum ProblemTypes {
     1.9 -	  ///Primal-dual feasible
    1.10 -	  PRIMAL_DUAL_FEASIBLE = 0,
    1.11 -	  ///Primal feasible dual infeasible
    1.12 -	  PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
    1.13 -	  ///Primal infeasible dual feasible
    1.14 -	  PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
    1.15 -	  ///Primal-dual infeasible
    1.16 -	  PRIMAL_DUAL_INFEASIBLE = 3,
    1.17 -	  ///Could not determine so far
    1.18 -	  UNKNOWN = 4
    1.19 -      };
    1.20 +    ///\e The type of the investigated LP problem
    1.21 +    enum ProblemTypes {
    1.22 +      ///Primal-dual feasible
    1.23 +      PRIMAL_DUAL_FEASIBLE = 0,
    1.24 +      ///Primal feasible dual infeasible
    1.25 +      PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
    1.26 +      ///Primal infeasible dual feasible
    1.27 +      PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
    1.28 +      ///Primal-dual infeasible
    1.29 +      PRIMAL_DUAL_INFEASIBLE = 3,
    1.30 +      ///Could not determine so far
    1.31 +      UNKNOWN = 4
    1.32 +    };
    1.33  
    1.34      ///The floating point type used by the solver
    1.35      typedef double Value;
    1.36 @@ -552,6 +552,8 @@
    1.37  
    1.38      virtual int _addCol() = 0;
    1.39      virtual int _addRow() = 0;
    1.40 +    virtual void _eraseCol(int col) = 0;
    1.41 +    virtual void _eraseRow(int row) = 0;
    1.42      virtual void _setRowCoeffs(int i, 
    1.43  			       int length,
    1.44                                 int  const * indices, 
    1.45 @@ -680,7 +682,7 @@
    1.46  
    1.47      ///\param c is the column to be modified
    1.48      ///\param e is a dual linear expression (see \ref DualExpr)
    1.49 -    ///\bug This is a temportary function. The interface will change to
    1.50 +    ///\bug This is a temporary function. The interface will change to
    1.51      ///a better one.
    1.52      void setCol(Col c,const DualExpr &e) {
    1.53        std::vector<int> indices;
    1.54 @@ -717,8 +719,8 @@
    1.55      ///\return The created row
    1.56      Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
    1.57  
    1.58 -    ///\brief Adds several new row
    1.59 -    ///(i.e a variables) at once
    1.60 +    ///\brief Add several new rows
    1.61 +    ///(i.e a constraints) at once
    1.62      ///
    1.63      ///This magic function takes a container as its argument
    1.64      ///and fills its elements
    1.65 @@ -843,6 +845,22 @@
    1.66        setRow(r,c);
    1.67        return r;
    1.68      }
    1.69 +    ///Erase a coloumn (i.e a variable) from the LP
    1.70 +
    1.71 +    ///\param c is the coloumn to be deleted
    1.72 +    ///\todo Please check this
    1.73 +    void eraseCol(Col c) {
    1.74 +      _eraseCol(cols.floatingId(c.id));
    1.75 +      cols.erase(c.id);
    1.76 +    }
    1.77 +    ///Erase a  row (i.e a constraint) from the LP
    1.78 +
    1.79 +    ///\param r is the row to be deleted
    1.80 +    ///\todo Please check this
    1.81 +    void eraseRow(Row r) {
    1.82 +      _eraseRow(rows.floatingId(r.id));
    1.83 +      rows.erase(r.id);
    1.84 +    }
    1.85  
    1.86      ///Set an element of the coefficient matrix of the LP
    1.87  
     2.1 --- a/lemon/lp_cplex.cc	Thu Jul 07 09:04:39 2005 +0000
     2.2 +++ b/lemon/lp_cplex.cc	Thu Jul 07 15:00:04 2005 +0000
     2.3 @@ -242,6 +242,7 @@
     2.4    {
     2.5      //CPX_PARAM_LPMETHOD 
     2.6      status = CPXlpopt(env, lp);
     2.7 +    //status = CPXprimopt(env, lp);
     2.8      if (status == 0){
     2.9        //We want to exclude some cases
    2.10        switch (CPXgetstat(env, lp)){
    2.11 @@ -327,11 +328,40 @@
    2.12  // ??case CPX_ABORT_DUAL_INFEAS           
    2.13  // ??case CPX_ABORT_CROSSOVER             
    2.14  // ??case CPX_INForUNBD                   
    2.15 -// ??case CPX_PIVOT                       
    2.16 +// ??case CPX_PIVOT              
    2.17 +         
    2.18 +//Some more interesting stuff:
    2.19 +
    2.20 +// CPX_PARAM_LPMETHOD  1062  int  LPMETHOD
    2.21 +// 0 Automatic 
    2.22 +// 1 Primal Simplex 
    2.23 +// 2 Dual Simplex 
    2.24 +// 3 Network Simplex 
    2.25 +// 4 Standard Barrier 
    2.26 +// Default: 0 
    2.27 +// Description: Method for linear optimization. 
    2.28 +// Determines which algorithm is used when CPXlpopt() (or "optimize" in the Interactive Optimizer) is called. Currently the behavior of the "Automatic" setting is that CPLEX simply invokes the dual simplex method, but this capability may be expanded in the future so that CPLEX chooses the method based on problem characteristics 
    2.29 +  //Hulye cplex
    2.30 +  void statusSwitch(CPXENVptr env,int& stat){
    2.31 +    int lpmethod;
    2.32 +    CPXgetintparam (env,CPX_PARAM_LPMETHOD,&lpmethod);
    2.33 +    if (lpmethod==2){
    2.34 +      if (stat==CPX_UNBOUNDED){
    2.35 +	stat=CPX_INFEASIBLE;
    2.36 +      }
    2.37 +      else{
    2.38 +	if (stat==CPX_INFEASIBLE)
    2.39 +	  stat=CPX_UNBOUNDED;
    2.40 +      }
    2.41 +    }
    2.42 +  }
    2.43  
    2.44    LpCplex::SolutionStatus LpCplex::_getPrimalStatus()
    2.45    {
    2.46 +    
    2.47      int stat = CPXgetstat(env, lp);
    2.48 +    statusSwitch(env,stat);
    2.49 +    //CPXgetstat(env, lp);
    2.50      //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL);
    2.51      switch (stat) {
    2.52      case 0:
    2.53 @@ -339,7 +369,8 @@
    2.54      case CPX_OPTIMAL://Optimal
    2.55        return OPTIMAL;
    2.56      case CPX_UNBOUNDED://Unbounded
    2.57 -      return INFINITE;
    2.58 +      return INFEASIBLE;//In case of dual simplex
    2.59 +      //return INFINITE;
    2.60      case CPX_INFEASIBLE://Infeasible 
    2.61   //    case CPX_IT_LIM_INFEAS:
    2.62  //     case CPX_TIME_LIM_INFEAS:
    2.63 @@ -348,7 +379,8 @@
    2.64  //     case CPX_ABORT_INFEAS:                
    2.65  //     case CPX_ABORT_PRIM_INFEAS:           
    2.66  //     case CPX_ABORT_PRIM_DUAL_INFEAS:      
    2.67 -      return INFEASIBLE;
    2.68 +      return INFINITE;//In case of dual simplex
    2.69 +      //return INFEASIBLE;
    2.70  //     case CPX_OBJ_LIM:                    
    2.71  //     case CPX_IT_LIM_FEAS:             
    2.72  //     case CPX_TIME_LIM_FEAS:                
    2.73 @@ -383,6 +415,7 @@
    2.74    LpCplex::SolutionStatus LpCplex::_getDualStatus()
    2.75    {
    2.76      int stat = CPXgetstat(env, lp);
    2.77 +    statusSwitch(env,stat);
    2.78      switch (stat) {
    2.79      case 0:
    2.80        return UNDEFINED; //Undefined
     3.1 --- a/test/lp_test.cc	Thu Jul 07 09:04:39 2005 +0000
     3.2 +++ b/test/lp_test.cc	Thu Jul 07 15:00:04 2005 +0000
     3.3 @@ -1,6 +1,7 @@
     3.4  #include<lemon/lp_skeleton.h>
     3.5  #include "test_tools.h"
     3.6  
     3.7 +
     3.8  #ifdef HAVE_CONFIG_H
     3.9  #include <config.h>
    3.10  #endif
    3.11 @@ -182,11 +183,26 @@
    3.12  #endif
    3.13  }
    3.14  
    3.15 +void solveAndCheck(LpSolverBase& lp, LpSolverBase::SolutionStatus stat, 
    3.16 +		   double exp_opt){
    3.17 +  lp.solve();
    3.18 +  //int decimal,sign;
    3.19 +  std::string buf1;
    3.20 +  //  itoa(stat,buf1, 10);
    3.21 +  check(lp.primalStatus()==stat,"Primalstatus should be "+buf1);
    3.22 +    
    3.23 +  if (stat ==  LpSolverBase::OPTIMAL){
    3.24 +    check(std::abs(lp.primalValue()-exp_opt)<1e-3,
    3.25 +	  "Wrong optimal value: the right optimum is ");
    3.26 +    //+ecvt(exp_opt,2)
    3.27 +  }
    3.28 +}
    3.29 + 
    3.30  void aTest(LpSolverBase & lp)
    3.31  {
    3.32    typedef LpSolverBase LP;
    3.33  
    3.34 - //The following example is taken from the book by Gáspár and Temesi, page 39.
    3.35 + //The following example is very simple
    3.36  
    3.37    typedef LpSolverBase::Row Row;
    3.38    typedef LpSolverBase::Col Col;
    3.39 @@ -197,38 +213,62 @@
    3.40  
    3.41  
    3.42    //Constraints
    3.43 -  lp.addRow(3*x1+2*x2 >=6);  
    3.44 -  lp.addRow(-1*x1+x2<=4);  
    3.45 -  lp.addRow(5*x1+8*x2<=40);  
    3.46 -  lp.addRow(x1-2*x2<=4);  
    3.47 +  Row upright=lp.addRow(x1+x2 <=1);  
    3.48 +  lp.addRow(x1+x2 >=-1);  
    3.49 +  lp.addRow(x1-x2 <=1);  
    3.50 +  lp.addRow(x1-x2 >=-1);  
    3.51    //Nonnegativity of the variables
    3.52    lp.colLowerBound(x1, 0);
    3.53    lp.colLowerBound(x2, 0);
    3.54    //Objective function
    3.55 -  lp.setObj(2*x1+x2);
    3.56 +  lp.setObj(x1+x2);
    3.57  
    3.58    lp.max();
    3.59 -  lp.solve();
    3.60  
    3.61 -  double opt=122.0/9.0;
    3.62 +  //Maximization of x1+x2
    3.63 +  //over the triangle with vertices (0,0) (0,1) (1,0)
    3.64 +  double expected_opt=1;
    3.65 +  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
    3.66    
    3.67 -  if (lp.primalStatus()==LpSolverBase::OPTIMAL){
    3.68 -    std::cout<< "Z = "<<lp.primalValue()
    3.69 -	     << " (error = " << lp.primalValue()-opt
    3.70 -	     << "); x1 = "<<lp.primal(x1)
    3.71 -	     << "; x2 = "<<lp.primal(x2)
    3.72 -	     <<std::endl;
    3.73 +  //Minimization
    3.74 +  lp.min();
    3.75 +  expected_opt=0;
    3.76 +  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
    3.77 +  
    3.78 +  //Vertex (-1,0) instead of (0,0)
    3.79 +  lp.colLowerBound(x1, -LpSolverBase::INF);
    3.80 +  expected_opt=-1;
    3.81 +  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
    3.82 +
    3.83 +  //Erase one constraint and return to maximization
    3.84 +  lp.eraseRow(upright);
    3.85 +  lp.max();
    3.86 +  expected_opt=LpSolverBase::INF;
    3.87 +  solveAndCheck(lp, LpSolverBase::INFINITE, expected_opt);
    3.88 +
    3.89 +  //Infeasibilty
    3.90 +  lp.addRow(x1+x2 <=-2);  
    3.91 +  solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt);
    3.92 +
    3.93 +  //Change problem and forget to solve
    3.94 +  lp.min();
    3.95 +  check(lp.primalStatus()==LpSolverBase::UNDEFINED,"Primalstatus should be UNDEFINED");
    3.96 +
    3.97 +//   lp.solve();
    3.98 +//   if (lp.primalStatus()==LpSolverBase::OPTIMAL){
    3.99 +//     std::cout<< "Z = "<<lp.primalValue()
   3.100 +// 	     << " (error = " << lp.primalValue()-expected_opt
   3.101 +// 	     << "); x1 = "<<lp.primal(x1)
   3.102 +// 	     << "; x2 = "<<lp.primal(x2)
   3.103 +// 	     <<std::endl;
   3.104      
   3.105 -  }
   3.106 -  else{
   3.107 -    std::cout<<"Optimal solution not found!"<<std::endl;
   3.108 -  }
   3.109 +//   }
   3.110 +//   else{
   3.111 +//     std::cout<<lp.primalStatus()<<std::endl;
   3.112 +//     std::cout<<"Optimal solution not found!"<<std::endl;
   3.113 +//   }
   3.114  
   3.115 -  check(lp.primalStatus()==LpSolverBase::OPTIMAL,"Primalstatus should be OPTIMAL");
   3.116 -
   3.117 -  check(std::abs(lp.primalValue()-opt)<1e-3,
   3.118 -	"Wrong optimal value: the right optimum is 122/9 (13.555555...)");
   3.119 -
   3.120 + 
   3.121  
   3.122  }
   3.123