test/lp_test.cc
author Alpar Juttner <alpar@cs.elte.hu>
Tue, 21 Apr 2009 15:18:54 +0100
changeset 611 85cb3aa71cce
parent 551 9d0d7e20f76d
child 627 20dac2104519
permissions -rw-r--r--
Merge and fix
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     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 
    24 #ifdef HAVE_CONFIG_H
    25 #include <lemon/config.h>
    26 #endif
    27 
    28 #ifdef HAVE_GLPK
    29 #include <lemon/glpk.h>
    30 #endif
    31 
    32 #ifdef HAVE_CPLEX
    33 #include <lemon/cplex.h>
    34 #endif
    35 
    36 #ifdef HAVE_SOPLEX
    37 #include <lemon/soplex.h>
    38 #endif
    39 
    40 #ifdef HAVE_CLP
    41 #include <lemon/clp.h>
    42 #endif
    43 
    44 using namespace lemon;
    45 
    46 void lpTest(LpSolver& lp)
    47 {
    48 
    49   typedef LpSolver LP;
    50 
    51   std::vector<LP::Col> x(10);
    52   //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
    53   lp.addColSet(x);
    54   lp.colLowerBound(x,1);
    55   lp.colUpperBound(x,1);
    56   lp.colBounds(x,1,2);
    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=12;
    91     e[p1]+=2;
    92     *e+=12;
    93     e[p1]-=2;
    94     *e-=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=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     e=((p1+p2)+(p1-0.99*p2));
   188     //e.prettyPrint(std::cout);
   189     //(e<=2).prettyPrint(std::cout);
   190     double tolerance=0.001;
   191     e.simplify(tolerance);
   192     buf << "Coeff. of p2 should be 0.01";
   193     check(e[p2]>0, buf.str());
   194 
   195     tolerance=0.02;
   196     e.simplify(tolerance);
   197     buf << "Coeff. of p2 should be 0";
   198     check(const_cast<const LpSolver::Expr&>(e)[p2]==0, buf.str());
   199 
   200     //Test for clone/new
   201     LP* lpnew = lp.newSolver();
   202     LP* lpclone = lp.cloneSolver();
   203     delete lpnew;
   204     delete lpclone;
   205 
   206   }
   207 
   208   {
   209     LP::DualExpr e,f,g;
   210     LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
   211       p4 = INVALID, p5 = INVALID;
   212 
   213     e[p1]=2;
   214     e[p1]+=2;
   215     e[p1]-=2;
   216 
   217     e=p1;
   218     e=f;
   219 
   220     e+=p1;
   221     e+=f;
   222 
   223     e-=p1;
   224     e-=f;
   225 
   226     e*=2;
   227     e*=2.2;
   228     e/=2;
   229     e/=2.2;
   230 
   231     e=((p1+p2)+(p1-p2)+
   232        (p1+f)+(f+p1)+(f+g)+
   233        (p1-f)+(f-p1)+(f-g)+
   234        2.2*f+f*2.2+f/2.2+
   235        2*f+f*2+f/2+
   236        2.2*p1+p1*2.2+p1/2.2+
   237        2*p1+p1*2+p1/2
   238        );
   239   }
   240 
   241 }
   242 
   243 void solveAndCheck(LpSolver& lp, LpSolver::ProblemType stat,
   244                    double exp_opt) {
   245   using std::string;
   246   lp.solve();
   247 
   248   std::ostringstream buf;
   249   buf << "PrimalType should be: " << int(stat) << int(lp.primalType());
   250 
   251   check(lp.primalType()==stat, buf.str());
   252 
   253   if (stat ==  LpSolver::OPTIMAL) {
   254     std::ostringstream sbuf;
   255     sbuf << "Wrong optimal value (" << lp.primal() <<") with "
   256          << lp.solverName() <<"\n     the right optimum is " << exp_opt;
   257     check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str());
   258   }
   259 }
   260 
   261 void aTest(LpSolver & lp)
   262 {
   263   typedef LpSolver LP;
   264 
   265  //The following example is very simple
   266 
   267   typedef LpSolver::Row Row;
   268   typedef LpSolver::Col Col;
   269 
   270 
   271   Col x1 = lp.addCol();
   272   Col x2 = lp.addCol();
   273 
   274 
   275   //Constraints
   276   Row upright=lp.addRow(x1+2*x2 <=1);
   277   lp.addRow(x1+x2 >=-1);
   278   lp.addRow(x1-x2 <=1);
   279   lp.addRow(x1-x2 >=-1);
   280   //Nonnegativity of the variables
   281   lp.colLowerBound(x1, 0);
   282   lp.colLowerBound(x2, 0);
   283   //Objective function
   284   lp.obj(x1+x2);
   285 
   286   lp.sense(lp.MAX);
   287 
   288   //Testing the problem retrieving routines
   289   check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
   290   check(lp.sense() == lp.MAX,"This is a maximization!");
   291   check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
   292   check(lp.colLowerBound(x1)==0,
   293         "The lower bound for variable x1 should be 0.");
   294   check(lp.colUpperBound(x1)==LpSolver::INF,
   295         "The upper bound for variable x1 should be infty.");
   296   check(lp.rowLowerBound(upright) == -LpSolver::INF,
   297         "The lower bound for the first row should be -infty.");
   298   check(lp.rowUpperBound(upright)==1,
   299         "The upper bound for the first row should be 1.");
   300   LpSolver::Expr e = lp.row(upright);
   301   check(e[x1] == 1, "The first coefficient should 1.");
   302   check(e[x2] == 2, "The second coefficient should 1.");
   303 
   304   lp.row(upright, x1+x2 <=1);
   305   e = lp.row(upright);
   306   check(e[x1] == 1, "The first coefficient should 1.");
   307   check(e[x2] == 1, "The second coefficient should 1.");
   308 
   309   LpSolver::DualExpr de = lp.col(x1);
   310   check(  de[upright] == 1, "The first coefficient should 1.");
   311 
   312   LpSolver* clp = lp.cloneSolver();
   313 
   314   //Testing the problem retrieving routines
   315   check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
   316   check(clp->sense() == clp->MAX,"This is a maximization!");
   317   check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
   318   //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   319   check(clp->colLowerBound(x1)==0,
   320         "The lower bound for variable x1 should be 0.");
   321   check(clp->colUpperBound(x1)==LpSolver::INF,
   322         "The upper bound for variable x1 should be infty.");
   323 
   324   check(lp.rowLowerBound(upright)==-LpSolver::INF,
   325         "The lower bound for the first row should be -infty.");
   326   check(lp.rowUpperBound(upright)==1,
   327         "The upper bound for the first row should be 1.");
   328   e = clp->row(upright);
   329   check(e[x1] == 1, "The first coefficient should 1.");
   330   check(e[x2] == 1, "The second coefficient should 1.");
   331 
   332   de = clp->col(x1);
   333   check(de[upright] == 1, "The first coefficient should 1.");
   334 
   335   delete clp;
   336 
   337   //Maximization of x1+x2
   338   //over the triangle with vertices (0,0) (0,1) (1,0)
   339   double expected_opt=1;
   340   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   341 
   342   //Minimization
   343   lp.sense(lp.MIN);
   344   expected_opt=0;
   345   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   346 
   347   //Vertex (-1,0) instead of (0,0)
   348   lp.colLowerBound(x1, -LpSolver::INF);
   349   expected_opt=-1;
   350   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   351 
   352   //Erase one constraint and return to maximization
   353   lp.erase(upright);
   354   lp.sense(lp.MAX);
   355   expected_opt=LpSolver::INF;
   356   solveAndCheck(lp, LpSolver::UNBOUNDED, expected_opt);
   357 
   358   //Infeasibilty
   359   lp.addRow(x1+x2 <=-2);
   360   solveAndCheck(lp, LpSolver::INFEASIBLE, expected_opt);
   361 
   362 }
   363 
   364 template<class LP>
   365 void cloneTest()
   366 {
   367   //Test for clone/new
   368 
   369   LP* lp = new LP();
   370   LP* lpnew = lp->newSolver();
   371   LP* lpclone = lp->cloneSolver();
   372   delete lp;
   373   delete lpnew;
   374   delete lpclone;
   375 }
   376 
   377 int main()
   378 {
   379   LpSkeleton lp_skel;
   380   lpTest(lp_skel);
   381 
   382 #ifdef HAVE_GLPK
   383   {
   384     GlpkLp lp_glpk1,lp_glpk2;
   385     lpTest(lp_glpk1);
   386     aTest(lp_glpk2);
   387     cloneTest<GlpkLp>();
   388   }
   389 #endif
   390 
   391 #ifdef HAVE_CPLEX
   392   try {
   393     CplexLp lp_cplex1,lp_cplex2;
   394     lpTest(lp_cplex1);
   395     aTest(lp_cplex2);
   396     cloneTest<CplexLp>();
   397   } catch (CplexEnv::LicenseError& error) {
   398     check(false, error.what());
   399   }
   400 #endif
   401 
   402 #ifdef HAVE_SOPLEX
   403   {
   404     SoplexLp lp_soplex1,lp_soplex2;
   405     lpTest(lp_soplex1);
   406     aTest(lp_soplex2);
   407     cloneTest<SoplexLp>();
   408   }
   409 #endif
   410 
   411 #ifdef HAVE_CLP
   412   {
   413     ClpLp lp_clp1,lp_clp2;
   414     lpTest(lp_clp1);
   415     aTest(lp_clp2);
   416     cloneTest<ClpLp>();
   417   }
   418 #endif
   419 
   420   return 0;
   421 }