test/lp_test.cc
author Balazs Dezso <deba@inf.elte.hu>
Mon, 23 Feb 2009 23:44:29 +0100
changeset 538 ba124394367a
parent 537 0fec6a017ead
child 542 fc6c7aab4b8d
permissions -rw-r--r--
Automatic GLPK env deallocation (#213)
     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-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 
    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 
   201   }
   202 
   203   {
   204     LP::DualExpr e,f,g;
   205     LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
   206       p4 = INVALID, p5 = INVALID;
   207 
   208     e[p1]=2;
   209     e[p1]+=2;
   210     e[p1]-=2;
   211 
   212     e=p1;
   213     e=f;
   214 
   215     e+=p1;
   216     e+=f;
   217 
   218     e-=p1;
   219     e-=f;
   220 
   221     e*=2;
   222     e*=2.2;
   223     e/=2;
   224     e/=2.2;
   225 
   226     e=((p1+p2)+(p1-p2)+
   227        (p1+f)+(f+p1)+(f+g)+
   228        (p1-f)+(f-p1)+(f-g)+
   229        2.2*f+f*2.2+f/2.2+
   230        2*f+f*2+f/2+
   231        2.2*p1+p1*2.2+p1/2.2+
   232        2*p1+p1*2+p1/2
   233        );
   234   }
   235 
   236 }
   237 
   238 void solveAndCheck(LpSolver& lp, LpSolver::ProblemType stat,
   239                    double exp_opt) {
   240   using std::string;
   241   lp.solve();
   242 
   243   std::ostringstream buf;
   244   buf << "PrimalType should be: " << int(stat) << int(lp.primalType());
   245 
   246   check(lp.primalType()==stat, buf.str());
   247 
   248   if (stat ==  LpSolver::OPTIMAL) {
   249     std::ostringstream sbuf;
   250     sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
   251     check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str());
   252   }
   253 }
   254 
   255 void aTest(LpSolver & lp)
   256 {
   257   typedef LpSolver LP;
   258 
   259  //The following example is very simple
   260 
   261   typedef LpSolver::Row Row;
   262   typedef LpSolver::Col Col;
   263 
   264 
   265   Col x1 = lp.addCol();
   266   Col x2 = lp.addCol();
   267 
   268 
   269   //Constraints
   270   Row upright=lp.addRow(x1+2*x2 <=1);
   271   lp.addRow(x1+x2 >=-1);
   272   lp.addRow(x1-x2 <=1);
   273   lp.addRow(x1-x2 >=-1);
   274   //Nonnegativity of the variables
   275   lp.colLowerBound(x1, 0);
   276   lp.colLowerBound(x2, 0);
   277   //Objective function
   278   lp.obj(x1+x2);
   279 
   280   lp.sense(lp.MAX);
   281 
   282   //Testing the problem retrieving routines
   283   check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
   284   check(lp.sense() == lp.MAX,"This is a maximization!");
   285   check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
   286   check(lp.colLowerBound(x1)==0,
   287         "The lower bound for variable x1 should be 0.");
   288   check(lp.colUpperBound(x1)==LpSolver::INF,
   289         "The upper bound for variable x1 should be infty.");
   290   check(lp.rowLowerBound(upright) == -LpSolver::INF,
   291         "The lower bound for the first row should be -infty.");
   292   check(lp.rowUpperBound(upright)==1,
   293         "The upper bound for the first row should be 1.");
   294   LpSolver::Expr e = lp.row(upright);
   295   check(e[x1] == 1, "The first coefficient should 1.");
   296   check(e[x2] == 2, "The second coefficient should 1.");
   297 
   298   lp.row(upright, x1+x2 <=1);
   299   e = lp.row(upright);
   300   check(e[x1] == 1, "The first coefficient should 1.");
   301   check(e[x2] == 1, "The second coefficient should 1.");
   302 
   303   LpSolver::DualExpr de = lp.col(x1);
   304   check(  de[upright] == 1, "The first coefficient should 1.");
   305 
   306   LpSolver* clp = lp.cloneSolver();
   307 
   308   //Testing the problem retrieving routines
   309   check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
   310   check(clp->sense() == clp->MAX,"This is a maximization!");
   311   check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
   312   //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   313   check(clp->colLowerBound(x1)==0,
   314         "The lower bound for variable x1 should be 0.");
   315   check(clp->colUpperBound(x1)==LpSolver::INF,
   316         "The upper bound for variable x1 should be infty.");
   317 
   318   check(lp.rowLowerBound(upright)==-LpSolver::INF,
   319         "The lower bound for the first row should be -infty.");
   320   check(lp.rowUpperBound(upright)==1,
   321         "The upper bound for the first row should be 1.");
   322   e = clp->row(upright);
   323   check(e[x1] == 1, "The first coefficient should 1.");
   324   check(e[x2] == 1, "The second coefficient should 1.");
   325 
   326   de = clp->col(x1);
   327   check(de[upright] == 1, "The first coefficient should 1.");
   328 
   329   delete clp;
   330 
   331   //Maximization of x1+x2
   332   //over the triangle with vertices (0,0) (0,1) (1,0)
   333   double expected_opt=1;
   334   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   335 
   336   //Minimization
   337   lp.sense(lp.MIN);
   338   expected_opt=0;
   339   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   340 
   341   //Vertex (-1,0) instead of (0,0)
   342   lp.colLowerBound(x1, -LpSolver::INF);
   343   expected_opt=-1;
   344   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   345 
   346   //Erase one constraint and return to maximization
   347   lp.erase(upright);
   348   lp.sense(lp.MAX);
   349   expected_opt=LpSolver::INF;
   350   solveAndCheck(lp, LpSolver::UNBOUNDED, expected_opt);
   351 
   352   //Infeasibilty
   353   lp.addRow(x1+x2 <=-2);
   354   solveAndCheck(lp, LpSolver::INFEASIBLE, expected_opt);
   355 
   356 }
   357 
   358 int main()
   359 {
   360   LpSkeleton lp_skel;
   361   lpTest(lp_skel);
   362 
   363 #ifdef HAVE_GLPK
   364   {
   365     GlpkLp lp_glpk1,lp_glpk2;
   366     lpTest(lp_glpk1);
   367     aTest(lp_glpk2);
   368   }
   369 #endif
   370 
   371 #ifdef HAVE_CPLEX
   372   try {
   373     CplexLp lp_cplex1,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   }
   384 #endif
   385 
   386 #ifdef HAVE_SOPLEX
   387   {
   388     SoplexLp lp_soplex1,lp_soplex2;
   389     lpTest(lp_soplex1);
   390     aTest(lp_soplex2);
   391   }
   392 #endif
   393 
   394 #ifdef HAVE_CLP
   395   {
   396     ClpLp lp_clp1,lp_clp2;
   397     lpTest(lp_clp1);
   398     aTest(lp_clp2);
   399   }
   400 #endif
   401 
   402   return 0;
   403 }