test/lp_test.cc
changeset 867 994c7df296c9
parent 674 20dac2104519
child 1092 2eebc8f7dca5
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/test/lp_test.cc	Thu Dec 10 17:05:35 2009 +0100
     1.3 @@ -0,0 +1,419 @@
     1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     1.5 + *
     1.6 + * This file is a part of LEMON, a generic C++ optimization library.
     1.7 + *
     1.8 + * Copyright (C) 2003-2009
     1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    1.11 + *
    1.12 + * Permission to use, modify and distribute this software is granted
    1.13 + * provided that this copyright notice appears in all copies. For
    1.14 + * precise terms see the accompanying LICENSE file.
    1.15 + *
    1.16 + * This software is provided "AS IS" with no warranty of any kind,
    1.17 + * express or implied, and with no claim as to its suitability for any
    1.18 + * purpose.
    1.19 + *
    1.20 + */
    1.21 +
    1.22 +#include <sstream>
    1.23 +#include <lemon/lp_skeleton.h>
    1.24 +#include "test_tools.h"
    1.25 +#include <lemon/tolerance.h>
    1.26 +
    1.27 +#include <lemon/config.h>
    1.28 +
    1.29 +#ifdef LEMON_HAVE_GLPK
    1.30 +#include <lemon/glpk.h>
    1.31 +#endif
    1.32 +
    1.33 +#ifdef LEMON_HAVE_CPLEX
    1.34 +#include <lemon/cplex.h>
    1.35 +#endif
    1.36 +
    1.37 +#ifdef LEMON_HAVE_SOPLEX
    1.38 +#include <lemon/soplex.h>
    1.39 +#endif
    1.40 +
    1.41 +#ifdef LEMON_HAVE_CLP
    1.42 +#include <lemon/clp.h>
    1.43 +#endif
    1.44 +
    1.45 +using namespace lemon;
    1.46 +
    1.47 +void lpTest(LpSolver& lp)
    1.48 +{
    1.49 +
    1.50 +  typedef LpSolver LP;
    1.51 +
    1.52 +  std::vector<LP::Col> x(10);
    1.53 +  //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
    1.54 +  lp.addColSet(x);
    1.55 +  lp.colLowerBound(x,1);
    1.56 +  lp.colUpperBound(x,1);
    1.57 +  lp.colBounds(x,1,2);
    1.58 +
    1.59 +  std::vector<LP::Col> y(10);
    1.60 +  lp.addColSet(y);
    1.61 +
    1.62 +  lp.colLowerBound(y,1);
    1.63 +  lp.colUpperBound(y,1);
    1.64 +  lp.colBounds(y,1,2);
    1.65 +
    1.66 +  std::map<int,LP::Col> z;
    1.67 +
    1.68 +  z.insert(std::make_pair(12,INVALID));
    1.69 +  z.insert(std::make_pair(2,INVALID));
    1.70 +  z.insert(std::make_pair(7,INVALID));
    1.71 +  z.insert(std::make_pair(5,INVALID));
    1.72 +
    1.73 +  lp.addColSet(z);
    1.74 +
    1.75 +  lp.colLowerBound(z,1);
    1.76 +  lp.colUpperBound(z,1);
    1.77 +  lp.colBounds(z,1,2);
    1.78 +
    1.79 +  {
    1.80 +    LP::Expr e,f,g;
    1.81 +    LP::Col p1,p2,p3,p4,p5;
    1.82 +    LP::Constr c;
    1.83 +
    1.84 +    p1=lp.addCol();
    1.85 +    p2=lp.addCol();
    1.86 +    p3=lp.addCol();
    1.87 +    p4=lp.addCol();
    1.88 +    p5=lp.addCol();
    1.89 +
    1.90 +    e[p1]=2;
    1.91 +    *e=12;
    1.92 +    e[p1]+=2;
    1.93 +    *e+=12;
    1.94 +    e[p1]-=2;
    1.95 +    *e-=12;
    1.96 +
    1.97 +    e=2;
    1.98 +    e=2.2;
    1.99 +    e=p1;
   1.100 +    e=f;
   1.101 +
   1.102 +    e+=2;
   1.103 +    e+=2.2;
   1.104 +    e+=p1;
   1.105 +    e+=f;
   1.106 +
   1.107 +    e-=2;
   1.108 +    e-=2.2;
   1.109 +    e-=p1;
   1.110 +    e-=f;
   1.111 +
   1.112 +    e*=2;
   1.113 +    e*=2.2;
   1.114 +    e/=2;
   1.115 +    e/=2.2;
   1.116 +
   1.117 +    e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
   1.118 +       (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
   1.119 +       (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
   1.120 +       2.2*f+f*2.2+f/2.2+
   1.121 +       2*f+f*2+f/2+
   1.122 +       2.2*p1+p1*2.2+p1/2.2+
   1.123 +       2*p1+p1*2+p1/2
   1.124 +       );
   1.125 +
   1.126 +
   1.127 +    c = (e  <= f  );
   1.128 +    c = (e  <= 2.2);
   1.129 +    c = (e  <= 2  );
   1.130 +    c = (e  <= p1 );
   1.131 +    c = (2.2<= f  );
   1.132 +    c = (2  <= f  );
   1.133 +    c = (p1 <= f  );
   1.134 +    c = (p1 <= p2 );
   1.135 +    c = (p1 <= 2.2);
   1.136 +    c = (p1 <= 2  );
   1.137 +    c = (2.2<= p2 );
   1.138 +    c = (2  <= p2 );
   1.139 +
   1.140 +    c = (e  >= f  );
   1.141 +    c = (e  >= 2.2);
   1.142 +    c = (e  >= 2  );
   1.143 +    c = (e  >= p1 );
   1.144 +    c = (2.2>= f  );
   1.145 +    c = (2  >= f  );
   1.146 +    c = (p1 >= f  );
   1.147 +    c = (p1 >= p2 );
   1.148 +    c = (p1 >= 2.2);
   1.149 +    c = (p1 >= 2  );
   1.150 +    c = (2.2>= p2 );
   1.151 +    c = (2  >= p2 );
   1.152 +
   1.153 +    c = (e  == f  );
   1.154 +    c = (e  == 2.2);
   1.155 +    c = (e  == 2  );
   1.156 +    c = (e  == p1 );
   1.157 +    c = (2.2== f  );
   1.158 +    c = (2  == f  );
   1.159 +    c = (p1 == f  );
   1.160 +    //c = (p1 == p2 );
   1.161 +    c = (p1 == 2.2);
   1.162 +    c = (p1 == 2  );
   1.163 +    c = (2.2== p2 );
   1.164 +    c = (2  == p2 );
   1.165 +
   1.166 +    c = ((2 <= e) <= 3);
   1.167 +    c = ((2 <= p1) <= 3);
   1.168 +
   1.169 +    c = ((2 >= e) >= 3);
   1.170 +    c = ((2 >= p1) >= 3);
   1.171 +
   1.172 +    e[x[3]]=2;
   1.173 +    e[x[3]]=4;
   1.174 +    e[x[3]]=1;
   1.175 +    *e=12;
   1.176 +
   1.177 +    lp.addRow(-LP::INF,e,23);
   1.178 +    lp.addRow(-LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
   1.179 +    lp.addRow(-LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
   1.180 +
   1.181 +    lp.addRow(x[1]+x[3]<=x[5]-3);
   1.182 +    lp.addRow((-7<=x[1]+x[3]-12)<=3);
   1.183 +    lp.addRow(x[1]<=x[5]);
   1.184 +
   1.185 +    std::ostringstream buf;
   1.186 +
   1.187 +
   1.188 +    e=((p1+p2)+(p1-0.99*p2));
   1.189 +    //e.prettyPrint(std::cout);
   1.190 +    //(e<=2).prettyPrint(std::cout);
   1.191 +    double tolerance=0.001;
   1.192 +    e.simplify(tolerance);
   1.193 +    buf << "Coeff. of p2 should be 0.01";
   1.194 +    check(e[p2]>0, buf.str());
   1.195 +
   1.196 +    tolerance=0.02;
   1.197 +    e.simplify(tolerance);
   1.198 +    buf << "Coeff. of p2 should be 0";
   1.199 +    check(const_cast<const LpSolver::Expr&>(e)[p2]==0, buf.str());
   1.200 +
   1.201 +    //Test for clone/new
   1.202 +    LP* lpnew = lp.newSolver();
   1.203 +    LP* lpclone = lp.cloneSolver();
   1.204 +    delete lpnew;
   1.205 +    delete lpclone;
   1.206 +
   1.207 +  }
   1.208 +
   1.209 +  {
   1.210 +    LP::DualExpr e,f,g;
   1.211 +    LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
   1.212 +      p4 = INVALID, p5 = INVALID;
   1.213 +
   1.214 +    e[p1]=2;
   1.215 +    e[p1]+=2;
   1.216 +    e[p1]-=2;
   1.217 +
   1.218 +    e=p1;
   1.219 +    e=f;
   1.220 +
   1.221 +    e+=p1;
   1.222 +    e+=f;
   1.223 +
   1.224 +    e-=p1;
   1.225 +    e-=f;
   1.226 +
   1.227 +    e*=2;
   1.228 +    e*=2.2;
   1.229 +    e/=2;
   1.230 +    e/=2.2;
   1.231 +
   1.232 +    e=((p1+p2)+(p1-p2)+
   1.233 +       (p1+f)+(f+p1)+(f+g)+
   1.234 +       (p1-f)+(f-p1)+(f-g)+
   1.235 +       2.2*f+f*2.2+f/2.2+
   1.236 +       2*f+f*2+f/2+
   1.237 +       2.2*p1+p1*2.2+p1/2.2+
   1.238 +       2*p1+p1*2+p1/2
   1.239 +       );
   1.240 +  }
   1.241 +
   1.242 +}
   1.243 +
   1.244 +void solveAndCheck(LpSolver& lp, LpSolver::ProblemType stat,
   1.245 +                   double exp_opt) {
   1.246 +  using std::string;
   1.247 +  lp.solve();
   1.248 +
   1.249 +  std::ostringstream buf;
   1.250 +  buf << "PrimalType should be: " << int(stat) << int(lp.primalType());
   1.251 +
   1.252 +  check(lp.primalType()==stat, buf.str());
   1.253 +
   1.254 +  if (stat ==  LpSolver::OPTIMAL) {
   1.255 +    std::ostringstream sbuf;
   1.256 +    sbuf << "Wrong optimal value (" << lp.primal() <<") with "
   1.257 +         << lp.solverName() <<"\n     the right optimum is " << exp_opt;
   1.258 +    check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str());
   1.259 +  }
   1.260 +}
   1.261 +
   1.262 +void aTest(LpSolver & lp)
   1.263 +{
   1.264 +  typedef LpSolver LP;
   1.265 +
   1.266 + //The following example is very simple
   1.267 +
   1.268 +  typedef LpSolver::Row Row;
   1.269 +  typedef LpSolver::Col Col;
   1.270 +
   1.271 +
   1.272 +  Col x1 = lp.addCol();
   1.273 +  Col x2 = lp.addCol();
   1.274 +
   1.275 +
   1.276 +  //Constraints
   1.277 +  Row upright=lp.addRow(x1+2*x2 <=1);
   1.278 +  lp.addRow(x1+x2 >=-1);
   1.279 +  lp.addRow(x1-x2 <=1);
   1.280 +  lp.addRow(x1-x2 >=-1);
   1.281 +  //Nonnegativity of the variables
   1.282 +  lp.colLowerBound(x1, 0);
   1.283 +  lp.colLowerBound(x2, 0);
   1.284 +  //Objective function
   1.285 +  lp.obj(x1+x2);
   1.286 +
   1.287 +  lp.sense(lp.MAX);
   1.288 +
   1.289 +  //Testing the problem retrieving routines
   1.290 +  check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
   1.291 +  check(lp.sense() == lp.MAX,"This is a maximization!");
   1.292 +  check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
   1.293 +  check(lp.colLowerBound(x1)==0,
   1.294 +        "The lower bound for variable x1 should be 0.");
   1.295 +  check(lp.colUpperBound(x1)==LpSolver::INF,
   1.296 +        "The upper bound for variable x1 should be infty.");
   1.297 +  check(lp.rowLowerBound(upright) == -LpSolver::INF,
   1.298 +        "The lower bound for the first row should be -infty.");
   1.299 +  check(lp.rowUpperBound(upright)==1,
   1.300 +        "The upper bound for the first row should be 1.");
   1.301 +  LpSolver::Expr e = lp.row(upright);
   1.302 +  check(e[x1] == 1, "The first coefficient should 1.");
   1.303 +  check(e[x2] == 2, "The second coefficient should 1.");
   1.304 +
   1.305 +  lp.row(upright, x1+x2 <=1);
   1.306 +  e = lp.row(upright);
   1.307 +  check(e[x1] == 1, "The first coefficient should 1.");
   1.308 +  check(e[x2] == 1, "The second coefficient should 1.");
   1.309 +
   1.310 +  LpSolver::DualExpr de = lp.col(x1);
   1.311 +  check(  de[upright] == 1, "The first coefficient should 1.");
   1.312 +
   1.313 +  LpSolver* clp = lp.cloneSolver();
   1.314 +
   1.315 +  //Testing the problem retrieving routines
   1.316 +  check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
   1.317 +  check(clp->sense() == clp->MAX,"This is a maximization!");
   1.318 +  check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
   1.319 +  //  std::cout<<lp.colLowerBound(x1)<<std::endl;
   1.320 +  check(clp->colLowerBound(x1)==0,
   1.321 +        "The lower bound for variable x1 should be 0.");
   1.322 +  check(clp->colUpperBound(x1)==LpSolver::INF,
   1.323 +        "The upper bound for variable x1 should be infty.");
   1.324 +
   1.325 +  check(lp.rowLowerBound(upright)==-LpSolver::INF,
   1.326 +        "The lower bound for the first row should be -infty.");
   1.327 +  check(lp.rowUpperBound(upright)==1,
   1.328 +        "The upper bound for the first row should be 1.");
   1.329 +  e = clp->row(upright);
   1.330 +  check(e[x1] == 1, "The first coefficient should 1.");
   1.331 +  check(e[x2] == 1, "The second coefficient should 1.");
   1.332 +
   1.333 +  de = clp->col(x1);
   1.334 +  check(de[upright] == 1, "The first coefficient should 1.");
   1.335 +
   1.336 +  delete clp;
   1.337 +
   1.338 +  //Maximization of x1+x2
   1.339 +  //over the triangle with vertices (0,0) (0,1) (1,0)
   1.340 +  double expected_opt=1;
   1.341 +  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   1.342 +
   1.343 +  //Minimization
   1.344 +  lp.sense(lp.MIN);
   1.345 +  expected_opt=0;
   1.346 +  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   1.347 +
   1.348 +  //Vertex (-1,0) instead of (0,0)
   1.349 +  lp.colLowerBound(x1, -LpSolver::INF);
   1.350 +  expected_opt=-1;
   1.351 +  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
   1.352 +
   1.353 +  //Erase one constraint and return to maximization
   1.354 +  lp.erase(upright);
   1.355 +  lp.sense(lp.MAX);
   1.356 +  expected_opt=LpSolver::INF;
   1.357 +  solveAndCheck(lp, LpSolver::UNBOUNDED, expected_opt);
   1.358 +
   1.359 +  //Infeasibilty
   1.360 +  lp.addRow(x1+x2 <=-2);
   1.361 +  solveAndCheck(lp, LpSolver::INFEASIBLE, expected_opt);
   1.362 +
   1.363 +}
   1.364 +
   1.365 +template<class LP>
   1.366 +void cloneTest()
   1.367 +{
   1.368 +  //Test for clone/new
   1.369 +
   1.370 +  LP* lp = new LP();
   1.371 +  LP* lpnew = lp->newSolver();
   1.372 +  LP* lpclone = lp->cloneSolver();
   1.373 +  delete lp;
   1.374 +  delete lpnew;
   1.375 +  delete lpclone;
   1.376 +}
   1.377 +
   1.378 +int main()
   1.379 +{
   1.380 +  LpSkeleton lp_skel;
   1.381 +  lpTest(lp_skel);
   1.382 +
   1.383 +#ifdef LEMON_HAVE_GLPK
   1.384 +  {
   1.385 +    GlpkLp lp_glpk1,lp_glpk2;
   1.386 +    lpTest(lp_glpk1);
   1.387 +    aTest(lp_glpk2);
   1.388 +    cloneTest<GlpkLp>();
   1.389 +  }
   1.390 +#endif
   1.391 +
   1.392 +#ifdef LEMON_HAVE_CPLEX
   1.393 +  try {
   1.394 +    CplexLp lp_cplex1,lp_cplex2;
   1.395 +    lpTest(lp_cplex1);
   1.396 +    aTest(lp_cplex2);
   1.397 +    cloneTest<CplexLp>();
   1.398 +  } catch (CplexEnv::LicenseError& error) {
   1.399 +    check(false, error.what());
   1.400 +  }
   1.401 +#endif
   1.402 +
   1.403 +#ifdef LEMON_HAVE_SOPLEX
   1.404 +  {
   1.405 +    SoplexLp lp_soplex1,lp_soplex2;
   1.406 +    lpTest(lp_soplex1);
   1.407 +    aTest(lp_soplex2);
   1.408 +    cloneTest<SoplexLp>();
   1.409 +  }
   1.410 +#endif
   1.411 +
   1.412 +#ifdef LEMON_HAVE_CLP
   1.413 +  {
   1.414 +    ClpLp lp_clp1,lp_clp2;
   1.415 +    lpTest(lp_clp1);
   1.416 +    aTest(lp_clp2);
   1.417 +    cloneTest<ClpLp>();
   1.418 +  }
   1.419 +#endif
   1.420 +
   1.421 +  return 0;
   1.422 +}