1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/test/lp_test.cc Tue Dec 02 21:40:33 2008 +0100
1.3 @@ -0,0 +1,423 @@
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-2008
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 +#ifdef HAVE_CONFIG_H
1.28 +#include <lemon/config.h>
1.29 +#endif
1.30 +
1.31 +#ifdef HAVE_GLPK
1.32 +#include <lemon/lp_glpk.h>
1.33 +#endif
1.34 +
1.35 +#ifdef HAVE_CPLEX
1.36 +#include <lemon/lp_cplex.h>
1.37 +#endif
1.38 +
1.39 +#ifdef HAVE_SOPLEX
1.40 +#include <lemon/lp_soplex.h>
1.41 +#endif
1.42 +
1.43 +using namespace lemon;
1.44 +
1.45 +void lpTest(LpSolverBase & lp)
1.46 +{
1.47 +
1.48 +
1.49 +
1.50 + typedef LpSolverBase 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 +#ifndef GYORSITAS
1.59 +
1.60 + std::vector<LP::Col> y(10);
1.61 + lp.addColSet(y);
1.62 +
1.63 + lp.colLowerBound(y,1);
1.64 + lp.colUpperBound(y,1);
1.65 + lp.colBounds(y,1,2);
1.66 +
1.67 + std::map<int,LP::Col> z;
1.68 +
1.69 + z.insert(std::make_pair(12,INVALID));
1.70 + z.insert(std::make_pair(2,INVALID));
1.71 + z.insert(std::make_pair(7,INVALID));
1.72 + z.insert(std::make_pair(5,INVALID));
1.73 +
1.74 + lp.addColSet(z);
1.75 +
1.76 + lp.colLowerBound(z,1);
1.77 + lp.colUpperBound(z,1);
1.78 + lp.colBounds(z,1,2);
1.79 +
1.80 + {
1.81 + LP::Expr e,f,g;
1.82 + LP::Col p1,p2,p3,p4,p5;
1.83 + LP::Constr c;
1.84 +
1.85 + p1=lp.addCol();
1.86 + p2=lp.addCol();
1.87 + p3=lp.addCol();
1.88 + p4=lp.addCol();
1.89 + p5=lp.addCol();
1.90 +
1.91 + e[p1]=2;
1.92 + e.constComp()=12;
1.93 + e[p1]+=2;
1.94 + e.constComp()+=12;
1.95 + e[p1]-=2;
1.96 + e.constComp()-=12;
1.97 +
1.98 + e=2;
1.99 + e=2.2;
1.100 + e=p1;
1.101 + e=f;
1.102 +
1.103 + e+=2;
1.104 + e+=2.2;
1.105 + e+=p1;
1.106 + e+=f;
1.107 +
1.108 + e-=2;
1.109 + e-=2.2;
1.110 + e-=p1;
1.111 + e-=f;
1.112 +
1.113 + e*=2;
1.114 + e*=2.2;
1.115 + e/=2;
1.116 + e/=2.2;
1.117 +
1.118 + e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
1.119 + (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
1.120 + (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
1.121 + 2.2*f+f*2.2+f/2.2+
1.122 + 2*f+f*2+f/2+
1.123 + 2.2*p1+p1*2.2+p1/2.2+
1.124 + 2*p1+p1*2+p1/2
1.125 + );
1.126 +
1.127 +
1.128 + c = (e <= f );
1.129 + c = (e <= 2.2);
1.130 + c = (e <= 2 );
1.131 + c = (e <= p1 );
1.132 + c = (2.2<= f );
1.133 + c = (2 <= f );
1.134 + c = (p1 <= f );
1.135 + c = (p1 <= p2 );
1.136 + c = (p1 <= 2.2);
1.137 + c = (p1 <= 2 );
1.138 + c = (2.2<= p2 );
1.139 + c = (2 <= p2 );
1.140 +
1.141 + c = (e >= f );
1.142 + c = (e >= 2.2);
1.143 + c = (e >= 2 );
1.144 + c = (e >= p1 );
1.145 + c = (2.2>= f );
1.146 + c = (2 >= f );
1.147 + c = (p1 >= f );
1.148 + c = (p1 >= p2 );
1.149 + c = (p1 >= 2.2);
1.150 + c = (p1 >= 2 );
1.151 + c = (2.2>= p2 );
1.152 + c = (2 >= p2 );
1.153 +
1.154 + c = (e == f );
1.155 + c = (e == 2.2);
1.156 + c = (e == 2 );
1.157 + c = (e == p1 );
1.158 + c = (2.2== f );
1.159 + c = (2 == f );
1.160 + c = (p1 == f );
1.161 + //c = (p1 == p2 );
1.162 + c = (p1 == 2.2);
1.163 + c = (p1 == 2 );
1.164 + c = (2.2== p2 );
1.165 + c = (2 == p2 );
1.166 +
1.167 + c = (2 <= e <= 3);
1.168 + c = (2 <= p1<= 3);
1.169 +
1.170 + c = (2 >= e >= 3);
1.171 + c = (2 >= p1>= 3);
1.172 +
1.173 + e[x[3]]=2;
1.174 + e[x[3]]=4;
1.175 + e[x[3]]=1;
1.176 + e.constComp()=12;
1.177 +
1.178 + lp.addRow(LP::INF,e,23);
1.179 + lp.addRow(LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
1.180 + lp.addRow(LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
1.181 +
1.182 + lp.addRow(x[1]+x[3]<=x[5]-3);
1.183 + lp.addRow(-7<=x[1]+x[3]-12<=3);
1.184 + lp.addRow(x[1]<=x[5]);
1.185 +
1.186 + std::ostringstream buf;
1.187 +
1.188 +
1.189 + //Checking the simplify function
1.190 +
1.191 +// //How to check the simplify function? A map gives no information
1.192 +// //on the question whether a given key is or is not stored in it, or
1.193 +// //it does?
1.194 +// Yes, it does, using the find() function.
1.195 + e=((p1+p2)+(p1-p2));
1.196 + e.simplify();
1.197 + buf << "Coeff. of p2 should be 0";
1.198 + // std::cout<<e[p1]<<e[p2]<<e[p3]<<std::endl;
1.199 + check(e.find(p2)==e.end(), buf.str());
1.200 +
1.201 +
1.202 +
1.203 +
1.204 + e=((p1+p2)+(p1-0.99*p2));
1.205 + //e.prettyPrint(std::cout);
1.206 + //(e<=2).prettyPrint(std::cout);
1.207 + double tolerance=0.001;
1.208 + e.simplify(tolerance);
1.209 + buf << "Coeff. of p2 should be 0.01";
1.210 + check(e[p2]>0, buf.str());
1.211 +
1.212 + tolerance=0.02;
1.213 + e.simplify(tolerance);
1.214 + buf << "Coeff. of p2 should be 0";
1.215 + check(e.find(p2)==e.end(), buf.str());
1.216 +
1.217 +
1.218 + }
1.219 +
1.220 + {
1.221 + LP::DualExpr e,f,g;
1.222 + LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
1.223 + p4 = INVALID, p5 = INVALID;
1.224 +
1.225 + e[p1]=2;
1.226 + e[p1]+=2;
1.227 + e[p1]-=2;
1.228 +
1.229 + e=p1;
1.230 + e=f;
1.231 +
1.232 + e+=p1;
1.233 + e+=f;
1.234 +
1.235 + e-=p1;
1.236 + e-=f;
1.237 +
1.238 + e*=2;
1.239 + e*=2.2;
1.240 + e/=2;
1.241 + e/=2.2;
1.242 +
1.243 + e=((p1+p2)+(p1-p2)+
1.244 + (p1+f)+(f+p1)+(f+g)+
1.245 + (p1-f)+(f-p1)+(f-g)+
1.246 + 2.2*f+f*2.2+f/2.2+
1.247 + 2*f+f*2+f/2+
1.248 + 2.2*p1+p1*2.2+p1/2.2+
1.249 + 2*p1+p1*2+p1/2
1.250 + );
1.251 + }
1.252 +
1.253 +#endif
1.254 +}
1.255 +
1.256 +void solveAndCheck(LpSolverBase& lp, LpSolverBase::SolutionStatus stat,
1.257 + double exp_opt) {
1.258 + using std::string;
1.259 + lp.solve();
1.260 + //int decimal,sign;
1.261 + std::ostringstream buf;
1.262 + buf << "Primalstatus should be: " << int(stat);
1.263 +
1.264 + // itoa(stat,buf1, 10);
1.265 + check(lp.primalStatus()==stat, buf.str());
1.266 +
1.267 + if (stat == LpSolverBase::OPTIMAL) {
1.268 + std::ostringstream sbuf;
1.269 + sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
1.270 + check(std::abs(lp.primalValue()-exp_opt) < 1e-3, sbuf.str());
1.271 + //+ecvt(exp_opt,2)
1.272 + }
1.273 +}
1.274 +
1.275 +void aTest(LpSolverBase & lp)
1.276 +{
1.277 + typedef LpSolverBase LP;
1.278 +
1.279 + //The following example is very simple
1.280 +
1.281 + typedef LpSolverBase::Row Row;
1.282 + typedef LpSolverBase::Col Col;
1.283 +
1.284 +
1.285 + Col x1 = lp.addCol();
1.286 + Col x2 = lp.addCol();
1.287 +
1.288 +
1.289 + //Constraints
1.290 + Row upright=lp.addRow(x1+x2 <=1);
1.291 + lp.addRow(x1+x2 >=-1);
1.292 + lp.addRow(x1-x2 <=1);
1.293 + lp.addRow(x1-x2 >=-1);
1.294 + //Nonnegativity of the variables
1.295 + lp.colLowerBound(x1, 0);
1.296 + lp.colLowerBound(x2, 0);
1.297 + //Objective function
1.298 + lp.obj(x1+x2);
1.299 +
1.300 + lp.max();
1.301 +
1.302 + //Testing the problem retrieving routines
1.303 + check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
1.304 + check(lp.isMax(),"This is a maximization!");
1.305 + check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
1.306 + // std::cout<<lp.colLowerBound(x1)<<std::endl;
1.307 + check( lp.colLowerBound(x1)==0,
1.308 + "The lower bound for variable x1 should be 0.");
1.309 + check( lp.colUpperBound(x1)==LpSolverBase::INF,
1.310 + "The upper bound for variable x1 should be infty.");
1.311 + LpSolverBase::Value lb,ub;
1.312 + lp.getRowBounds(upright,lb,ub);
1.313 + check( lb==-LpSolverBase::INF,
1.314 + "The lower bound for the first row should be -infty.");
1.315 + check( ub==1,"The upper bound for the first row should be 1.");
1.316 + LpSolverBase::Expr e = lp.row(upright);
1.317 + check( e.size() == 2, "The row retrieval gives back wrong expression.");
1.318 + check( e[x1] == 1, "The first coefficient should 1.");
1.319 + check( e[x2] == 1, "The second coefficient should 1.");
1.320 +
1.321 + LpSolverBase::DualExpr de = lp.col(x1);
1.322 + check( de.size() == 4, "The col retrieval gives back wrong expression.");
1.323 + check( de[upright] == 1, "The first coefficient should 1.");
1.324 +
1.325 + LpSolverBase* clp = lp.copyLp();
1.326 +
1.327 + //Testing the problem retrieving routines
1.328 + check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
1.329 + check(clp->isMax(),"This is a maximization!");
1.330 + check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
1.331 + // std::cout<<lp.colLowerBound(x1)<<std::endl;
1.332 + check( clp->colLowerBound(x1)==0,
1.333 + "The lower bound for variable x1 should be 0.");
1.334 + check( clp->colUpperBound(x1)==LpSolverBase::INF,
1.335 + "The upper bound for variable x1 should be infty.");
1.336 +
1.337 + clp->getRowBounds(upright,lb,ub);
1.338 + check( lb==-LpSolverBase::INF,
1.339 + "The lower bound for the first row should be -infty.");
1.340 + check( ub==1,"The upper bound for the first row should be 1.");
1.341 + e = clp->row(upright);
1.342 + check( e.size() == 2, "The row retrieval gives back wrong expression.");
1.343 + check( e[x1] == 1, "The first coefficient should 1.");
1.344 + check( e[x2] == 1, "The second coefficient should 1.");
1.345 +
1.346 + de = clp->col(x1);
1.347 + check( de.size() == 4, "The col retrieval gives back wrong expression.");
1.348 + check( de[upright] == 1, "The first coefficient should 1.");
1.349 +
1.350 + delete clp;
1.351 +
1.352 + //Maximization of x1+x2
1.353 + //over the triangle with vertices (0,0) (0,1) (1,0)
1.354 + double expected_opt=1;
1.355 + solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
1.356 +
1.357 + //Minimization
1.358 + lp.min();
1.359 + expected_opt=0;
1.360 + solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
1.361 +
1.362 + //Vertex (-1,0) instead of (0,0)
1.363 + lp.colLowerBound(x1, -LpSolverBase::INF);
1.364 + expected_opt=-1;
1.365 + solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
1.366 +
1.367 + //Erase one constraint and return to maximization
1.368 + lp.eraseRow(upright);
1.369 + lp.max();
1.370 + expected_opt=LpSolverBase::INF;
1.371 + solveAndCheck(lp, LpSolverBase::INFINITE, expected_opt);
1.372 +
1.373 + //Infeasibilty
1.374 + lp.addRow(x1+x2 <=-2);
1.375 + solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt);
1.376 +
1.377 + //Change problem and forget to solve
1.378 + lp.min();
1.379 + check(lp.primalStatus()==LpSolverBase::UNDEFINED,
1.380 + "Primalstatus should be UNDEFINED");
1.381 +
1.382 +
1.383 +// lp.solve();
1.384 +// if (lp.primalStatus()==LpSolverBase::OPTIMAL){
1.385 +// std::cout<< "Z = "<<lp.primalValue()
1.386 +// << " (error = " << lp.primalValue()-expected_opt
1.387 +// << "); x1 = "<<lp.primal(x1)
1.388 +// << "; x2 = "<<lp.primal(x2)
1.389 +// <<std::endl;
1.390 +
1.391 +// }
1.392 +// else{
1.393 +// std::cout<<lp.primalStatus()<<std::endl;
1.394 +// std::cout<<"Optimal solution not found!"<<std::endl;
1.395 +// }
1.396 +
1.397 +
1.398 +
1.399 +}
1.400 +
1.401 +
1.402 +int main()
1.403 +{
1.404 + LpSkeleton lp_skel;
1.405 + lpTest(lp_skel);
1.406 +
1.407 +#ifdef HAVE_GLPK
1.408 + LpGlpk lp_glpk1,lp_glpk2;
1.409 + lpTest(lp_glpk1);
1.410 + aTest(lp_glpk2);
1.411 +#endif
1.412 +
1.413 +#ifdef HAVE_CPLEX
1.414 + LpCplex lp_cplex1,lp_cplex2;
1.415 + lpTest(lp_cplex1);
1.416 + aTest(lp_cplex2);
1.417 +#endif
1.418 +
1.419 +#ifdef HAVE_SOPLEX
1.420 + LpSoplex lp_soplex1,lp_soplex2;
1.421 + lpTest(lp_soplex1);
1.422 + aTest(lp_soplex2);
1.423 +#endif
1.424 +
1.425 + return 0;
1.426 +}