test/lp_test.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Thu, 17 Oct 2013 09:30:57 +0200
changeset 1297 c0c2f5c87aa6
parent 1092 2eebc8f7dca5
child 1141 c3a7ca108705
child 1232 fc3854d936f7
permissions -rw-r--r--
Rename field in min cost flow codes (#478)
     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 #include <lemon/config.h>
    25 
    26 #ifdef LEMON_HAVE_GLPK
    27 #include <lemon/glpk.h>
    28 #endif
    29 
    30 #ifdef LEMON_HAVE_CPLEX
    31 #include <lemon/cplex.h>
    32 #endif
    33 
    34 #ifdef LEMON_HAVE_SOPLEX
    35 #include <lemon/soplex.h>
    36 #endif
    37 
    38 #ifdef LEMON_HAVE_CLP
    39 #include <lemon/clp.h>
    40 #endif
    41 
    42 using namespace lemon;
    43 
    44 int countCols(LpBase & lp) {
    45   int count=0;
    46   for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
    47   return count;
    48 }
    49 
    50 int countRows(LpBase & lp) {
    51   int count=0;
    52   for (LpBase::RowIt r(lp); r!=INVALID; ++r) ++count;
    53   return count;
    54 }
    55 
    56 
    57 void lpTest(LpSolver& lp)
    58 {
    59 
    60   typedef LpSolver LP;
    61 
    62   // Test LpBase::clear()
    63   check(countRows(lp)==0, "Wrong number of rows");
    64   check(countCols(lp)==0, "Wrong number of cols");
    65   lp.addCol(); lp.addRow(); lp.addRow();
    66   check(countRows(lp)==2, "Wrong number of rows");
    67   check(countCols(lp)==1, "Wrong number of cols");
    68   lp.clear();
    69   check(countRows(lp)==0, "Wrong number of rows");
    70   check(countCols(lp)==0, "Wrong number of cols");
    71   lp.addCol(); lp.addCol(); lp.addCol(); lp.addRow();
    72   check(countRows(lp)==1, "Wrong number of rows");
    73   check(countCols(lp)==3, "Wrong number of cols");
    74   lp.clear();
    75 
    76   std::vector<LP::Col> x(10);
    77   //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
    78   lp.addColSet(x);
    79   lp.colLowerBound(x,1);
    80   lp.colUpperBound(x,1);
    81   lp.colBounds(x,1,2);
    82 
    83   std::vector<LP::Col> y(10);
    84   lp.addColSet(y);
    85 
    86   lp.colLowerBound(y,1);
    87   lp.colUpperBound(y,1);
    88   lp.colBounds(y,1,2);
    89 
    90   std::map<int,LP::Col> z;
    91 
    92   z.insert(std::make_pair(12,INVALID));
    93   z.insert(std::make_pair(2,INVALID));
    94   z.insert(std::make_pair(7,INVALID));
    95   z.insert(std::make_pair(5,INVALID));
    96 
    97   lp.addColSet(z);
    98 
    99   lp.colLowerBound(z,1);
   100   lp.colUpperBound(z,1);
   101   lp.colBounds(z,1,2);
   102 
   103   {
   104     LP::Expr e,f,g;
   105     LP::Col p1,p2,p3,p4,p5;
   106     LP::Constr c;
   107 
   108     p1=lp.addCol();
   109     p2=lp.addCol();
   110     p3=lp.addCol();
   111     p4=lp.addCol();
   112     p5=lp.addCol();
   113 
   114     e[p1]=2;
   115     *e=12;
   116     e[p1]+=2;
   117     *e+=12;
   118     e[p1]-=2;
   119     *e-=12;
   120 
   121     e=2;
   122     e=2.2;
   123     e=p1;
   124     e=f;
   125 
   126     e+=2;
   127     e+=2.2;
   128     e+=p1;
   129     e+=f;
   130 
   131     e-=2;
   132     e-=2.2;
   133     e-=p1;
   134     e-=f;
   135 
   136     e*=2;
   137     e*=2.2;
   138     e/=2;
   139     e/=2.2;
   140 
   141     e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
   142        (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
   143        (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
   144        2.2*f+f*2.2+f/2.2+
   145        2*f+f*2+f/2+
   146        2.2*p1+p1*2.2+p1/2.2+
   147        2*p1+p1*2+p1/2
   148        );
   149 
   150 
   151     c = (e  <= f  );
   152     c = (e  <= 2.2);
   153     c = (e  <= 2  );
   154     c = (e  <= p1 );
   155     c = (2.2<= f  );
   156     c = (2  <= f  );
   157     c = (p1 <= f  );
   158     c = (p1 <= p2 );
   159     c = (p1 <= 2.2);
   160     c = (p1 <= 2  );
   161     c = (2.2<= p2 );
   162     c = (2  <= p2 );
   163 
   164     c = (e  >= f  );
   165     c = (e  >= 2.2);
   166     c = (e  >= 2  );
   167     c = (e  >= p1 );
   168     c = (2.2>= f  );
   169     c = (2  >= f  );
   170     c = (p1 >= f  );
   171     c = (p1 >= p2 );
   172     c = (p1 >= 2.2);
   173     c = (p1 >= 2  );
   174     c = (2.2>= p2 );
   175     c = (2  >= p2 );
   176 
   177     c = (e  == f  );
   178     c = (e  == 2.2);
   179     c = (e  == 2  );
   180     c = (e  == p1 );
   181     c = (2.2== f  );
   182     c = (2  == f  );
   183     c = (p1 == f  );
   184     //c = (p1 == p2 );
   185     c = (p1 == 2.2);
   186     c = (p1 == 2  );
   187     c = (2.2== p2 );
   188     c = (2  == p2 );
   189 
   190     c = ((2 <= e) <= 3);
   191     c = ((2 <= p1) <= 3);
   192 
   193     c = ((2 >= e) >= 3);
   194     c = ((2 >= p1) >= 3);
   195 
   196     { //Tests for #430
   197       LP::Col v=lp.addCol();
   198       LP::Constr c = v >= -3;
   199       c = c <= 4;
   200       LP::Constr c2;
   201       c2 = -3 <= v <= 4;
   202     }
   203 
   204     e[x[3]]=2;
   205     e[x[3]]=4;
   206     e[x[3]]=1;
   207     *e=12;
   208 
   209     lp.addRow(-LP::INF,e,23);
   210     lp.addRow(-LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
   211     lp.addRow(-LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
   212 
   213     lp.addRow(x[1]+x[3]<=x[5]-3);
   214     lp.addRow((-7<=x[1]+x[3]-12)<=3);
   215     lp.addRow(x[1]<=x[5]);
   216 
   217     std::ostringstream buf;
   218 
   219 
   220     e=((p1+p2)+(p1-0.99*p2));
   221     //e.prettyPrint(std::cout);
   222     //(e<=2).prettyPrint(std::cout);
   223     double tolerance=0.001;
   224     e.simplify(tolerance);
   225     buf << "Coeff. of p2 should be 0.01";
   226     check(e[p2]>0, buf.str());
   227 
   228     tolerance=0.02;
   229     e.simplify(tolerance);
   230     buf << "Coeff. of p2 should be 0";
   231     check(const_cast<const LpSolver::Expr&>(e)[p2]==0, buf.str());
   232 
   233     //Test for clone/new
   234     LP* lpnew = lp.newSolver();
   235     LP* lpclone = lp.cloneSolver();
   236     delete lpnew;
   237     delete lpclone;
   238 
   239   }
   240 
   241   {
   242     LP::DualExpr e,f,g;
   243     LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
   244       p4 = INVALID, p5 = INVALID;
   245 
   246     e[p1]=2;
   247     e[p1]+=2;
   248     e[p1]-=2;
   249 
   250     e=p1;
   251     e=f;
   252 
   253     e+=p1;
   254     e+=f;
   255 
   256     e-=p1;
   257     e-=f;
   258 
   259     e*=2;
   260     e*=2.2;
   261     e/=2;
   262     e/=2.2;
   263 
   264     e=((p1+p2)+(p1-p2)+
   265        (p1+f)+(f+p1)+(f+g)+
   266        (p1-f)+(f-p1)+(f-g)+
   267        2.2*f+f*2.2+f/2.2+
   268        2*f+f*2+f/2+
   269        2.2*p1+p1*2.2+p1/2.2+
   270        2*p1+p1*2+p1/2
   271        );
   272   }
   273 
   274 }
   275 
   276 void solveAndCheck(LpSolver& lp, LpSolver::ProblemType stat,
   277                    double exp_opt) {
   278   using std::string;
   279   lp.solve();
   280 
   281   std::ostringstream buf;
   282   buf << "PrimalType should be: " << int(stat) << int(lp.primalType());
   283 
   284   check(lp.primalType()==stat, buf.str());
   285 
   286   if (stat ==  LpSolver::OPTIMAL) {
   287     std::ostringstream sbuf;
   288     sbuf << "Wrong optimal value (" << lp.primal() <<") with "
   289          << lp.solverName() <<"\n     the right optimum is " << exp_opt;
   290     check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str());
   291   }
   292 }
   293 
   294 void aTest(LpSolver & lp)
   295 {
   296   typedef LpSolver LP;
   297 
   298  //The following example is very simple
   299 
   300   typedef LpSolver::Row Row;
   301   typedef LpSolver::Col Col;
   302 
   303 
   304   Col x1 = lp.addCol();
   305   Col x2 = lp.addCol();
   306 
   307 
   308   //Constraints
   309   Row upright=lp.addRow(x1+2*x2 <=1);
   310   lp.addRow(x1+x2 >=-1);
   311   lp.addRow(x1-x2 <=1);
   312   lp.addRow(x1-x2 >=-1);
   313   //Nonnegativity of the variables
   314   lp.colLowerBound(x1, 0);
   315   lp.colLowerBound(x2, 0);
   316   //Objective function
   317   lp.obj(x1+x2);
   318 
   319   lp.sense(lp.MAX);
   320 
   321   //Testing the problem retrieving routines
   322   check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
   323   check(lp.sense() == lp.MAX,"This is a maximization!");
   324   check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
   325   check(lp.colLowerBound(x1)==0,
   326         "The lower bound for variable x1 should be 0.");
   327   check(lp.colUpperBound(x1)==LpSolver::INF,
   328         "The upper bound for variable x1 should be infty.");
   329   check(lp.rowLowerBound(upright) == -LpSolver::INF,
   330         "The lower bound for the first row should be -infty.");
   331   check(lp.rowUpperBound(upright)==1,
   332         "The upper bound for the first row should be 1.");
   333   LpSolver::Expr e = lp.row(upright);
   334   check(e[x1] == 1, "The first coefficient should 1.");
   335   check(e[x2] == 2, "The second coefficient should 1.");
   336 
   337   lp.row(upright, x1+x2 <=1);
   338   e = lp.row(upright);
   339   check(e[x1] == 1, "The first coefficient should 1.");
   340   check(e[x2] == 1, "The second coefficient should 1.");
   341 
   342   LpSolver::DualExpr de = lp.col(x1);
   343   check(  de[upright] == 1, "The first coefficient should 1.");
   344 
   345   LpSolver* clp = lp.cloneSolver();
   346 
   347   //Testing the problem retrieving routines
   348   check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
   349   check(clp->sense() == clp->MAX,"This is a maximization!");
   350   check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
   351   //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   352   check(clp->colLowerBound(x1)==0,
   353         "The lower bound for variable x1 should be 0.");
   354   check(clp->colUpperBound(x1)==LpSolver::INF,
   355         "The upper bound for variable x1 should be infty.");
   356 
   357   check(lp.rowLowerBound(upright)==-LpSolver::INF,
   358         "The lower bound for the first row should be -infty.");
   359   check(lp.rowUpperBound(upright)==1,
   360         "The upper bound for the first row should be 1.");
   361   e = clp->row(upright);
   362   check(e[x1] == 1, "The first coefficient should 1.");
   363   check(e[x2] == 1, "The second coefficient should 1.");
   364 
   365   de = clp->col(x1);
   366   check(de[upright] == 1, "The first coefficient should 1.");
   367 
   368   delete clp;
   369 
   370   //Maximization of x1+x2
   371   //over the triangle with vertices (0,0) (0,1) (1,0)
   372   double expected_opt=1;
   373   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   374 
   375   //Minimization
   376   lp.sense(lp.MIN);
   377   expected_opt=0;
   378   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   379 
   380   //Vertex (-1,0) instead of (0,0)
   381   lp.colLowerBound(x1, -LpSolver::INF);
   382   expected_opt=-1;
   383   solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   384 
   385   //Erase one constraint and return to maximization
   386   lp.erase(upright);
   387   lp.sense(lp.MAX);
   388   expected_opt=LpSolver::INF;
   389   solveAndCheck(lp, LpSolver::UNBOUNDED, expected_opt);
   390 
   391   //Infeasibilty
   392   lp.addRow(x1+x2 <=-2);
   393   solveAndCheck(lp, LpSolver::INFEASIBLE, expected_opt);
   394 
   395 }
   396 
   397 template<class LP>
   398 void cloneTest()
   399 {
   400   //Test for clone/new
   401 
   402   LP* lp = new LP();
   403   LP* lpnew = lp->newSolver();
   404   LP* lpclone = lp->cloneSolver();
   405   delete lp;
   406   delete lpnew;
   407   delete lpclone;
   408 }
   409 
   410 int main()
   411 {
   412   LpSkeleton lp_skel;
   413   lpTest(lp_skel);
   414 
   415 #ifdef LEMON_HAVE_GLPK
   416   {
   417     GlpkLp lp_glpk1,lp_glpk2;
   418     lpTest(lp_glpk1);
   419     aTest(lp_glpk2);
   420     cloneTest<GlpkLp>();
   421   }
   422 #endif
   423 
   424 #ifdef LEMON_HAVE_CPLEX
   425   try {
   426     CplexLp lp_cplex1,lp_cplex2;
   427     lpTest(lp_cplex1);
   428     aTest(lp_cplex2);
   429     cloneTest<CplexLp>();
   430   } catch (CplexEnv::LicenseError& error) {
   431     check(false, error.what());
   432   }
   433 #endif
   434 
   435 #ifdef LEMON_HAVE_SOPLEX
   436   {
   437     SoplexLp lp_soplex1,lp_soplex2;
   438     lpTest(lp_soplex1);
   439     aTest(lp_soplex2);
   440     cloneTest<SoplexLp>();
   441   }
   442 #endif
   443 
   444 #ifdef LEMON_HAVE_CLP
   445   {
   446     ClpLp lp_clp1,lp_clp2;
   447     lpTest(lp_clp1);
   448     aTest(lp_clp2);
   449     cloneTest<ClpLp>();
   450   }
   451 #endif
   452 
   453   return 0;
   454 }