test/lp_test.cc
changeset 459 ed54c0d13df0
parent 458 7afc121e0689
child 460 76ec7bd57026
equal deleted inserted replaced
0:5d99988629cd 1:170744daf33d
    35 
    35 
    36 #ifdef HAVE_SOPLEX
    36 #ifdef HAVE_SOPLEX
    37 #include <lemon/lp_soplex.h>
    37 #include <lemon/lp_soplex.h>
    38 #endif
    38 #endif
    39 
    39 
       
    40 #ifdef HAVE_CLP
       
    41 #include <lemon/lp_clp.h>
       
    42 #endif
       
    43 
    40 using namespace lemon;
    44 using namespace lemon;
    41 
    45 
    42 void lpTest(LpSolverBase & lp)
    46 void lpTest(LpSolver& lp)
    43 {
    47 {
    44 
    48 
    45 
    49   typedef LpSolver LP;
    46 
       
    47   typedef LpSolverBase LP;
       
    48 
    50 
    49   std::vector<LP::Col> x(10);
    51   std::vector<LP::Col> x(10);
    50   //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
    52   //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
    51   lp.addColSet(x);
    53   lp.addColSet(x);
    52   lp.colLowerBound(x,1);
    54   lp.colLowerBound(x,1);
    53   lp.colUpperBound(x,1);
    55   lp.colUpperBound(x,1);
    54   lp.colBounds(x,1,2);
    56   lp.colBounds(x,1,2);
    55 #ifndef GYORSITAS
       
    56 
    57 
    57   std::vector<LP::Col> y(10);
    58   std::vector<LP::Col> y(10);
    58   lp.addColSet(y);
    59   lp.addColSet(y);
    59 
    60 
    60   lp.colLowerBound(y,1);
    61   lp.colLowerBound(y,1);
    84     p3=lp.addCol();
    85     p3=lp.addCol();
    85     p4=lp.addCol();
    86     p4=lp.addCol();
    86     p5=lp.addCol();
    87     p5=lp.addCol();
    87 
    88 
    88     e[p1]=2;
    89     e[p1]=2;
    89     e.constComp()=12;
    90     *e=12;
    90     e[p1]+=2;
    91     e[p1]+=2;
    91     e.constComp()+=12;
    92     *e+=12;
    92     e[p1]-=2;
    93     e[p1]-=2;
    93     e.constComp()-=12;
    94     *e-=12;
    94 
    95 
    95     e=2;
    96     e=2;
    96     e=2.2;
    97     e=2.2;
    97     e=p1;
    98     e=p1;
    98     e=f;
    99     e=f;
   168     c = (2 >= p1>= 3);
   169     c = (2 >= p1>= 3);
   169 
   170 
   170     e[x[3]]=2;
   171     e[x[3]]=2;
   171     e[x[3]]=4;
   172     e[x[3]]=4;
   172     e[x[3]]=1;
   173     e[x[3]]=1;
   173     e.constComp()=12;
   174     *e=12;
   174 
   175 
   175     lp.addRow(LP::INF,e,23);
   176     lp.addRow(-LP::INF,e,23);
   176     lp.addRow(LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
   177     lp.addRow(-LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
   177     lp.addRow(LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
   178     lp.addRow(-LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
   178 
   179 
   179     lp.addRow(x[1]+x[3]<=x[5]-3);
   180     lp.addRow(x[1]+x[3]<=x[5]-3);
   180     lp.addRow(-7<=x[1]+x[3]-12<=3);
   181     lp.addRow(-7<=x[1]+x[3]-12<=3);
   181     lp.addRow(x[1]<=x[5]);
   182     lp.addRow(x[1]<=x[5]);
   182 
   183 
   183     std::ostringstream buf;
   184     std::ostringstream buf;
   184 
       
   185 
       
   186     //Checking the simplify function
       
   187 
       
   188 //     //How to check the simplify function? A map gives no information
       
   189 //     //on the question whether a given key is or is not stored in it, or
       
   190 //     //it does?
       
   191 //   Yes, it does, using the find() function.
       
   192     e=((p1+p2)+(p1-p2));
       
   193     e.simplify();
       
   194     buf << "Coeff. of p2 should be 0";
       
   195     //    std::cout<<e[p1]<<e[p2]<<e[p3]<<std::endl;
       
   196     check(e.find(p2)==e.end(), buf.str());
       
   197 
       
   198 
       
   199 
   185 
   200 
   186 
   201     e=((p1+p2)+(p1-0.99*p2));
   187     e=((p1+p2)+(p1-0.99*p2));
   202     //e.prettyPrint(std::cout);
   188     //e.prettyPrint(std::cout);
   203     //(e<=2).prettyPrint(std::cout);
   189     //(e<=2).prettyPrint(std::cout);
   207     check(e[p2]>0, buf.str());
   193     check(e[p2]>0, buf.str());
   208 
   194 
   209     tolerance=0.02;
   195     tolerance=0.02;
   210     e.simplify(tolerance);
   196     e.simplify(tolerance);
   211     buf << "Coeff. of p2 should be 0";
   197     buf << "Coeff. of p2 should be 0";
   212     check(e.find(p2)==e.end(), buf.str());
   198     check(const_cast<const LpSolver::Expr&>(e)[p2]==0, buf.str());
   213 
   199 
   214 
   200 
   215   }
   201   }
   216 
   202 
   217   {
   203   {
   245        2.2*p1+p1*2.2+p1/2.2+
   231        2.2*p1+p1*2.2+p1/2.2+
   246        2*p1+p1*2+p1/2
   232        2*p1+p1*2+p1/2
   247        );
   233        );
   248   }
   234   }
   249 
   235 
   250 #endif
       
   251 }
   236 }
   252 
   237 
   253 void solveAndCheck(LpSolverBase& lp, LpSolverBase::SolutionStatus stat,
   238 void solveAndCheck(LpSolver& lp, LpSolver::ProblemType stat,
   254                    double exp_opt) {
   239                    double exp_opt) {
   255   using std::string;
   240   using std::string;
   256   lp.solve();
   241   lp.solve();
   257   //int decimal,sign;
   242 
   258   std::ostringstream buf;
   243   std::ostringstream buf;
   259   buf << "Primalstatus should be: " << int(stat);
   244   buf << "PrimalType should be: " << int(stat) << int(lp.primalType());
   260 
   245 
   261   //  itoa(stat,buf1, 10);
   246   check(lp.primalType()==stat, buf.str());
   262   check(lp.primalStatus()==stat, buf.str());
   247 
   263 
   248   if (stat ==  LpSolver::OPTIMAL) {
   264   if (stat ==  LpSolverBase::OPTIMAL) {
       
   265     std::ostringstream sbuf;
   249     std::ostringstream sbuf;
   266     sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
   250     sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
   267     check(std::abs(lp.primalValue()-exp_opt) < 1e-3, sbuf.str());
   251     check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str());
   268     //+ecvt(exp_opt,2)
       
   269   }
   252   }
   270 }
   253 }
   271 
   254 
   272 void aTest(LpSolverBase & lp)
   255 void aTest(LpSolver & lp)
   273 {
   256 {
   274   typedef LpSolverBase LP;
   257   typedef LpSolver LP;
   275 
   258 
   276  //The following example is very simple
   259  //The following example is very simple
   277 
   260 
   278   typedef LpSolverBase::Row Row;
   261   typedef LpSolver::Row Row;
   279   typedef LpSolverBase::Col Col;
   262   typedef LpSolver::Col Col;
   280 
   263 
   281 
   264 
   282   Col x1 = lp.addCol();
   265   Col x1 = lp.addCol();
   283   Col x2 = lp.addCol();
   266   Col x2 = lp.addCol();
   284 
   267 
   285 
   268 
   286   //Constraints
   269   //Constraints
   287   Row upright=lp.addRow(x1+x2 <=1);
   270   Row upright=lp.addRow(x1+2*x2 <=1);
   288   lp.addRow(x1+x2 >=-1);
   271   lp.addRow(x1+x2 >=-1);
   289   lp.addRow(x1-x2 <=1);
   272   lp.addRow(x1-x2 <=1);
   290   lp.addRow(x1-x2 >=-1);
   273   lp.addRow(x1-x2 >=-1);
   291   //Nonnegativity of the variables
   274   //Nonnegativity of the variables
   292   lp.colLowerBound(x1, 0);
   275   lp.colLowerBound(x1, 0);
   293   lp.colLowerBound(x2, 0);
   276   lp.colLowerBound(x2, 0);
   294   //Objective function
   277   //Objective function
   295   lp.obj(x1+x2);
   278   lp.obj(x1+x2);
   296 
   279 
   297   lp.max();
   280   lp.sense(lp.MAX);
   298 
   281 
   299   //Testing the problem retrieving routines
   282   //Testing the problem retrieving routines
   300   check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
   283   check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
   301   check(lp.isMax(),"This is a maximization!");
   284   check(lp.sense() == lp.MAX,"This is a maximization!");
   302   check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
   285   check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
   303   //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   286   check(lp.colLowerBound(x1)==0,
   304   check(  lp.colLowerBound(x1)==0,
   287         "The lower bound for variable x1 should be 0.");
   305           "The lower bound for variable x1 should be 0.");
   288   check(lp.colUpperBound(x1)==LpSolver::INF,
   306   check(  lp.colUpperBound(x1)==LpSolverBase::INF,
   289         "The upper bound for variable x1 should be infty.");
   307           "The upper bound for variable x1 should be infty.");
   290   check(lp.rowLowerBound(upright) == -LpSolver::INF,
   308   LpSolverBase::Value lb,ub;
   291         "The lower bound for the first row should be -infty.");
   309   lp.getRowBounds(upright,lb,ub);
   292   check(lp.rowUpperBound(upright)==1,
   310   check(  lb==-LpSolverBase::INF,
   293         "The upper bound for the first row should be 1.");
   311           "The lower bound for the first row should be -infty.");
   294   LpSolver::Expr e = lp.row(upright);
   312   check(  ub==1,"The upper bound for the first row should be 1.");
   295   check(e[x1] == 1, "The first coefficient should 1.");
   313   LpSolverBase::Expr e = lp.row(upright);
   296   check(e[x2] == 2, "The second coefficient should 1.");
   314   check(  e.size() == 2, "The row retrieval gives back wrong expression.");
   297 
   315   check(  e[x1] == 1, "The first coefficient should 1.");
   298   lp.row(upright, x1+x2 <=1);
   316   check(  e[x2] == 1, "The second coefficient should 1.");
   299   e = lp.row(upright);
   317 
   300   check(e[x1] == 1, "The first coefficient should 1.");
   318   LpSolverBase::DualExpr de = lp.col(x1);
   301   check(e[x2] == 1, "The second coefficient should 1.");
   319   check(  de.size() == 4, "The col retrieval gives back wrong expression.");
   302 
       
   303   LpSolver::DualExpr de = lp.col(x1);
   320   check(  de[upright] == 1, "The first coefficient should 1.");
   304   check(  de[upright] == 1, "The first coefficient should 1.");
   321 
   305 
   322   LpSolverBase* clp = lp.copyLp();
   306   LpSolver* clp = lp.cloneSolver();
   323 
   307 
   324   //Testing the problem retrieving routines
   308   //Testing the problem retrieving routines
   325   check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
   309   check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
   326   check(clp->isMax(),"This is a maximization!");
   310   check(clp->sense() == clp->MAX,"This is a maximization!");
   327   check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
   311   check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
   328   //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   312   //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   329   check(  clp->colLowerBound(x1)==0,
   313   check(clp->colLowerBound(x1)==0,
   330           "The lower bound for variable x1 should be 0.");
   314         "The lower bound for variable x1 should be 0.");
   331   check(  clp->colUpperBound(x1)==LpSolverBase::INF,
   315   check(clp->colUpperBound(x1)==LpSolver::INF,
   332           "The upper bound for variable x1 should be infty.");
   316         "The upper bound for variable x1 should be infty.");
   333 
   317 
   334   clp->getRowBounds(upright,lb,ub);
   318   check(lp.rowLowerBound(upright)==-LpSolver::INF,
   335   check(  lb==-LpSolverBase::INF,
   319         "The lower bound for the first row should be -infty.");
   336           "The lower bound for the first row should be -infty.");
   320   check(lp.rowUpperBound(upright)==1,
   337   check(  ub==1,"The upper bound for the first row should be 1.");
   321         "The upper bound for the first row should be 1.");
   338   e = clp->row(upright);
   322   e = clp->row(upright);
   339   check(  e.size() == 2, "The row retrieval gives back wrong expression.");
   323   check(e[x1] == 1, "The first coefficient should 1.");
   340   check(  e[x1] == 1, "The first coefficient should 1.");
   324   check(e[x2] == 1, "The second coefficient should 1.");
   341   check(  e[x2] == 1, "The second coefficient should 1.");
       
   342 
   325 
   343   de = clp->col(x1);
   326   de = clp->col(x1);
   344   check(  de.size() == 4, "The col retrieval gives back wrong expression.");
   327   check(de[upright] == 1, "The first coefficient should 1.");
   345   check(  de[upright] == 1, "The first coefficient should 1.");
       
   346 
   328 
   347   delete clp;
   329   delete clp;
   348 
   330 
   349   //Maximization of x1+x2
   331   //Maximization of x1+x2
   350   //over the triangle with vertices (0,0) (0,1) (1,0)
   332   //over the triangle with vertices (0,0) (0,1) (1,0)
   351   double expected_opt=1;
   333   double expected_opt=1;
   352   solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
   334   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   353 
   335 
   354   //Minimization
   336   //Minimization
   355   lp.min();
   337   lp.sense(lp.MIN);
   356   expected_opt=0;
   338   expected_opt=0;
   357   solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
   339   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   358 
   340 
   359   //Vertex (-1,0) instead of (0,0)
   341   //Vertex (-1,0) instead of (0,0)
   360   lp.colLowerBound(x1, -LpSolverBase::INF);
   342   lp.colLowerBound(x1, -LpSolver::INF);
   361   expected_opt=-1;
   343   expected_opt=-1;
   362   solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
   344   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   363 
   345 
   364   //Erase one constraint and return to maximization
   346   //Erase one constraint and return to maximization
   365   lp.eraseRow(upright);
   347   lp.erase(upright);
   366   lp.max();
   348   lp.sense(lp.MAX);
   367   expected_opt=LpSolverBase::INF;
   349   expected_opt=LpSolver::INF;
   368   solveAndCheck(lp, LpSolverBase::INFINITE, expected_opt);
   350   solveAndCheck(lp, LpSolver::UNBOUNDED, expected_opt);
   369 
   351 
   370   //Infeasibilty
   352   //Infeasibilty
   371   lp.addRow(x1+x2 <=-2);
   353   lp.addRow(x1+x2 <=-2);
   372   solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt);
   354   solveAndCheck(lp, LpSolver::INFEASIBLE, expected_opt);
   373 
       
   374   //Change problem and forget to solve
       
   375   lp.min();
       
   376   check(lp.primalStatus()==LpSolverBase::UNDEFINED,
       
   377         "Primalstatus should be UNDEFINED");
       
   378 
       
   379 
       
   380 //   lp.solve();
       
   381 //   if (lp.primalStatus()==LpSolverBase::OPTIMAL){
       
   382 //     std::cout<< "Z = "<<lp.primalValue()
       
   383 //              << " (error = " << lp.primalValue()-expected_opt
       
   384 //              << "); x1 = "<<lp.primal(x1)
       
   385 //              << "; x2 = "<<lp.primal(x2)
       
   386 //              <<std::endl;
       
   387 
       
   388 //   }
       
   389 //   else{
       
   390 //     std::cout<<lp.primalStatus()<<std::endl;
       
   391 //     std::cout<<"Optimal solution not found!"<<std::endl;
       
   392 //   }
       
   393 
       
   394 
       
   395 
   355 
   396 }
   356 }
   397 
       
   398 
   357 
   399 int main()
   358 int main()
   400 {
   359 {
   401   LpSkeleton lp_skel;
   360   LpSkeleton lp_skel;
   402   lpTest(lp_skel);
   361   lpTest(lp_skel);
   403 
   362 
   404 #ifdef HAVE_GLPK
   363 #ifdef HAVE_GLPK
   405   LpGlpk lp_glpk1,lp_glpk2;
   364   {
   406   lpTest(lp_glpk1);
   365     LpGlpk lp_glpk1,lp_glpk2;
   407   aTest(lp_glpk2);
   366     lpTest(lp_glpk1);
       
   367     aTest(lp_glpk2);
       
   368   }
   408 #endif
   369 #endif
   409 
   370 
   410 #ifdef HAVE_CPLEX
   371 #ifdef HAVE_CPLEX
   411   LpCplex lp_cplex1,lp_cplex2;
   372   try {
   412   lpTest(lp_cplex1);
   373     LpCplex lp_cplex1,lp_cplex2;
   413   aTest(lp_cplex2);
   374     lpTest(lp_cplex1);
       
   375     aTest(lp_cplex2);
       
   376   } catch (CplexEnv::LicenseError& error) {
       
   377 #ifdef LEMON_FORCE_CPLEX_CHECK
       
   378     check(false, error.what());
       
   379 #else
       
   380     std::cerr << error.what() << std::endl;
       
   381     std::cerr << "Cplex license check failed, lp check skipped" << std::endl;
       
   382 #endif
       
   383   }
   414 #endif
   384 #endif
   415 
   385 
   416 #ifdef HAVE_SOPLEX
   386 #ifdef HAVE_SOPLEX
   417   LpSoplex lp_soplex1,lp_soplex2;
   387   {
   418   lpTest(lp_soplex1);
   388     LpSoplex lp_soplex1,lp_soplex2;
   419   aTest(lp_soplex2);
   389     lpTest(lp_soplex1);
       
   390     aTest(lp_soplex2);
       
   391   }
       
   392 #endif
       
   393 
       
   394 #ifdef HAVE_CLP
       
   395   {
       
   396     LpClp lp_clp1,lp_clp2;
       
   397     lpTest(lp_clp1);
       
   398     aTest(lp_clp2);
       
   399   }
   420 #endif
   400 #endif
   421 
   401 
   422   return 0;
   402   return 0;
   423 }
   403 }