test/lp_test.cc
author deba
Fri, 16 May 2008 12:28:23 +0000
changeset 2609 c36f00f19f2b
parent 2571 f5c6e63f51d1
child 2637 bafe730864da
permissions -rw-r--r--
Fixing lp constraints (solving ticket #92)
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2008
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #include <sstream>
    20 #include <lemon/lp_skeleton.h>
    21 #include "test_tools.h"
    22 #include <lemon/tolerance.h>
    23 #include <lemon/lp_utils.h>
    24 
    25 #ifdef HAVE_CONFIG_H
    26 #include <lemon/config.h>
    27 #endif
    28 
    29 #ifdef HAVE_GLPK
    30 #include <lemon/lp_glpk.h>
    31 #endif
    32 
    33 #ifdef HAVE_CPLEX
    34 #include <lemon/lp_cplex.h>
    35 #endif
    36 
    37 #ifdef HAVE_SOPLEX
    38 #include <lemon/lp_soplex.h>
    39 #endif
    40 
    41 using namespace lemon;
    42 
    43 void lpTest(LpSolverBase & lp)
    44 {
    45 
    46 
    47 
    48   typedef LpSolverBase LP;
    49 
    50   std::vector<LP::Col> x(10);
    51   //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
    52   lp.addColSet(x);
    53   lp.colLowerBound(x,1);
    54   lp.colUpperBound(x,1);
    55   lp.colBounds(x,1,2);
    56 #ifndef GYORSITAS
    57 
    58   std::vector<LP::Col> y(10);
    59   lp.addColSet(y);
    60 
    61   lp.colLowerBound(y,1);
    62   lp.colUpperBound(y,1);
    63   lp.colBounds(y,1,2);
    64 
    65   std::map<int,LP::Col> z;
    66   
    67   z.insert(std::make_pair(12,INVALID));
    68   z.insert(std::make_pair(2,INVALID));
    69   z.insert(std::make_pair(7,INVALID));
    70   z.insert(std::make_pair(5,INVALID));
    71 
    72   lp.addColSet(z);
    73 
    74   lp.colLowerBound(z,1);
    75   lp.colUpperBound(z,1);
    76   lp.colBounds(z,1,2);
    77 
    78   {
    79     LP::Expr e,f,g;
    80     LP::Col p1,p2,p3,p4,p5;
    81     LP::Constr c;
    82     
    83     p1=lp.addCol();
    84     p2=lp.addCol();
    85     p3=lp.addCol();
    86     p4=lp.addCol();
    87     p5=lp.addCol();
    88     
    89     e[p1]=2;
    90     e.constComp()=12;
    91     e[p1]+=2;
    92     e.constComp()+=12;
    93     e[p1]-=2;
    94     e.constComp()-=12;
    95     
    96     e=2;
    97     e=2.2;
    98     e=p1;
    99     e=f;
   100     
   101     e+=2;
   102     e+=2.2;
   103     e+=p1;
   104     e+=f;
   105     
   106     e-=2;
   107     e-=2.2;
   108     e-=p1;
   109     e-=f;
   110     
   111     e*=2;
   112     e*=2.2;
   113     e/=2;
   114     e/=2.2;
   115     
   116     e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
   117        (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
   118        (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
   119        2.2*f+f*2.2+f/2.2+
   120        2*f+f*2+f/2+
   121        2.2*p1+p1*2.2+p1/2.2+
   122        2*p1+p1*2+p1/2
   123        );
   124 
   125 
   126     c = (e  <= f  );
   127     c = (e  <= 2.2);
   128     c = (e  <= 2  );
   129     c = (e  <= p1 );
   130     c = (2.2<= f  );
   131     c = (2  <= f  );
   132     c = (p1 <= f  );
   133     c = (p1 <= p2 );
   134     c = (p1 <= 2.2);
   135     c = (p1 <= 2  );
   136     c = (2.2<= p2 );
   137     c = (2  <= p2 );
   138     
   139     c = (e  >= f  );
   140     c = (e  >= 2.2);
   141     c = (e  >= 2  );
   142     c = (e  >= p1 );
   143     c = (2.2>= f  );
   144     c = (2  >= f  );
   145     c = (p1 >= f  );
   146     c = (p1 >= p2 );
   147     c = (p1 >= 2.2);
   148     c = (p1 >= 2  );
   149     c = (2.2>= p2 );
   150     c = (2  >= p2 );
   151     
   152     c = (e  == f  );
   153     c = (e  == 2.2);
   154     c = (e  == 2  );
   155     c = (e  == p1 );
   156     c = (2.2== f  );
   157     c = (2  == f  );
   158     c = (p1 == f  );
   159     //c = (p1 == p2 );
   160     c = (p1 == 2.2);
   161     c = (p1 == 2  );
   162     c = (2.2== p2 );
   163     c = (2  == p2 );
   164     
   165     c = (2 <= e <= 3);
   166     c = (2 <= p1<= 3);
   167     
   168     c = (2 >= e >= 3);
   169     c = (2 >= p1>= 3);
   170     
   171     e[x[3]]=2;
   172     e[x[3]]=4;
   173     e[x[3]]=1;
   174     e.constComp()=12;
   175     
   176     lp.addRow(LP::INF,e,23);
   177     lp.addRow(LP::INF,3.0*(x[1]+x[2]/2)-x[3],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);
   179     
   180     lp.addRow(x[1]+x[3]<=x[5]-3);
   181     lp.addRow(-7<=x[1]+x[3]-12<=3);
   182     lp.addRow(x[1]<=x[5]);
   183 
   184     std::ostringstream buf;
   185 
   186 
   187     //Checking the simplify function
   188 
   189 //     //How to check the simplify function? A map gives no information
   190 //     //on the question whether a given key is or is not stored in it, or
   191 //     //it does?
   192 //   Yes, it does, using the find() function.
   193     e=((p1+p2)+(p1-p2));
   194     e.simplify();
   195     buf << "Coeff. of p2 should be 0";
   196     //    std::cout<<e[p1]<<e[p2]<<e[p3]<<std::endl;
   197     check(e.find(p2)==e.end(), buf.str());
   198     
   199      
   200 
   201 
   202     e=((p1+p2)+(p1-0.99*p2));
   203     //e.prettyPrint(std::cout);
   204     //(e<=2).prettyPrint(std::cout);
   205     double tolerance=0.001;
   206     e.simplify(tolerance);
   207     buf << "Coeff. of p2 should be 0.01";
   208     check(e[p2]>0, buf.str());
   209     
   210     tolerance=0.02;
   211     e.simplify(tolerance);
   212     buf << "Coeff. of p2 should be 0";
   213     check(e.find(p2)==e.end(), buf.str());
   214     
   215 
   216   }
   217   
   218   {
   219     LP::DualExpr e,f,g;
   220     LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID, 
   221       p4 = INVALID, p5 = INVALID;
   222     
   223     e[p1]=2;
   224     e[p1]+=2;
   225     e[p1]-=2;
   226     
   227     e=p1;
   228     e=f;
   229     
   230     e+=p1;
   231     e+=f;
   232     
   233     e-=p1;
   234     e-=f;
   235     
   236     e*=2;
   237     e*=2.2;
   238     e/=2;
   239     e/=2.2;
   240     
   241     e=((p1+p2)+(p1-p2)+
   242        (p1+f)+(f+p1)+(f+g)+
   243        (p1-f)+(f-p1)+(f-g)+
   244        2.2*f+f*2.2+f/2.2+
   245        2*f+f*2+f/2+
   246        2.2*p1+p1*2.2+p1/2.2+
   247        2*p1+p1*2+p1/2
   248        );
   249   }
   250   
   251 #endif
   252 }
   253 
   254 void solveAndCheck(LpSolverBase& lp, LpSolverBase::SolutionStatus stat, 
   255 		   double exp_opt) {
   256   using std::string;
   257   lp.solve();
   258   //int decimal,sign;
   259   std::ostringstream buf;
   260   buf << "Primalstatus should be: " << int(stat);
   261 
   262   //  itoa(stat,buf1, 10);
   263   check(lp.primalStatus()==stat, buf.str());
   264 
   265   if (stat ==  LpSolverBase::OPTIMAL) {
   266     std::ostringstream sbuf;
   267     sbuf << "Wrong optimal value: the right optimum is " << exp_opt; 
   268     check(std::abs(lp.primalValue()-exp_opt) < 1e-3, sbuf.str());
   269     //+ecvt(exp_opt,2)
   270   }
   271 }
   272  
   273 void aTest(LpSolverBase & lp)
   274 {
   275   typedef LpSolverBase LP;
   276 
   277  //The following example is very simple
   278 
   279   typedef LpSolverBase::Row Row;
   280   typedef LpSolverBase::Col Col;
   281 
   282 
   283   Col x1 = lp.addCol();
   284   Col x2 = lp.addCol();
   285 
   286 
   287   //Constraints
   288   Row upright=lp.addRow(x1+x2 <=1);  
   289   lp.addRow(x1+x2 >=-1);  
   290   lp.addRow(x1-x2 <=1);  
   291   lp.addRow(x1-x2 >=-1);  
   292   //Nonnegativity of the variables
   293   lp.colLowerBound(x1, 0);
   294   lp.colLowerBound(x2, 0);
   295   //Objective function
   296   lp.obj(x1+x2);
   297 
   298   lp.max();
   299 
   300   //Testing the problem retrieving routines
   301   check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
   302   check(lp.isMax(),"This is a maximization!");
   303   check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
   304   //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   305   check(  lp.colLowerBound(x1)==0,"The lower bound for variable x1 should be 0.");
   306   check(  lp.colUpperBound(x1)==LpSolverBase::INF,"The upper bound for variable x1 should be infty.");
   307   LpSolverBase::Value lb,ub;
   308   lp.getRowBounds(upright,lb,ub);
   309   check(  lb==-LpSolverBase::INF,"The lower bound for the first row should be -infty.");
   310   check(  ub==1,"The upper bound for the first row should be 1.");
   311   LpSolverBase::Expr e = lp.row(upright);
   312   check(  e.size() == 2, "The row retrieval gives back wrong expression.");
   313   check(  e[x1] == 1, "The first coefficient should 1.");
   314   check(  e[x2] == 1, "The second coefficient should 1.");
   315 
   316   LpSolverBase::DualExpr de = lp.col(x1);
   317   check(  de.size() == 4, "The col retrieval gives back wrong expression.");
   318   check(  de[upright] == 1, "The first coefficient should 1.");
   319 
   320   LpSolverBase* clp = lp.copyLp();
   321 
   322   //Testing the problem retrieving routines
   323   check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
   324   check(clp->isMax(),"This is a maximization!");
   325   check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
   326   //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   327   check(  clp->colLowerBound(x1)==0,"The lower bound for variable x1 should be 0.");
   328   std::cerr << clp->colUpperBound(x1) << std::endl;
   329   check(  clp->colUpperBound(x1)==LpSolverBase::INF,"The upper bound for variable x1 should be infty.");
   330 
   331   clp->getRowBounds(upright,lb,ub);
   332   check(  lb==-LpSolverBase::INF,"The lower bound for the first row should be -infty.");
   333   check(  ub==1,"The upper bound for the first row should be 1.");
   334   e = clp->row(upright);
   335   check(  e.size() == 2, "The row retrieval gives back wrong expression.");
   336   check(  e[x1] == 1, "The first coefficient should 1.");
   337   check(  e[x2] == 1, "The second coefficient should 1.");
   338 
   339   de = clp->col(x1);
   340   check(  de.size() == 4, "The col retrieval gives back wrong expression.");
   341   check(  de[upright] == 1, "The first coefficient should 1.");
   342   
   343   delete clp;
   344 
   345   //Maximization of x1+x2
   346   //over the triangle with vertices (0,0) (0,1) (1,0)
   347   double expected_opt=1;
   348   solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
   349   
   350   //Minimization
   351   lp.min();
   352   expected_opt=0;
   353   solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
   354   
   355   //Vertex (-1,0) instead of (0,0)
   356   lp.colLowerBound(x1, -LpSolverBase::INF);
   357   expected_opt=-1;
   358   solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
   359 
   360   //Erase one constraint and return to maximization
   361   lp.eraseRow(upright);
   362   lp.max();
   363   expected_opt=LpSolverBase::INF;
   364   solveAndCheck(lp, LpSolverBase::INFINITE, expected_opt);
   365 
   366   //Infeasibilty
   367   lp.addRow(x1+x2 <=-2);  
   368   solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt);
   369 
   370   //Change problem and forget to solve
   371   lp.min();
   372   check(lp.primalStatus()==LpSolverBase::UNDEFINED,"Primalstatus should be UNDEFINED");
   373 
   374 
   375 //   lp.solve();
   376 //   if (lp.primalStatus()==LpSolverBase::OPTIMAL){
   377 //     std::cout<< "Z = "<<lp.primalValue()
   378 // 	     << " (error = " << lp.primalValue()-expected_opt
   379 // 	     << "); x1 = "<<lp.primal(x1)
   380 // 	     << "; x2 = "<<lp.primal(x2)
   381 // 	     <<std::endl;
   382     
   383 //   }
   384 //   else{
   385 //     std::cout<<lp.primalStatus()<<std::endl;
   386 //     std::cout<<"Optimal solution not found!"<<std::endl;
   387 //   }
   388 
   389  
   390 
   391 }
   392 
   393 
   394 int main() 
   395 {
   396   LpSkeleton lp_skel;
   397   lpTest(lp_skel);
   398 
   399 #ifdef HAVE_GLPK
   400   LpGlpk lp_glpk1,lp_glpk2;
   401   lpTest(lp_glpk1);
   402   aTest(lp_glpk2);
   403 #endif
   404 
   405 #ifdef HAVE_CPLEX
   406   LpCplex lp_cplex1,lp_cplex2;
   407   lpTest(lp_cplex1);
   408   aTest(lp_cplex2);
   409 #endif
   410 
   411 #ifdef HAVE_SOPLEX
   412   LpSoplex lp_soplex1,lp_soplex2;
   413   lpTest(lp_soplex1);
   414   aTest(lp_soplex2);
   415 #endif
   416 
   417   return 0;
   418 }