1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/test/lp_test.cc Thu Nov 05 15:50:01 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 +}