test/lp_test.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 03 Apr 2009 18:59:15 +0200
changeset 608 6ac5d9ae1d3d
parent 461 08d495d48089
child 537 0fec6a017ead
child 540 9db62975c32b
permissions -rw-r--r--
Support real types + numerical stability fix in NS (#254)

- Real types are supported by appropriate inicialization.
- A feature of the XTI spanning tree structure is removed to ensure
numerical stability (could cause problems using integer types).
The node potentials are updated always on the lower subtree,
in order to prevent overflow problems.
The former method isn't notably faster during to our tests.
     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 }