diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -50,9 +50,9 @@ AC_SUBST([WARNINGCXXFLAGS]) dnl Checks for libraries. -#LX_CHECK_GLPK -#LX_CHECK_CPLEX -#LX_CHECK_SOPLEX +LX_CHECK_GLPK +LX_CHECK_CPLEX +LX_CHECK_SOPLEX AM_CONDITIONAL([HAVE_LP], [test x"$lx_lp_found" = x"yes"]) AM_CONDITIONAL([HAVE_MIP], [test x"$lx_mip_found" = x"yes"]) @@ -117,10 +117,10 @@ echo C++ compiler.................. : $CXX echo C++ compiles flags............ : $WARNINGCXXFLAGS $CXXFLAGS echo -#echo GLPK support.................. : $lx_glpk_found -#echo CPLEX support................. : $lx_cplex_found -#echo SOPLEX support................ : $lx_soplex_found -#echo +echo GLPK support.................. : $lx_glpk_found +echo CPLEX support................. : $lx_cplex_found +echo SOPLEX support................ : $lx_soplex_found +echo echo Build demo programs........... : $enable_demo echo Build additional tools........ : $enable_tools echo diff --git a/lemon/CMakeLists.txt b/lemon/CMakeLists.txt --- a/lemon/CMakeLists.txt +++ b/lemon/CMakeLists.txt @@ -4,6 +4,9 @@ arg_parser.cc base.cc color.cc + lp_base.cc + lp_skeleton.cc + lp_utils.cc random.cc) INSTALL( diff --git a/lemon/Makefile.am b/lemon/Makefile.am --- a/lemon/Makefile.am +++ b/lemon/Makefile.am @@ -10,10 +10,32 @@ lemon/arg_parser.cc \ lemon/base.cc \ lemon/color.cc \ + lemon/lp_base.cc \ + lemon/lp_skeleton.cc \ lemon/random.cc -#lemon_libemon_la_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS) $(SOPLEX_CXXFLAGS) $(AM_CXXFLAGS) -#lemon_libemon_la_LDFLAGS = $(GLPK_LIBS) $(CPLEX_LIBS) $(SOPLEX_LIBS) + +lemon_libemon_la_CXXFLAGS = \ + $(GLPK_CFLAGS) \ + $(CPLEX_CFLAGS) \ + $(SOPLEX_CXXFLAGS) + +lemon_libemon_la_LDFLAGS = \ + $(GLPK_LIBS) \ + $(CPLEX_LIBS) \ + $(SOPLEX_LIBS) + +if HAVE_GLPK +lemon_libemon_la_SOURCES += lemon/lp_glpk.cc lemon/mip_glpk.cc +endif + +if HAVE_CPLEX +lemon_libemon_la_SOURCES += lemon/lp_cplex.cc lemon/mip_cplex.cc +endif + +if HAVE_SOPLEX +lemon_libemon_la_SOURCES += lemon/lp_soplex.cc +endif lemon_HEADERS += \ lemon/adaptors.h \ @@ -41,6 +63,14 @@ lemon/lgf_reader.h \ lemon/lgf_writer.h \ lemon/list_graph.h \ + lemon/lp.h \ + lemon/lp_base.h \ + lemon/lp_cplex.h \ + lemon/lp_glpk.h \ + lemon/lp_skeleton.h \ + lemon/lp_soplex.h \ + lemon/mip_cplex.h \ + lemon/mip_glpk.h \ lemon/maps.h \ lemon/math.h \ lemon/max_matching.h \ @@ -64,6 +94,7 @@ lemon/bits/enable_if.h \ lemon/bits/graph_adaptor_extender.h \ lemon/bits/graph_extender.h \ + lemon/bits/lp_id.h \ lemon/bits/map_extender.h \ lemon/bits/path_dump.h \ lemon/bits/traits.h \ diff --git a/lemon/bits/lp_id.h b/lemon/bits/lp_id.h new file mode 100644 --- /dev/null +++ b/lemon/bits/lp_id.h @@ -0,0 +1,157 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_BITS_LP_SOLVER_ID_H +#define LEMON_BITS_LP_SOLVER_ID_H + +namespace lemon { + + namespace _lp_bits { + + struct LpIdImpl { + std::vector index; + std::vector cross; + int first_index; + int first_free; + }; + + class LpId { + public: + + class IdHandler { + public: + virtual ~IdHandler() {} + virtual int addId(LpIdImpl&) = 0; + virtual void eraseId(LpIdImpl&, int xn) = 0; + }; + + LpId(int min_index = 0) { + id_handler = 0; + impl.first_free = -1; + impl.first_index = min_index; + impl.cross.resize(impl.first_index); + } + + LpId(const LpId& li) { + id_handler = 0; + impl = li.impl; + } + + LpId& operator=(const LpId& li) { + id_handler = 0; + impl = li.impl; + return *this; + } + + void setIdHandler(IdHandler& ih) { + id_handler = &ih; + } + + int fixId(int fn) const {return impl.cross[fn];} + int floatingId(int xn) const {return impl.index[xn];} + + int addId() { + if (id_handler == 0) { + int xn, fn = impl.cross.size(); + if (impl.first_free == -1) { + xn = impl.index.size(); + impl.index.push_back(fn); + } else { + xn = impl.first_free; + impl.first_free = impl.index[impl.first_free]; + impl.index[xn] = fn; + } + impl.cross.push_back(xn); + return xn; + } else { + return id_handler->addId(impl); + } + } + + void eraseId(int xn) { + if (id_handler == 0) { + int fn = impl.index[xn]; + impl.index[xn] = impl.first_free; + impl.first_free = xn; + for(int i = fn + 1; i < int(impl.cross.size()); ++i) { + impl.cross[i - 1] = impl.cross[i]; + impl.index[impl.cross[i]]--; + } + impl.cross.pop_back(); + } else { + id_handler->eraseId(impl, xn); + } + } + + void firstFloating(int& fn) const { + fn = impl.first_index; + if (fn == int(impl.cross.size())) fn = -1; + } + + void nextFloating(int& fn) const { + ++fn; + if (fn == int(impl.cross.size())) fn = -1; + } + + void firstFix(int& xn) const { + int fn; + firstFloating(fn); + xn = fn != -1 ? fixId(fn) : -1; + } + + void nextFix(int& xn) const { + int fn = floatingId(xn); + nextFloating(fn); + xn = fn != -1 ? fixId(fn) : -1; + } + + protected: + LpIdImpl impl; + IdHandler *id_handler; + }; + + class RelocateIdHandler : public LpId::IdHandler { + public: + + virtual int addId(LpIdImpl& impl) { + int xn, fn = impl.cross.size(); + if (impl.first_free == -1) { + xn = impl.index.size(); + impl.index.push_back(fn); + } else { + xn = impl.first_free; + impl.first_free = impl.index[impl.first_free]; + impl.index[xn] = fn; + } + impl.cross.push_back(xn); + return xn; + } + + virtual void eraseId(LpIdImpl& impl, int xn) { + int fn = impl.index[xn]; + impl.index[xn] = impl.first_free; + impl.first_free = xn; + impl.cross[fn] = impl.cross.back(); + impl.index[impl.cross.back()] = fn; + impl.cross.pop_back(); + } + }; + } +} + +#endif diff --git a/lemon/config.h.in b/lemon/config.h.in --- a/lemon/config.h.in +++ b/lemon/config.h.in @@ -9,3 +9,6 @@ /* Define to 1 if you have GLPK. */ #undef HAVE_GLPK + +/* Define to 1 if you have SOPLEX */ +#undef HAVE_SOPLEX \ No newline at end of file diff --git a/lemon/lp.h b/lemon/lp.h new file mode 100644 --- /dev/null +++ b/lemon/lp.h @@ -0,0 +1,90 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_LP_H +#define LEMON_LP_H + +#include + + +#ifdef HAVE_GLPK +#include +#include +#elif HAVE_CPLEX +#include +#include +#elif HAVE_SOPLEX +#include +#endif + +///\file +///\brief Defines a default LP solver +///\ingroup lp_group +namespace lemon { + +#ifdef DOXYGEN + ///The default LP solver identifier + + ///The default LP solver identifier. + ///\ingroup lp_group + /// + ///Currently, the possible values are \c GLPK or \c CPLEX +#define DEFAULT_LP SOLVER + ///The default LP solver + + ///The default LP solver. + ///\ingroup lp_group + /// + ///Currently, it is either \c LpGlpk or \c LpCplex + typedef LpGlpk Lp; + ///The default LP solver identifier string + + ///The default LP solver identifier string. + ///\ingroup lp_group + /// + ///Currently, the possible values are "GLPK" or "CPLEX" + const char default_solver_name[]="SOLVER"; + + ///The default ILP solver. + + ///The default ILP solver. + ///\ingroup lp_group + /// + ///Currently, it is either \c LpGlpk or \c LpCplex + typedef MipGlpk Mip; +#else +#ifdef HAVE_GLPK +#define DEFAULT_LP GLPK + typedef LpGlpk Lp; + typedef MipGlpk Mip; + const char default_solver_name[]="GLPK"; +#elif HAVE_CPLEX +#define DEFAULT_LP CPLEX + typedef LpCplex Lp; + typedef MipCplex Mip; + const char default_solver_name[]="CPLEX"; +#elif HAVE_SOPLEX +#define DEFAULT_LP SOPLEX + typedef LpSoplex Lp; + const char default_solver_name[]="SOPLEX"; +#endif +#endif + +} //namespace lemon + +#endif //LEMON_LP_H diff --git a/lemon/lp_base.cc b/lemon/lp_base.cc new file mode 100644 --- /dev/null +++ b/lemon/lp_base.cc @@ -0,0 +1,35 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +///\file +///\brief The implementation of the LP solver interface. + +#include +namespace lemon { + + const LpSolverBase::Value + LpSolverBase::INF = std::numeric_limits::infinity(); + const LpSolverBase::Value + LpSolverBase::NaN = std::numeric_limits::quiet_NaN(); + +// const LpSolverBase::Constr::Value +// LpSolverBase::Constr::INF = std::numeric_limits::infinity(); +// const LpSolverBase::Constr::Value +// LpSolverBase::Constr::NaN = std::numeric_limits::quiet_NaN(); + +} //namespace lemon diff --git a/lemon/lp_base.h b/lemon/lp_base.h new file mode 100644 --- /dev/null +++ b/lemon/lp_base.h @@ -0,0 +1,1705 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_LP_BASE_H +#define LEMON_LP_BASE_H + +#include +#include +#include +#include +#include + +#include +#include + +///\file +///\brief The interface of the LP solver interface. +///\ingroup lp_group +namespace lemon { + + /// Function to decide whether a floating point value is finite or not. + + /// Retruns true if the argument is not infinity, minus infinity or NaN. + /// It does the same as the isfinite() function defined by C99. + template + bool isFinite(T value) + { + typedef std::numeric_limits Lim; + if ((Lim::has_infinity && (value == Lim::infinity() || value == + -Lim::infinity())) || + ((Lim::has_quiet_NaN || Lim::has_signaling_NaN) && value != value)) + { + return false; + } + return true; + } + + ///Common base class for LP solvers + + ///\todo Much more docs + ///\ingroup lp_group + class LpSolverBase { + + protected: + + _lp_bits::LpId rows; + _lp_bits::LpId cols; + + public: + + ///Possible outcomes of an LP solving procedure + enum SolveExitStatus { + ///This means that the problem has been successfully solved: either + ///an optimal solution has been found or infeasibility/unboundedness + ///has been proved. + SOLVED = 0, + ///Any other case (including the case when some user specified + ///limit has been exceeded) + UNSOLVED = 1 + }; + + ///\e + enum SolutionStatus { + ///Feasible solution hasn't been found (but may exist). + + ///\todo NOTFOUND might be a better name. + /// + UNDEFINED = 0, + ///The problem has no feasible solution + INFEASIBLE = 1, + ///Feasible solution found + FEASIBLE = 2, + ///Optimal solution exists and found + OPTIMAL = 3, + ///The cost function is unbounded + + ///\todo Give a feasible solution and an infinite ray (and the + ///corresponding bases) + INFINITE = 4 + }; + + ///\e The type of the investigated LP problem + enum ProblemTypes { + ///Primal-dual feasible + PRIMAL_DUAL_FEASIBLE = 0, + ///Primal feasible dual infeasible + PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1, + ///Primal infeasible dual feasible + PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2, + ///Primal-dual infeasible + PRIMAL_DUAL_INFEASIBLE = 3, + ///Could not determine so far + UNKNOWN = 4 + }; + + ///The floating point type used by the solver + typedef double Value; + ///The infinity constant + static const Value INF; + ///The not a number constant + static const Value NaN; + + static inline bool isNaN(const Value& v) { return v!=v; } + + friend class Col; + friend class ColIt; + friend class Row; + + ///Refer to a column of the LP. + + ///This type is used to refer to a column of the LP. + /// + ///Its value remains valid and correct even after the addition or erase of + ///other columns. + /// + ///\todo Document what can one do with a Col (INVALID, comparing, + ///it is similar to Node/Edge) + class Col { + protected: + int id; + friend class LpSolverBase; + friend class MipSolverBase; + explicit Col(int _id) : id(_id) {} + public: + typedef Value ExprValue; + typedef True LpSolverCol; + Col() {} + Col(const Invalid&) : id(-1) {} + bool operator< (Col c) const {return id< c.id;} + bool operator> (Col c) const {return id> c.id;} + bool operator==(Col c) const {return id==c.id;} + bool operator!=(Col c) const {return id!=c.id;} + }; + + class ColIt : public Col { + const LpSolverBase *_lp; + public: + ColIt() {} + ColIt(const LpSolverBase &lp) : _lp(&lp) + { + _lp->cols.firstFix(id); + } + ColIt(const Invalid&) : Col(INVALID) {} + ColIt &operator++() + { + _lp->cols.nextFix(id); + return *this; + } + }; + + static int id(const Col& col) { return col.id; } + + + ///Refer to a row of the LP. + + ///This type is used to refer to a row of the LP. + /// + ///Its value remains valid and correct even after the addition or erase of + ///other rows. + /// + ///\todo Document what can one do with a Row (INVALID, comparing, + ///it is similar to Node/Edge) + class Row { + protected: + int id; + friend class LpSolverBase; + explicit Row(int _id) : id(_id) {} + public: + typedef Value ExprValue; + typedef True LpSolverRow; + Row() {} + Row(const Invalid&) : id(-1) {} + + bool operator< (Row c) const {return id< c.id;} + bool operator> (Row c) const {return id> c.id;} + bool operator==(Row c) const {return id==c.id;} + bool operator!=(Row c) const {return id!=c.id;} + }; + + class RowIt : public Row { + const LpSolverBase *_lp; + public: + RowIt() {} + RowIt(const LpSolverBase &lp) : _lp(&lp) + { + _lp->rows.firstFix(id); + } + RowIt(const Invalid&) : Row(INVALID) {} + RowIt &operator++() + { + _lp->rows.nextFix(id); + return *this; + } + }; + + static int id(const Row& row) { return row.id; } + + protected: + + int _lpId(const Col& c) const { + return cols.floatingId(id(c)); + } + + int _lpId(const Row& r) const { + return rows.floatingId(id(r)); + } + + Col _item(int i, Col) const { + return Col(cols.fixId(i)); + } + + Row _item(int i, Row) const { + return Row(rows.fixId(i)); + } + + + public: + + ///Linear expression of variables and a constant component + + ///This data structure stores a linear expression of the variables + ///(\ref Col "Col"s) and also has a constant component. + /// + ///There are several ways to access and modify the contents of this + ///container. + ///- Its it fully compatible with \c std::map, so for expamle + ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can + ///read and modify the coefficients like + ///these. + ///\code + ///e[v]=5; + ///e[v]+=12; + ///e.erase(v); + ///\endcode + ///or you can also iterate through its elements. + ///\code + ///double s=0; + ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i) + /// s+=i->second; + ///\endcode + ///(This code computes the sum of all coefficients). + ///- Numbers (double's) + ///and variables (\ref Col "Col"s) directly convert to an + ///\ref Expr and the usual linear operations are defined, so + ///\code + ///v+w + ///2*v-3.12*(v-w/2)+2 + ///v*2.1+(3*v+(v*12+w+6)*3)/2 + ///\endcode + ///are valid \ref Expr "Expr"essions. + ///The usual assignment operations are also defined. + ///\code + ///e=v+w; + ///e+=2*v-3.12*(v-w/2)+2; + ///e*=3.4; + ///e/=5; + ///\endcode + ///- The constant member can be set and read by \ref constComp() + ///\code + ///e.constComp()=12; + ///double c=e.constComp(); + ///\endcode + /// + ///\note \ref clear() not only sets all coefficients to 0 but also + ///clears the constant components. + /// + ///\sa Constr + /// + class Expr : public std::map + { + public: + typedef LpSolverBase::Col Key; + typedef LpSolverBase::Value Value; + + protected: + typedef std::map Base; + + Value const_comp; + public: + typedef True IsLinExpression; + ///\e + Expr() : Base(), const_comp(0) { } + ///\e + Expr(const Key &v) : const_comp(0) { + Base::insert(std::make_pair(v, 1)); + } + ///\e + Expr(const Value &v) : const_comp(v) {} + ///\e + void set(const Key &v,const Value &c) { + Base::insert(std::make_pair(v, c)); + } + ///\e + Value &constComp() { return const_comp; } + ///\e + const Value &constComp() const { return const_comp; } + + ///Removes the components with zero coefficient. + void simplify() { + for (Base::iterator i=Base::begin(); i!=Base::end();) { + Base::iterator j=i; + ++j; + if ((*i).second==0) Base::erase(i); + i=j; + } + } + + void simplify() const { + const_cast(this)->simplify(); + } + + ///Removes the coefficients closer to zero than \c tolerance. + void simplify(double &tolerance) { + for (Base::iterator i=Base::begin(); i!=Base::end();) { + Base::iterator j=i; + ++j; + if (std::fabs((*i).second)first]+=j->second; + const_comp+=e.const_comp; + return *this; + } + ///\e + Expr &operator-=(const Expr &e) { + for (Base::const_iterator j=e.begin(); j!=e.end(); ++j) + (*this)[j->first]-=j->second; + const_comp-=e.const_comp; + return *this; + } + ///\e + Expr &operator*=(const Value &c) { + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) + j->second*=c; + const_comp*=c; + return *this; + } + ///\e + Expr &operator/=(const Value &c) { + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) + j->second/=c; + const_comp/=c; + return *this; + } + + }; + + ///Linear constraint + + ///This data stucture represents a linear constraint in the LP. + ///Basically it is a linear expression with a lower or an upper bound + ///(or both). These parts of the constraint can be obtained by the member + ///functions \ref expr(), \ref lowerBound() and \ref upperBound(), + ///respectively. + ///There are two ways to construct a constraint. + ///- You can set the linear expression and the bounds directly + /// by the functions above. + ///- The operators \<=, == and \>= + /// are defined between expressions, or even between constraints whenever + /// it makes sense. Therefore if \c e and \c f are linear expressions and + /// \c s and \c t are numbers, then the followings are valid expressions + /// and thus they can be used directly e.g. in \ref addRow() whenever + /// it makes sense. + ///\code + /// e<=s + /// e<=f + /// e==f + /// s<=e<=t + /// e>=t + ///\endcode + ///\warning The validity of a constraint is checked only at run time, so + ///e.g. \ref addRow(x[1]\<=x[2]<=5) will compile, but will throw + ///an assertion. + class Constr + { + public: + typedef LpSolverBase::Expr Expr; + typedef Expr::Key Key; + typedef Expr::Value Value; + + protected: + Expr _expr; + Value _lb,_ub; + public: + ///\e + Constr() : _expr(), _lb(NaN), _ub(NaN) {} + ///\e + Constr(Value lb,const Expr &e,Value ub) : + _expr(e), _lb(lb), _ub(ub) {} + ///\e + Constr(const Expr &e,Value ub) : + _expr(e), _lb(NaN), _ub(ub) {} + ///\e + Constr(Value lb,const Expr &e) : + _expr(e), _lb(lb), _ub(NaN) {} + ///\e + Constr(const Expr &e) : + _expr(e), _lb(NaN), _ub(NaN) {} + ///\e + void clear() + { + _expr.clear(); + _lb=_ub=NaN; + } + + ///Reference to the linear expression + Expr &expr() { return _expr; } + ///Cont reference to the linear expression + const Expr &expr() const { return _expr; } + ///Reference to the lower bound. + + ///\return + ///- \ref INF "INF": the constraint is lower unbounded. + ///- \ref NaN "NaN": lower bound has not been set. + ///- finite number: the lower bound + Value &lowerBound() { return _lb; } + ///The const version of \ref lowerBound() + const Value &lowerBound() const { return _lb; } + ///Reference to the upper bound. + + ///\return + ///- \ref INF "INF": the constraint is upper unbounded. + ///- \ref NaN "NaN": upper bound has not been set. + ///- finite number: the upper bound + Value &upperBound() { return _ub; } + ///The const version of \ref upperBound() + const Value &upperBound() const { return _ub; } + ///Is the constraint lower bounded? + bool lowerBounded() const { + return isFinite(_lb); + } + ///Is the constraint upper bounded? + bool upperBounded() const { + return isFinite(_ub); + } + + }; + + ///Linear expression of rows + + ///This data structure represents a column of the matrix, + ///thas is it strores a linear expression of the dual variables + ///(\ref Row "Row"s). + /// + ///There are several ways to access and modify the contents of this + ///container. + ///- Its it fully compatible with \c std::map, so for expamle + ///if \c e is an DualExpr and \c v + ///and \c w are of type \ref Row, then you can + ///read and modify the coefficients like + ///these. + ///\code + ///e[v]=5; + ///e[v]+=12; + ///e.erase(v); + ///\endcode + ///or you can also iterate through its elements. + ///\code + ///double s=0; + ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i) + /// s+=i->second; + ///\endcode + ///(This code computes the sum of all coefficients). + ///- Numbers (double's) + ///and variables (\ref Row "Row"s) directly convert to an + ///\ref DualExpr and the usual linear operations are defined, so + ///\code + ///v+w + ///2*v-3.12*(v-w/2) + ///v*2.1+(3*v+(v*12+w)*3)/2 + ///\endcode + ///are valid \ref DualExpr "DualExpr"essions. + ///The usual assignment operations are also defined. + ///\code + ///e=v+w; + ///e+=2*v-3.12*(v-w/2); + ///e*=3.4; + ///e/=5; + ///\endcode + /// + ///\sa Expr + /// + class DualExpr : public std::map + { + public: + typedef LpSolverBase::Row Key; + typedef LpSolverBase::Value Value; + + protected: + typedef std::map Base; + + public: + typedef True IsLinExpression; + ///\e + DualExpr() : Base() { } + ///\e + DualExpr(const Key &v) { + Base::insert(std::make_pair(v, 1)); + } + ///\e + void set(const Key &v,const Value &c) { + Base::insert(std::make_pair(v, c)); + } + + ///Removes the components with zero coefficient. + void simplify() { + for (Base::iterator i=Base::begin(); i!=Base::end();) { + Base::iterator j=i; + ++j; + if ((*i).second==0) Base::erase(i); + i=j; + } + } + + void simplify() const { + const_cast(this)->simplify(); + } + + ///Removes the coefficients closer to zero than \c tolerance. + void simplify(double &tolerance) { + for (Base::iterator i=Base::begin(); i!=Base::end();) { + Base::iterator j=i; + ++j; + if (std::fabs((*i).second)first]+=j->second; + return *this; + } + ///\e + DualExpr &operator-=(const DualExpr &e) { + for (Base::const_iterator j=e.begin(); j!=e.end(); ++j) + (*this)[j->first]-=j->second; + return *this; + } + ///\e + DualExpr &operator*=(const Value &c) { + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) + j->second*=c; + return *this; + } + ///\e + DualExpr &operator/=(const Value &c) { + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) + j->second/=c; + return *this; + } + }; + + + private: + + template + class MappedOutputIterator { + public: + + typedef std::insert_iterator<_Expr> Base; + + typedef std::output_iterator_tag iterator_category; + typedef void difference_type; + typedef void value_type; + typedef void reference; + typedef void pointer; + + MappedOutputIterator(const Base& _base, const LpSolverBase& _lp) + : base(_base), lp(_lp) {} + + MappedOutputIterator& operator*() { + return *this; + } + + MappedOutputIterator& operator=(const std::pair& value) { + *base = std::make_pair(lp._item(value.first, typename _Expr::Key()), + value.second); + return *this; + } + + MappedOutputIterator& operator++() { + ++base; + return *this; + } + + MappedOutputIterator operator++(int) { + MappedOutputIterator tmp(*this); + ++base; + return tmp; + } + + bool operator==(const MappedOutputIterator& it) const { + return base == it.base; + } + + bool operator!=(const MappedOutputIterator& it) const { + return base != it.base; + } + + private: + Base base; + const LpSolverBase& lp; + }; + + template + class MappedInputIterator { + public: + + typedef typename Expr::const_iterator Base; + + typedef typename Base::iterator_category iterator_category; + typedef typename Base::difference_type difference_type; + typedef const std::pair value_type; + typedef value_type reference; + class pointer { + public: + pointer(value_type& _value) : value(_value) {} + value_type* operator->() { return &value; } + private: + value_type value; + }; + + MappedInputIterator(const Base& _base, const LpSolverBase& _lp) + : base(_base), lp(_lp) {} + + reference operator*() { + return std::make_pair(lp._lpId(base->first), base->second); + } + + pointer operator->() { + return pointer(operator*()); + } + + MappedInputIterator& operator++() { + ++base; + return *this; + } + + MappedInputIterator operator++(int) { + MappedInputIterator tmp(*this); + ++base; + return tmp; + } + + bool operator==(const MappedInputIterator& it) const { + return base == it.base; + } + + bool operator!=(const MappedInputIterator& it) const { + return base != it.base; + } + + private: + Base base; + const LpSolverBase& lp; + }; + + protected: + + /// STL compatible iterator for lp col + typedef MappedInputIterator ConstRowIterator; + /// STL compatible iterator for lp row + typedef MappedInputIterator ConstColIterator; + + /// STL compatible iterator for lp col + typedef MappedOutputIterator RowIterator; + /// STL compatible iterator for lp row + typedef MappedOutputIterator ColIterator; + + //Abstract virtual functions + virtual LpSolverBase* _newLp() = 0; + virtual LpSolverBase* _copyLp(){ + LpSolverBase* newlp = _newLp(); + + std::map ref; + for (LpSolverBase::ColIt it(*this); it != INVALID; ++it) { + Col ccol = newlp->addCol(); + ref[it] = ccol; + newlp->colName(ccol, colName(it)); + newlp->colLowerBound(ccol, colLowerBound(it)); + newlp->colUpperBound(ccol, colUpperBound(it)); + } + + for (LpSolverBase::RowIt it(*this); it != INVALID; ++it) { + Expr e = row(it), ce; + for (Expr::iterator jt = e.begin(); jt != e.end(); ++jt) { + ce[ref[jt->first]] = jt->second; + } + ce += e.constComp(); + Row r = newlp->addRow(ce); + + double lower, upper; + getRowBounds(it, lower, upper); + newlp->rowBounds(r, lower, upper); + } + + return newlp; + }; + + virtual int _addCol() = 0; + virtual int _addRow() = 0; + + virtual void _eraseCol(int col) = 0; + virtual void _eraseRow(int row) = 0; + + virtual void _getColName(int col, std::string & name) const = 0; + virtual void _setColName(int col, const std::string & name) = 0; + virtual int _colByName(const std::string& name) const = 0; + + virtual void _setRowCoeffs(int i, ConstRowIterator b, + ConstRowIterator e) = 0; + virtual void _getRowCoeffs(int i, RowIterator b) const = 0; + virtual void _setColCoeffs(int i, ConstColIterator b, + ConstColIterator e) = 0; + virtual void _getColCoeffs(int i, ColIterator b) const = 0; + virtual void _setCoeff(int row, int col, Value value) = 0; + virtual Value _getCoeff(int row, int col) const = 0; + virtual void _setColLowerBound(int i, Value value) = 0; + virtual Value _getColLowerBound(int i) const = 0; + virtual void _setColUpperBound(int i, Value value) = 0; + virtual Value _getColUpperBound(int i) const = 0; + virtual void _setRowBounds(int i, Value lower, Value upper) = 0; + virtual void _getRowBounds(int i, Value &lower, Value &upper) const = 0; + + virtual void _setObjCoeff(int i, Value obj_coef) = 0; + virtual Value _getObjCoeff(int i) const = 0; + virtual void _clearObj()=0; + + virtual SolveExitStatus _solve() = 0; + virtual Value _getPrimal(int i) const = 0; + virtual Value _getDual(int i) const = 0; + virtual Value _getPrimalValue() const = 0; + virtual bool _isBasicCol(int i) const = 0; + virtual SolutionStatus _getPrimalStatus() const = 0; + virtual SolutionStatus _getDualStatus() const = 0; + virtual ProblemTypes _getProblemType() const = 0; + + virtual void _setMax() = 0; + virtual void _setMin() = 0; + + + virtual bool _isMax() const = 0; + + //Own protected stuff + + //Constant component of the objective function + Value obj_const_comp; + + public: + + ///\e + LpSolverBase() : obj_const_comp(0) {} + + ///\e + virtual ~LpSolverBase() {} + + ///Creates a new LP problem + LpSolverBase* newLp() {return _newLp();} + ///Makes a copy of the LP problem + LpSolverBase* copyLp() {return _copyLp();} + + ///\name Build up and modify the LP + + ///@{ + + ///Add a new empty column (i.e a new variable) to the LP + Col addCol() { Col c; _addCol(); c.id = cols.addId(); return c;} + + ///\brief Adds several new columns + ///(i.e a variables) at once + /// + ///This magic function takes a container as its argument + ///and fills its elements + ///with new columns (i.e. variables) + ///\param t can be + ///- a standard STL compatible iterable container with + ///\ref Col as its \c values_type + ///like + ///\code + ///std::vector + ///std::list + ///\endcode + ///- a standard STL compatible iterable container with + ///\ref Col as its \c mapped_type + ///like + ///\code + ///std::map + ///\endcode + ///- an iterable lemon \ref concepts::WriteMap "write map" like + ///\code + ///ListGraph::NodeMap + ///ListGraph::EdgeMap + ///\endcode + ///\return The number of the created column. +#ifdef DOXYGEN + template + int addColSet(T &t) { return 0;} +#else + template + typename enable_if::type + addColSet(T &t,dummy<0> = 0) { + int s=0; + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;} + return s; + } + template + typename enable_if::type + addColSet(T &t,dummy<1> = 1) { + int s=0; + for(typename T::iterator i=t.begin();i!=t.end();++i) { + i->second=addCol(); + s++; + } + return s; + } + template + typename enable_if::type + addColSet(T &t,dummy<2> = 2) { + int s=0; + for(typename T::MapIt i(t); i!=INVALID; ++i) + { + i.set(addCol()); + s++; + } + return s; + } +#endif + + ///Set a column (i.e a dual constraint) of the LP + + ///\param c is the column to be modified + ///\param e is a dual linear expression (see \ref DualExpr) + ///a better one. + void col(Col c,const DualExpr &e) { + e.simplify(); + _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this), + ConstColIterator(e.end(), *this)); + } + + ///Get a column (i.e a dual constraint) of the LP + + ///\param r is the column to get + ///\return the dual expression associated to the column + DualExpr col(Col c) const { + DualExpr e; + _getColCoeffs(_lpId(c), ColIterator(std::inserter(e, e.end()), *this)); + return e; + } + + ///Add a new column to the LP + + ///\param e is a dual linear expression (see \ref DualExpr) + ///\param obj is the corresponding component of the objective + ///function. It is 0 by default. + ///\return The created column. + Col addCol(const DualExpr &e, Value o = 0) { + Col c=addCol(); + col(c,e); + objCoeff(c,o); + return c; + } + + ///Add a new empty row (i.e a new constraint) to the LP + + ///This function adds a new empty row (i.e a new constraint) to the LP. + ///\return The created row + Row addRow() { Row r; _addRow(); r.id = rows.addId(); return r;} + + ///\brief Add several new rows + ///(i.e a constraints) at once + /// + ///This magic function takes a container as its argument + ///and fills its elements + ///with new row (i.e. variables) + ///\param t can be + ///- a standard STL compatible iterable container with + ///\ref Row as its \c values_type + ///like + ///\code + ///std::vector + ///std::list + ///\endcode + ///- a standard STL compatible iterable container with + ///\ref Row as its \c mapped_type + ///like + ///\code + ///std::map + ///\endcode + ///- an iterable lemon \ref concepts::WriteMap "write map" like + ///\code + ///ListGraph::NodeMap + ///ListGraph::EdgeMap + ///\endcode + ///\return The number of rows created. +#ifdef DOXYGEN + template + int addRowSet(T &t) { return 0;} +#else + template + typename enable_if::type + addRowSet(T &t,dummy<0> = 0) { + int s=0; + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;} + return s; + } + template + typename enable_if::type + addRowSet(T &t,dummy<1> = 1) { + int s=0; + for(typename T::iterator i=t.begin();i!=t.end();++i) { + i->second=addRow(); + s++; + } + return s; + } + template + typename enable_if::type + addRowSet(T &t,dummy<2> = 2) { + int s=0; + for(typename T::MapIt i(t); i!=INVALID; ++i) + { + i.set(addRow()); + s++; + } + return s; + } +#endif + + ///Set a row (i.e a constraint) of the LP + + ///\param r is the row to be modified + ///\param l is lower bound (-\ref INF means no bound) + ///\param e is a linear expression (see \ref Expr) + ///\param u is the upper bound (\ref INF means no bound) + ///\bug This is a temporary function. The interface will change to + ///a better one. + ///\todo Option to control whether a constraint with a single variable is + ///added or not. + void row(Row r, Value l, const Expr &e, Value u) { + e.simplify(); + _setRowCoeffs(_lpId(r), ConstRowIterator(e.begin(), *this), + ConstRowIterator(e.end(), *this)); + _setRowBounds(_lpId(r),l-e.constComp(),u-e.constComp()); + } + + ///Set a row (i.e a constraint) of the LP + + ///\param r is the row to be modified + ///\param c is a linear expression (see \ref Constr) + void row(Row r, const Constr &c) { + row(r, c.lowerBounded()?c.lowerBound():-INF, + c.expr(), c.upperBounded()?c.upperBound():INF); + } + + + ///Get a row (i.e a constraint) of the LP + + ///\param r is the row to get + ///\return the expression associated to the row + Expr row(Row r) const { + Expr e; + _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this)); + return e; + } + + ///Add a new row (i.e a new constraint) to the LP + + ///\param l is the lower bound (-\ref INF means no bound) + ///\param e is a linear expression (see \ref Expr) + ///\param u is the upper bound (\ref INF means no bound) + ///\return The created row. + ///\bug This is a temporary function. The interface will change to + ///a better one. + Row addRow(Value l,const Expr &e, Value u) { + Row r=addRow(); + row(r,l,e,u); + return r; + } + + ///Add a new row (i.e a new constraint) to the LP + + ///\param c is a linear expression (see \ref Constr) + ///\return The created row. + Row addRow(const Constr &c) { + Row r=addRow(); + row(r,c); + return r; + } + ///Erase a coloumn (i.e a variable) from the LP + + ///\param c is the coloumn to be deleted + ///\todo Please check this + void eraseCol(Col c) { + _eraseCol(_lpId(c)); + cols.eraseId(c.id); + } + ///Erase a row (i.e a constraint) from the LP + + ///\param r is the row to be deleted + ///\todo Please check this + void eraseRow(Row r) { + _eraseRow(_lpId(r)); + rows.eraseId(r.id); + } + + /// Get the name of a column + + ///\param c is the coresponding coloumn + ///\return The name of the colunm + std::string colName(Col c) const { + std::string name; + _getColName(_lpId(c), name); + return name; + } + + /// Set the name of a column + + ///\param c is the coresponding coloumn + ///\param name The name to be given + void colName(Col c, const std::string& name) { + _setColName(_lpId(c), name); + } + + /// Get the column by its name + + ///\param name The name of the column + ///\return the proper column or \c INVALID + Col colByName(const std::string& name) const { + int k = _colByName(name); + return k != -1 ? Col(cols.fixId(k)) : Col(INVALID); + } + + /// Set an element of the coefficient matrix of the LP + + ///\param r is the row of the element to be modified + ///\param c is the coloumn of the element to be modified + ///\param val is the new value of the coefficient + + void coeff(Row r, Col c, Value val) { + _setCoeff(_lpId(r),_lpId(c), val); + } + + /// Get an element of the coefficient matrix of the LP + + ///\param r is the row of the element in question + ///\param c is the coloumn of the element in question + ///\return the corresponding coefficient + + Value coeff(Row r, Col c) const { + return _getCoeff(_lpId(r),_lpId(c)); + } + + /// Set the lower bound of a column (i.e a variable) + + /// The lower bound of a variable (column) has to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or -\ref INF. + void colLowerBound(Col c, Value value) { + _setColLowerBound(_lpId(c),value); + } + + /// Get the lower bound of a column (i.e a variable) + + /// This function returns the lower bound for column (variable) \t c + /// (this might be -\ref INF as well). + ///\return The lower bound for coloumn \t c + Value colLowerBound(Col c) const { + return _getColLowerBound(_lpId(c)); + } + + ///\brief Set the lower bound of several columns + ///(i.e a variables) at once + /// + ///This magic function takes a container as its argument + ///and applies the function on all of its elements. + /// The lower bound of a variable (column) has to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or -\ref INF. +#ifdef DOXYGEN + template + void colLowerBound(T &t, Value value) { return 0;} +#else + template + typename enable_if::type + colLowerBound(T &t, Value value,dummy<0> = 0) { + for(typename T::iterator i=t.begin();i!=t.end();++i) { + colLowerBound(*i, value); + } + } + template + typename enable_if::type + colLowerBound(T &t, Value value,dummy<1> = 1) { + for(typename T::iterator i=t.begin();i!=t.end();++i) { + colLowerBound(i->second, value); + } + } + template + typename enable_if::type + colLowerBound(T &t, Value value,dummy<2> = 2) { + for(typename T::MapIt i(t); i!=INVALID; ++i){ + colLowerBound(*i, value); + } + } +#endif + + /// Set the upper bound of a column (i.e a variable) + + /// The upper bound of a variable (column) has to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or \ref INF. + void colUpperBound(Col c, Value value) { + _setColUpperBound(_lpId(c),value); + }; + + /// Get the upper bound of a column (i.e a variable) + + /// This function returns the upper bound for column (variable) \t c + /// (this might be \ref INF as well). + ///\return The upper bound for coloumn \t c + Value colUpperBound(Col c) const { + return _getColUpperBound(_lpId(c)); + } + + ///\brief Set the upper bound of several columns + ///(i.e a variables) at once + /// + ///This magic function takes a container as its argument + ///and applies the function on all of its elements. + /// The upper bound of a variable (column) has to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or \ref INF. +#ifdef DOXYGEN + template + void colUpperBound(T &t, Value value) { return 0;} +#else + template + typename enable_if::type + colUpperBound(T &t, Value value,dummy<0> = 0) { + for(typename T::iterator i=t.begin();i!=t.end();++i) { + colUpperBound(*i, value); + } + } + template + typename enable_if::type + colUpperBound(T &t, Value value,dummy<1> = 1) { + for(typename T::iterator i=t.begin();i!=t.end();++i) { + colUpperBound(i->second, value); + } + } + template + typename enable_if::type + colUpperBound(T &t, Value value,dummy<2> = 2) { + for(typename T::MapIt i(t); i!=INVALID; ++i){ + colUpperBound(*i, value); + } + } +#endif + + /// Set the lower and the upper bounds of a column (i.e a variable) + + /// The lower and the upper bounds of + /// a variable (column) have to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value, -\ref INF or \ref INF. + void colBounds(Col c, Value lower, Value upper) { + _setColLowerBound(_lpId(c),lower); + _setColUpperBound(_lpId(c),upper); + } + + ///\brief Set the lower and the upper bound of several columns + ///(i.e a variables) at once + /// + ///This magic function takes a container as its argument + ///and applies the function on all of its elements. + /// The lower and the upper bounds of + /// a variable (column) have to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value, -\ref INF or \ref INF. +#ifdef DOXYGEN + template + void colBounds(T &t, Value lower, Value upper) { return 0;} +#else + template + typename enable_if::type + colBounds(T &t, Value lower, Value upper,dummy<0> = 0) { + for(typename T::iterator i=t.begin();i!=t.end();++i) { + colBounds(*i, lower, upper); + } + } + template + typename enable_if::type + colBounds(T &t, Value lower, Value upper,dummy<1> = 1) { + for(typename T::iterator i=t.begin();i!=t.end();++i) { + colBounds(i->second, lower, upper); + } + } + template + typename enable_if::type + colBounds(T &t, Value lower, Value upper,dummy<2> = 2) { + for(typename T::MapIt i(t); i!=INVALID; ++i){ + colBounds(*i, lower, upper); + } + } +#endif + + + /// Set the lower and the upper bounds of a row (i.e a constraint) + + /// The lower and the upper bound of a constraint (row) have to be + /// given by an extended number of type Value, i.e. a finite + /// number of type Value, -\ref INF or \ref INF. There is no + /// separate function for the lower and the upper bound because + /// that would have been hard to implement for CPLEX. + void rowBounds(Row c, Value lower, Value upper) { + _setRowBounds(_lpId(c),lower, upper); + } + + /// Get the lower and the upper bounds of a row (i.e a constraint) + + /// The lower and the upper bound of + /// a constraint (row) are + /// extended numbers of type Value, i.e. finite numbers of type + /// Value, -\ref INF or \ref INF. + /// \todo There is no separate function for the + /// lower and the upper bound because we had problems with the + /// implementation of the setting functions for CPLEX: + /// check out whether this can be done for these functions. + void getRowBounds(Row c, Value &lower, Value &upper) const { + _getRowBounds(_lpId(c),lower, upper); + } + + ///Set an element of the objective function + void objCoeff(Col c, Value v) {_setObjCoeff(_lpId(c),v); }; + + ///Get an element of the objective function + Value objCoeff(Col c) const { return _getObjCoeff(_lpId(c)); }; + + ///Set the objective function + + ///\param e is a linear expression of type \ref Expr. + void obj(Expr e) { + _clearObj(); + for (Expr::iterator i=e.begin(); i!=e.end(); ++i) + objCoeff((*i).first,(*i).second); + obj_const_comp=e.constComp(); + } + + ///Get the objective function + + ///\return the objective function as a linear expression of type \ref Expr. + Expr obj() const { + Expr e; + for (ColIt it(*this); it != INVALID; ++it) { + double c = objCoeff(it); + if (c != 0.0) { + e.insert(std::make_pair(it, c)); + } + } + return e; + } + + + ///Maximize + void max() { _setMax(); } + ///Minimize + void min() { _setMin(); } + + ///Query function: is this a maximization problem? + bool isMax() const {return _isMax(); } + + ///Query function: is this a minimization problem? + bool isMin() const {return !isMax(); } + + ///@} + + + ///\name Solve the LP + + ///@{ + + ///\e Solve the LP problem at hand + /// + ///\return The result of the optimization procedure. Possible + ///values and their meanings can be found in the documentation of + ///\ref SolveExitStatus. + /// + ///\todo Which method is used to solve the problem + SolveExitStatus solve() { return _solve(); } + + ///@} + + ///\name Obtain the solution + + ///@{ + + /// The status of the primal problem (the original LP problem) + SolutionStatus primalStatus() const { + return _getPrimalStatus(); + } + + /// The status of the dual (of the original LP) problem + SolutionStatus dualStatus() const { + return _getDualStatus(); + } + + ///The type of the original LP problem + ProblemTypes problemType() const { + return _getProblemType(); + } + + ///\e + Value primal(Col c) const { return _getPrimal(_lpId(c)); } + ///\e + Value primal(const Expr& e) const { + double res = e.constComp(); + for (std::map::const_iterator it = e.begin(); + it != e.end(); ++it) { + res += _getPrimal(_lpId(it->first)) * it->second; + } + return res; + } + + ///\e + Value dual(Row r) const { return _getDual(_lpId(r)); } + ///\e + Value dual(const DualExpr& e) const { + double res = 0.0; + for (std::map::const_iterator it = e.begin(); + it != e.end(); ++it) { + res += _getPrimal(_lpId(it->first)) * it->second; + } + return res; + } + + ///\e + bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); } + + ///\e + + ///\return + ///- \ref INF or -\ref INF means either infeasibility or unboundedness + /// of the primal problem, depending on whether we minimize or maximize. + ///- \ref NaN if no primal solution is found. + ///- The (finite) objective value if an optimal solution is found. + Value primalValue() const { return _getPrimalValue()+obj_const_comp;} + ///@} + + }; + + + /// \ingroup lp_group + /// + /// \brief Common base class for MIP solvers + /// \todo Much more docs + class MipSolverBase : virtual public LpSolverBase{ + public: + + ///Possible variable (coloumn) types (e.g. real, integer, binary etc.) + enum ColTypes { + ///Continuous variable + REAL = 0, + ///Integer variable + + ///Unfortunately, cplex 7.5 somewhere writes something like + ///#define INTEGER 'I' + INT = 1 + ///\todo No support for other types yet. + }; + + ///Sets the type of the given coloumn to the given type + /// + ///Sets the type of the given coloumn to the given type. + void colType(Col c, ColTypes col_type) { + _colType(_lpId(c),col_type); + } + + ///Gives back the type of the column. + /// + ///Gives back the type of the column. + ColTypes colType(Col c) const { + return _colType(_lpId(c)); + } + + ///Sets the type of the given Col to integer or remove that property. + /// + ///Sets the type of the given Col to integer or remove that property. + void integer(Col c, bool enable) { + if (enable) + colType(c,INT); + else + colType(c,REAL); + } + + ///Gives back whether the type of the column is integer or not. + /// + ///Gives back the type of the column. + ///\return true if the column has integer type and false if not. + bool integer(Col c) const { + return (colType(c)==INT); + } + + /// The status of the MIP problem + SolutionStatus mipStatus() const { + return _getMipStatus(); + } + + protected: + + virtual ColTypes _colType(int col) const = 0; + virtual void _colType(int col, ColTypes col_type) = 0; + virtual SolutionStatus _getMipStatus() const = 0; + + }; + + ///\relates LpSolverBase::Expr + /// + inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a, + const LpSolverBase::Expr &b) + { + LpSolverBase::Expr tmp(a); + tmp+=b; + return tmp; + } + ///\e + + ///\relates LpSolverBase::Expr + /// + inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a, + const LpSolverBase::Expr &b) + { + LpSolverBase::Expr tmp(a); + tmp-=b; + return tmp; + } + ///\e + + ///\relates LpSolverBase::Expr + /// + inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a, + const LpSolverBase::Value &b) + { + LpSolverBase::Expr tmp(a); + tmp*=b; + return tmp; + } + + ///\e + + ///\relates LpSolverBase::Expr + /// + inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a, + const LpSolverBase::Expr &b) + { + LpSolverBase::Expr tmp(b); + tmp*=a; + return tmp; + } + ///\e + + ///\relates LpSolverBase::Expr + /// + inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a, + const LpSolverBase::Value &b) + { + LpSolverBase::Expr tmp(a); + tmp/=b; + return tmp; + } + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e, + const LpSolverBase::Expr &f) + { + return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0); + } + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e, + const LpSolverBase::Expr &f) + { + return LpSolverBase::Constr(e,f); + } + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e, + const LpSolverBase::Value &f) + { + return LpSolverBase::Constr(-LpSolverBase::INF,e,f); + } + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e, + const LpSolverBase::Expr &f) + { + return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0); + } + + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e, + const LpSolverBase::Expr &f) + { + return LpSolverBase::Constr(f,e); + } + + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e, + const LpSolverBase::Value &f) + { + return LpSolverBase::Constr(f,e,LpSolverBase::INF); + } + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e, + const LpSolverBase::Value &f) + { + return LpSolverBase::Constr(f,e,f); + } + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e, + const LpSolverBase::Expr &f) + { + return LpSolverBase::Constr(0,e-f,0); + } + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n, + const LpSolverBase::Constr&c) + { + LpSolverBase::Constr tmp(c); + LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint"); + tmp.lowerBound()=n; + return tmp; + } + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c, + const LpSolverBase::Value &n) + { + LpSolverBase::Constr tmp(c); + LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint"); + tmp.upperBound()=n; + return tmp; + } + + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n, + const LpSolverBase::Constr&c) + { + LpSolverBase::Constr tmp(c); + LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint"); + tmp.upperBound()=n; + return tmp; + } + ///\e + + ///\relates LpSolverBase::Constr + /// + inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c, + const LpSolverBase::Value &n) + { + LpSolverBase::Constr tmp(c); + LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint"); + tmp.lowerBound()=n; + return tmp; + } + + ///\e + + ///\relates LpSolverBase::DualExpr + /// + inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a, + const LpSolverBase::DualExpr &b) + { + LpSolverBase::DualExpr tmp(a); + tmp+=b; + return tmp; + } + ///\e + + ///\relates LpSolverBase::DualExpr + /// + inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a, + const LpSolverBase::DualExpr &b) + { + LpSolverBase::DualExpr tmp(a); + tmp-=b; + return tmp; + } + ///\e + + ///\relates LpSolverBase::DualExpr + /// + inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a, + const LpSolverBase::Value &b) + { + LpSolverBase::DualExpr tmp(a); + tmp*=b; + return tmp; + } + + ///\e + + ///\relates LpSolverBase::DualExpr + /// + inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a, + const LpSolverBase::DualExpr &b) + { + LpSolverBase::DualExpr tmp(b); + tmp*=a; + return tmp; + } + ///\e + + ///\relates LpSolverBase::DualExpr + /// + inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a, + const LpSolverBase::Value &b) + { + LpSolverBase::DualExpr tmp(a); + tmp/=b; + return tmp; + } + + +} //namespace lemon + +#endif //LEMON_LP_BASE_H diff --git a/lemon/lp_cplex.cc b/lemon/lp_cplex.cc new file mode 100644 --- /dev/null +++ b/lemon/lp_cplex.cc @@ -0,0 +1,699 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#include +#include +#include + +extern "C" { +#include +} + + +///\file +///\brief Implementation of the LEMON-CPLEX lp solver interface. +namespace lemon { + + LpCplex::LpCplex() { + // env = CPXopenCPLEXdevelop(&status); + env = CPXopenCPLEX(&status); + lp = CPXcreateprob(env, &status, "LP problem"); + } + + LpCplex::LpCplex(const LpCplex& cplex) : LpSolverBase() { + env = CPXopenCPLEX(&status); + lp = CPXcloneprob(env, cplex.lp, &status); + rows = cplex.rows; + cols = cplex.cols; + } + + LpCplex::~LpCplex() { + CPXfreeprob(env,&lp); + CPXcloseCPLEX(&env); + } + + LpSolverBase* LpCplex::_newLp() + { + //The first approach opens a new environment + return new LpCplex(); + } + + LpSolverBase* LpCplex::_copyLp() { + return new LpCplex(*this); + } + + int LpCplex::_addCol() + { + int i = CPXgetnumcols(env, lp); + Value lb[1],ub[1]; + lb[0]=-INF; + ub[0]=INF; + status = CPXnewcols(env, lp, 1, NULL, lb, ub, NULL, NULL); + return i; + } + + + int LpCplex::_addRow() + { + //We want a row that is not constrained + char sense[1]; + sense[0]='L';//<= constraint + Value rhs[1]; + rhs[0]=INF; + int i = CPXgetnumrows(env, lp); + status = CPXnewrows(env, lp, 1, rhs, sense, NULL, NULL); + return i; + } + + + void LpCplex::_eraseCol(int i) { + CPXdelcols(env, lp, i, i); + } + + void LpCplex::_eraseRow(int i) { + CPXdelrows(env, lp, i, i); + } + + void LpCplex::_getColName(int col, std::string &name) const + { + ///\bug Untested + int storespace; + CPXgetcolname(env, lp, 0, 0, 0, &storespace, col, col); + if (storespace == 0) { + name.clear(); + return; + } + + storespace *= -1; + std::vector buf(storespace); + char *names[1]; + int dontcare; + ///\bug return code unchecked for error + CPXgetcolname(env, lp, names, &*buf.begin(), storespace, + &dontcare, col, col); + name = names[0]; + } + + void LpCplex::_setColName(int col, const std::string &name) + { + ///\bug Untested + char *names[1]; + names[0] = const_cast(name.c_str()); + ///\bug return code unchecked for error + CPXchgcolname(env, lp, 1, &col, names); + } + + int LpCplex::_colByName(const std::string& name) const + { + int index; + if (CPXgetcolindex(env, lp, + const_cast(name.c_str()), &index) == 0) { + return index; + } + return -1; + } + + ///\warning Data at index 0 is ignored in the arrays. + void LpCplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e) + { + std::vector indices; + std::vector rowlist; + std::vector values; + + for(ConstRowIterator it=b; it!=e; ++it) { + indices.push_back(it->first); + values.push_back(it->second); + rowlist.push_back(i); + } + + status = CPXchgcoeflist(env, lp, values.size(), + &rowlist[0], &indices[0], &values[0]); + } + + void LpCplex::_getRowCoeffs(int i, RowIterator b) const { + int tmp1, tmp2, tmp3, length; + CPXgetrows(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i); + + length = -length; + std::vector indices(length); + std::vector values(length); + + CPXgetrows(env, lp, &tmp1, &tmp2, &indices[0], &values[0], + length, &tmp3, i, i); + + for (int i = 0; i < length; ++i) { + *b = std::make_pair(indices[i], values[i]); + ++b; + } + + /// \todo implement + } + + void LpCplex::_setColCoeffs(int i, ConstColIterator b, ConstColIterator e) + { + std::vector indices; + std::vector collist; + std::vector values; + + for(ConstColIterator it=b; it!=e; ++it) { + indices.push_back(it->first); + values.push_back(it->second); + collist.push_back(i); + } + + status = CPXchgcoeflist(env, lp, values.size(), + &indices[0], &collist[0], &values[0]); + } + + void LpCplex::_getColCoeffs(int i, ColIterator b) const { + + int tmp1, tmp2, tmp3, length; + CPXgetcols(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i); + + length = -length; + std::vector indices(length); + std::vector values(length); + + CPXgetcols(env, lp, &tmp1, &tmp2, &indices[0], &values[0], + length, &tmp3, i, i); + + for (int i = 0; i < length; ++i) { + *b = std::make_pair(indices[i], values[i]); + ++b; + } + + } + + void LpCplex::_setCoeff(int row, int col, Value value) + { + CPXchgcoef(env, lp, row, col, value); + } + + LpCplex::Value LpCplex::_getCoeff(int row, int col) const + { + LpCplex::Value value; + CPXgetcoef(env, lp, row, col, &value); + return value; + } + + void LpCplex::_setColLowerBound(int i, Value value) + { + int indices[1]; + indices[0]=i; + char lu[1]; + lu[0]='L'; + Value bd[1]; + bd[0]=value; + status = CPXchgbds(env, lp, 1, indices, lu, bd); + + } + + LpCplex::Value LpCplex::_getColLowerBound(int i) const + { + LpCplex::Value x; + CPXgetlb (env, lp, &x, i, i); + if (x <= -CPX_INFBOUND) x = -INF; + return x; + } + + void LpCplex::_setColUpperBound(int i, Value value) + { + int indices[1]; + indices[0]=i; + char lu[1]; + lu[0]='U'; + Value bd[1]; + bd[0]=value; + status = CPXchgbds(env, lp, 1, indices, lu, bd); + } + + LpCplex::Value LpCplex::_getColUpperBound(int i) const + { + LpCplex::Value x; + CPXgetub (env, lp, &x, i, i); + if (x >= CPX_INFBOUND) x = INF; + return x; + } + + //This will be easier to implement + void LpCplex::_setRowBounds(int i, Value lb, Value ub) + { + //Bad parameter + if (lb==INF || ub==-INF) { + //FIXME error + } + + int cnt=1; + int indices[1]; + indices[0]=i; + char sense[1]; + + if (lb==-INF){ + sense[0]='L'; + CPXchgsense(env, lp, cnt, indices, sense); + CPXchgcoef(env, lp, i, -1, ub); + + } + else{ + if (ub==INF){ + sense[0]='G'; + CPXchgsense(env, lp, cnt, indices, sense); + CPXchgcoef(env, lp, i, -1, lb); + } + else{ + if (lb == ub){ + sense[0]='E'; + CPXchgsense(env, lp, cnt, indices, sense); + CPXchgcoef(env, lp, i, -1, lb); + } + else{ + sense[0]='R'; + CPXchgsense(env, lp, cnt, indices, sense); + CPXchgcoef(env, lp, i, -1, lb); + CPXchgcoef(env, lp, i, -2, ub-lb); + } + } + } + } + +// void LpCplex::_setRowLowerBound(int i, Value value) +// { +// //Not implemented, obsolete +// } + +// void LpCplex::_setRowUpperBound(int i, Value value) +// { +// //Not implemented, obsolete +// // //TODO Ezt kell meg megirni +// // //type of the problem +// // char sense[1]; +// // status = CPXgetsense(env, lp, sense, i, i); +// // Value rhs[1]; +// // status = CPXgetrhs(env, lp, rhs, i, i); + +// // switch (sense[0]) { +// // case 'L'://<= constraint +// // break; +// // case 'E'://= constraint +// // break; +// // case 'G'://>= constraint +// // break; +// // case 'R'://ranged constraint +// // break; +// // default: ; +// // //FIXME error +// // } + +// // status = CPXchgcoef(env, lp, i, -2, value_rng); +// } + + void LpCplex::_getRowBounds(int i, Value &lb, Value &ub) const + { + char sense; + CPXgetsense(env, lp, &sense,i,i); + lb=-INF; + ub=INF; + switch (sense) + { + case 'L': + CPXgetcoef(env, lp, i, -1, &ub); + break; + case 'G': + CPXgetcoef(env, lp, i, -1, &lb); + break; + case 'E': + CPXgetcoef(env, lp, i, -1, &lb); + ub=lb; + break; + case 'R': + CPXgetcoef(env, lp, i, -1, &lb); + Value x; + CPXgetcoef(env, lp, i, -2, &x); + ub=lb+x; + break; + } + } + + void LpCplex::_setObjCoeff(int i, Value obj_coef) + { + CPXchgcoef(env, lp, -1, i, obj_coef); + } + + LpCplex::Value LpCplex::_getObjCoeff(int i) const + { + Value x; + CPXgetcoef(env, lp, -1, i, &x); + return x; + } + + void LpCplex::_clearObj() + { + for (int i=0;i< CPXgetnumcols(env, lp);++i){ + CPXchgcoef(env, lp, -1, i, 0); + } + + } + // The routine returns zero unless an error occurred during the + // optimization. Examples of errors include exhausting available + // memory (CPXERR_NO_MEMORY) or encountering invalid data in the + // CPLEX problem object (CPXERR_NO_PROBLEM). Exceeding a + // user-specified CPLEX limit, or proving the model infeasible or + // unbounded, are not considered errors. Note that a zero return + // value does not necessarily mean that a solution exists. Use query + // routines CPXsolninfo, CPXgetstat, and CPXsolution to obtain + // further information about the status of the optimization. + LpCplex::SolveExitStatus LpCplex::_solve() + { + //CPX_PARAM_LPMETHOD + status = CPXlpopt(env, lp); + //status = CPXprimopt(env, lp); +#if CPX_VERSION >= 800 + if (status) + { + return UNSOLVED; + } + else + { + switch (CPXgetstat(env, lp)) + { + case CPX_STAT_OPTIMAL: + case CPX_STAT_INFEASIBLE: + case CPX_STAT_UNBOUNDED: + return SOLVED; + default: + return UNSOLVED; + } + } +#else + if (status == 0){ + //We want to exclude some cases + switch (CPXgetstat(env, lp)){ + case CPX_OBJ_LIM: + case CPX_IT_LIM_FEAS: + case CPX_IT_LIM_INFEAS: + case CPX_TIME_LIM_FEAS: + case CPX_TIME_LIM_INFEAS: + return UNSOLVED; + default: + return SOLVED; + } + } + else{ + return UNSOLVED; + } +#endif + } + + LpCplex::Value LpCplex::_getPrimal(int i) const + { + Value x; + CPXgetx(env, lp, &x, i, i); + return x; + } + + LpCplex::Value LpCplex::_getDual(int i) const + { + Value y; + CPXgetpi(env, lp, &y, i, i); + return y; + } + + LpCplex::Value LpCplex::_getPrimalValue() const + { + Value objval; + //method = CPXgetmethod (env, lp); + //printf("CPXgetprobtype %d \n",CPXgetprobtype(env,lp)); + CPXgetobjval(env, lp, &objval); + //printf("Objective value: %g \n",objval); + return objval; + } + bool LpCplex::_isBasicCol(int i) const + { + std::vector cstat(CPXgetnumcols(env, lp)); + CPXgetbase(env, lp, &*cstat.begin(), NULL); + return (cstat[i]==CPX_BASIC); + } + +//7.5-os cplex statusai (Vigyazat: a 9.0-asei masok!) +// This table lists the statuses, returned by the CPXgetstat() +// routine, for solutions to LP problems or mixed integer problems. If +// no solution exists, the return value is zero. + +// For Simplex, Barrier +// 1 CPX_OPTIMAL +// Optimal solution found +// 2 CPX_INFEASIBLE +// Problem infeasible +// 3 CPX_UNBOUNDED +// Problem unbounded +// 4 CPX_OBJ_LIM +// Objective limit exceeded in Phase II +// 5 CPX_IT_LIM_FEAS +// Iteration limit exceeded in Phase II +// 6 CPX_IT_LIM_INFEAS +// Iteration limit exceeded in Phase I +// 7 CPX_TIME_LIM_FEAS +// Time limit exceeded in Phase II +// 8 CPX_TIME_LIM_INFEAS +// Time limit exceeded in Phase I +// 9 CPX_NUM_BEST_FEAS +// Problem non-optimal, singularities in Phase II +// 10 CPX_NUM_BEST_INFEAS +// Problem non-optimal, singularities in Phase I +// 11 CPX_OPTIMAL_INFEAS +// Optimal solution found, unscaled infeasibilities +// 12 CPX_ABORT_FEAS +// Aborted in Phase II +// 13 CPX_ABORT_INFEAS +// Aborted in Phase I +// 14 CPX_ABORT_DUAL_INFEAS +// Aborted in barrier, dual infeasible +// 15 CPX_ABORT_PRIM_INFEAS +// Aborted in barrier, primal infeasible +// 16 CPX_ABORT_PRIM_DUAL_INFEAS +// Aborted in barrier, primal and dual infeasible +// 17 CPX_ABORT_PRIM_DUAL_FEAS +// Aborted in barrier, primal and dual feasible +// 18 CPX_ABORT_CROSSOVER +// Aborted in crossover +// 19 CPX_INForUNBD +// Infeasible or unbounded +// 20 CPX_PIVOT +// User pivot used +// +// Ezeket hova tegyem: +// ??case CPX_ABORT_DUAL_INFEAS +// ??case CPX_ABORT_CROSSOVER +// ??case CPX_INForUNBD +// ??case CPX_PIVOT + +//Some more interesting stuff: + +// CPX_PARAM_LPMETHOD 1062 int LPMETHOD +// 0 Automatic +// 1 Primal Simplex +// 2 Dual Simplex +// 3 Network Simplex +// 4 Standard Barrier +// Default: 0 +// Description: Method for linear optimization. +// Determines which algorithm is used when CPXlpopt() (or "optimize" +// in the Interactive Optimizer) is called. Currently the behavior of +// the "Automatic" setting is that CPLEX simply invokes the dual +// simplex method, but this capability may be expanded in the future +// so that CPLEX chooses the method based on problem characteristics +#if CPX_VERSION < 900 + void statusSwitch(CPXENVptr env,int& stat){ + int lpmethod; + CPXgetintparam (env,CPX_PARAM_LPMETHOD,&lpmethod); + if (lpmethod==2){ + if (stat==CPX_UNBOUNDED){ + stat=CPX_INFEASIBLE; + } + else{ + if (stat==CPX_INFEASIBLE) + stat=CPX_UNBOUNDED; + } + } + } +#else + void statusSwitch(CPXENVptr,int&){} +#endif + + LpCplex::SolutionStatus LpCplex::_getPrimalStatus() const + { + //Unboundedness not treated well: the following is from cplex 9.0 doc + // About Unboundedness + + // The treatment of models that are unbounded involves a few + // subtleties. Specifically, a declaration of unboundedness means that + // ILOG CPLEX has determined that the model has an unbounded + // ray. Given any feasible solution x with objective z, a multiple of + // the unbounded ray can be added to x to give a feasible solution + // with objective z-1 (or z+1 for maximization models). Thus, if a + // feasible solution exists, then the optimal objective is + // unbounded. Note that ILOG CPLEX has not necessarily concluded that + // a feasible solution exists. Users can call the routine CPXsolninfo + // to determine whether ILOG CPLEX has also concluded that the model + // has a feasible solution. + + int stat = CPXgetstat(env, lp); +#if CPX_VERSION >= 800 + switch (stat) + { + case CPX_STAT_OPTIMAL: + return OPTIMAL; + case CPX_STAT_UNBOUNDED: + return INFINITE; + case CPX_STAT_INFEASIBLE: + return INFEASIBLE; + default: + return UNDEFINED; + } +#else + statusSwitch(env,stat); + //CPXgetstat(env, lp); + //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL); + switch (stat) { + case 0: + return UNDEFINED; //Undefined + case CPX_OPTIMAL://Optimal + return OPTIMAL; + case CPX_UNBOUNDED://Unbounded + return INFEASIBLE;//In case of dual simplex + //return INFINITE; + case CPX_INFEASIBLE://Infeasible + // case CPX_IT_LIM_INFEAS: +// case CPX_TIME_LIM_INFEAS: +// case CPX_NUM_BEST_INFEAS: +// case CPX_OPTIMAL_INFEAS: +// case CPX_ABORT_INFEAS: +// case CPX_ABORT_PRIM_INFEAS: +// case CPX_ABORT_PRIM_DUAL_INFEAS: + return INFINITE;//In case of dual simplex + //return INFEASIBLE; +// case CPX_OBJ_LIM: +// case CPX_IT_LIM_FEAS: +// case CPX_TIME_LIM_FEAS: +// case CPX_NUM_BEST_FEAS: +// case CPX_ABORT_FEAS: +// case CPX_ABORT_PRIM_DUAL_FEAS: +// return FEASIBLE; + default: + return UNDEFINED; //Everything else comes here + //FIXME error + } +#endif + } + +//9.0-as cplex verzio statusai +// CPX_STAT_ABORT_DUAL_OBJ_LIM +// CPX_STAT_ABORT_IT_LIM +// CPX_STAT_ABORT_OBJ_LIM +// CPX_STAT_ABORT_PRIM_OBJ_LIM +// CPX_STAT_ABORT_TIME_LIM +// CPX_STAT_ABORT_USER +// CPX_STAT_FEASIBLE_RELAXED +// CPX_STAT_INFEASIBLE +// CPX_STAT_INForUNBD +// CPX_STAT_NUM_BEST +// CPX_STAT_OPTIMAL +// CPX_STAT_OPTIMAL_FACE_UNBOUNDED +// CPX_STAT_OPTIMAL_INFEAS +// CPX_STAT_OPTIMAL_RELAXED +// CPX_STAT_UNBOUNDED + + LpCplex::SolutionStatus LpCplex::_getDualStatus() const + { + int stat = CPXgetstat(env, lp); +#if CPX_VERSION >= 800 + switch (stat) + { + case CPX_STAT_OPTIMAL: + return OPTIMAL; + case CPX_STAT_UNBOUNDED: + return INFEASIBLE; + default: + return UNDEFINED; + } +#else + statusSwitch(env,stat); + switch (stat) { + case 0: + return UNDEFINED; //Undefined + case CPX_OPTIMAL://Optimal + return OPTIMAL; + case CPX_UNBOUNDED: + return INFEASIBLE; + default: + return UNDEFINED; //Everything else comes here + //FIXME error + } +#endif + } + + LpCplex::ProblemTypes LpCplex::_getProblemType() const + { + int stat = CPXgetstat(env, lp); +#if CPX_VERSION >= 800 + switch (stat) + { + case CPX_STAT_OPTIMAL: + return PRIMAL_DUAL_FEASIBLE; + case CPX_STAT_UNBOUNDED: + return PRIMAL_FEASIBLE_DUAL_INFEASIBLE; + default: + return UNKNOWN; + } +#else + switch (stat) { + case CPX_OPTIMAL://Optimal + return PRIMAL_DUAL_FEASIBLE; + case CPX_UNBOUNDED: + return PRIMAL_FEASIBLE_DUAL_INFEASIBLE; +// return PRIMAL_INFEASIBLE_DUAL_FEASIBLE; +// return PRIMAL_DUAL_INFEASIBLE; + +//Seems to be that this is all we can say for sure + default: + //In all other cases + return UNKNOWN; + //FIXME error + } +#endif + } + + void LpCplex::_setMax() + { + CPXchgobjsen(env, lp, CPX_MAX); + } + void LpCplex::_setMin() + { + CPXchgobjsen(env, lp, CPX_MIN); + } + + bool LpCplex::_isMax() const + { + if (CPXgetobjsen(env, lp)==CPX_MAX) + return true; + else + return false; + } + +} //namespace lemon + diff --git a/lemon/lp_cplex.h b/lemon/lp_cplex.h new file mode 100644 --- /dev/null +++ b/lemon/lp_cplex.h @@ -0,0 +1,113 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_LP_CPLEX_H +#define LEMON_LP_CPLEX_H + +///\file +///\brief Header of the LEMON-CPLEX lp solver interface. + +#include + +struct cpxenv; +struct cpxlp; + +namespace lemon { + + + /// \brief Interface for the CPLEX solver + /// + /// This class implements an interface for the CPLEX LP solver. + class LpCplex :virtual public LpSolverBase { + + public: + + typedef LpSolverBase Parent; + + /// \e + int status; + cpxenv* env; + cpxlp* lp; + + + /// \e + LpCplex(); + /// \e + LpCplex(const LpCplex&); + /// \e + ~LpCplex(); + + protected: + virtual LpSolverBase* _newLp(); + virtual LpSolverBase* _copyLp(); + + + virtual int _addCol(); + virtual int _addRow(); + virtual void _eraseCol(int i); + virtual void _eraseRow(int i); + virtual void _getColName(int col, std::string & name) const; + virtual void _setColName(int col, const std::string & name); + virtual int _colByName(const std::string& name) const; + virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e); + virtual void _getRowCoeffs(int i, RowIterator b) const; + virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e); + virtual void _getColCoeffs(int i, ColIterator b) const; + virtual void _setCoeff(int row, int col, Value value); + virtual Value _getCoeff(int row, int col) const; + + virtual void _setColLowerBound(int i, Value value); + virtual Value _getColLowerBound(int i) const; + virtual void _setColUpperBound(int i, Value value); + virtual Value _getColUpperBound(int i) const; + +// virtual void _setRowLowerBound(int i, Value value); +// virtual void _setRowUpperBound(int i, Value value); + virtual void _setRowBounds(int i, Value lower, Value upper); + virtual void _getRowBounds(int i, Value &lb, Value &ub) const; + virtual void _setObjCoeff(int i, Value obj_coef); + virtual Value _getObjCoeff(int i) const; + virtual void _clearObj(); + + + virtual SolveExitStatus _solve(); + virtual Value _getPrimal(int i) const; + virtual Value _getDual(int i) const; + virtual Value _getPrimalValue() const; + virtual bool _isBasicCol(int i) const; + + virtual SolutionStatus _getPrimalStatus() const; + virtual SolutionStatus _getDualStatus() const; + virtual ProblemTypes _getProblemType() const; + + + virtual void _setMax(); + virtual void _setMin(); + + virtual bool _isMax() const; + + public: + + cpxenv* cplexEnv() { return env; } + cpxlp* cplexLp() { return lp; } + + }; +} //END OF NAMESPACE LEMON + +#endif //LEMON_LP_CPLEX_H + diff --git a/lemon/lp_glpk.cc b/lemon/lp_glpk.cc new file mode 100644 --- /dev/null +++ b/lemon/lp_glpk.cc @@ -0,0 +1,644 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +///\file +///\brief Implementation of the LEMON-GLPK lp solver interface. + +#include +//#include + +extern "C" { +#include +} + +#if GLP_MAJOR_VERSION > 4 || (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION > 15) +#define LEMON_glp(func) (glp_##func) +#define LEMON_lpx(func) (lpx_##func) + +#define LEMON_GLP(def) (GLP_##def) +#define LEMON_LPX(def) (LPX_##def) + +#else + +#define LEMON_glp(func) (lpx_##func) +#define LEMON_lpx(func) (lpx_##func) + +#define LEMON_GLP(def) (LPX_##def) +#define LEMON_LPX(def) (LPX_##def) + +#endif + +namespace lemon { + + LpGlpk::LpGlpk() : Parent() { + solved = false; + rows = _lp_bits::LpId(1); + cols = _lp_bits::LpId(1); + lp = LEMON_glp(create_prob)(); + LEMON_glp(create_index)(lp); + messageLevel(0); + } + + LpGlpk::LpGlpk(const LpGlpk &glp) : Parent() { + solved = false; + rows = _lp_bits::LpId(1); + cols = _lp_bits::LpId(1); + lp = LEMON_glp(create_prob)(); + LEMON_glp(create_index)(lp); + messageLevel(0); + //Coefficient matrix, row bounds + LEMON_glp(add_rows)(lp, LEMON_glp(get_num_rows)(glp.lp)); + LEMON_glp(add_cols)(lp, LEMON_glp(get_num_cols)(glp.lp)); + int len; + std::vector ind(1+LEMON_glp(get_num_cols)(glp.lp)); + std::vector val(1+LEMON_glp(get_num_cols)(glp.lp)); + for (int i=1;i<=LEMON_glp(get_num_rows)(glp.lp);++i) + { + len=LEMON_glp(get_mat_row)(glp.lp,i,&*ind.begin(),&*val.begin()); + LEMON_glp(set_mat_row)(lp, i,len,&*ind.begin(),&*val.begin()); + LEMON_glp(set_row_bnds)(lp,i, + LEMON_glp(get_row_type)(glp.lp,i), + LEMON_glp(get_row_lb)(glp.lp,i), + LEMON_glp(get_row_ub)(glp.lp,i)); + } + + //Objective function, coloumn bounds + LEMON_glp(set_obj_dir)(lp, LEMON_glp(get_obj_dir)(glp.lp)); + //Objectif function's constant term treated separately + LEMON_glp(set_obj_coef)(lp,0,LEMON_glp(get_obj_coef)(glp.lp,0)); + for (int i=1;i<=LEMON_glp(get_num_cols)(glp.lp);++i) + { + LEMON_glp(set_obj_coef)(lp,i, + LEMON_glp(get_obj_coef)(glp.lp,i)); + LEMON_glp(set_col_bnds)(lp,i, + LEMON_glp(get_col_type)(glp.lp,i), + LEMON_glp(get_col_lb)(glp.lp,i), + LEMON_glp(get_col_ub)(glp.lp,i)); + } + rows = glp.rows; + cols = glp.cols; + } + + LpGlpk::~LpGlpk() { + LEMON_glp(delete_prob)(lp); + } + + int LpGlpk::_addCol() { + int i=LEMON_glp(add_cols)(lp, 1); + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), 0.0, 0.0); + solved = false; + return i; + } + + ///\e + + + LpSolverBase* LpGlpk::_newLp() + { + LpGlpk* newlp = new LpGlpk; + return newlp; + } + + ///\e + + LpSolverBase* LpGlpk::_copyLp() + { + LpGlpk *newlp = new LpGlpk(*this); + return newlp; + } + + int LpGlpk::_addRow() { + int i=LEMON_glp(add_rows)(lp, 1); + solved = false; + return i; + } + + + void LpGlpk::_eraseCol(int i) { + int ca[2]; + ca[1]=i; + LEMON_glp(del_cols)(lp, 1, ca); + solved = false; + } + + void LpGlpk::_eraseRow(int i) { + int ra[2]; + ra[1]=i; + LEMON_glp(del_rows)(lp, 1, ra); + solved = false; + } + + void LpGlpk::_getColName(int c, std::string & name) const + { + + const char *n = LEMON_glp(get_col_name)(lp,c); + name = n?n:""; + } + + + void LpGlpk::_setColName(int c, const std::string & name) + { + LEMON_glp(set_col_name)(lp,c,const_cast(name.c_str())); + + } + + int LpGlpk::_colByName(const std::string& name) const + { + int k = LEMON_glp(find_col)(lp, const_cast(name.c_str())); + return k > 0 ? k : -1; + } + + + void LpGlpk::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e) + { + std::vector indices; + std::vector values; + + indices.push_back(0); + values.push_back(0); + + for(ConstRowIterator it=b; it!=e; ++it) { + indices.push_back(it->first); + values.push_back(it->second); + } + + LEMON_glp(set_mat_row)(lp, i, values.size() - 1, + &indices[0], &values[0]); + + solved = false; + } + + void LpGlpk::_getRowCoeffs(int ix, RowIterator b) const + { + int length = LEMON_glp(get_mat_row)(lp, ix, 0, 0); + + std::vector indices(length + 1); + std::vector values(length + 1); + + LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]); + + for (int i = 1; i <= length; ++i) { + *b = std::make_pair(indices[i], values[i]); + ++b; + } + } + + void LpGlpk::_setColCoeffs(int ix, ConstColIterator b, ConstColIterator e) { + + std::vector indices; + std::vector values; + + indices.push_back(0); + values.push_back(0); + + for(ConstColIterator it=b; it!=e; ++it) { + indices.push_back(it->first); + values.push_back(it->second); + } + + LEMON_glp(set_mat_col)(lp, ix, values.size() - 1, + &indices[0], &values[0]); + + solved = false; + } + + void LpGlpk::_getColCoeffs(int ix, ColIterator b) const + { + int length = LEMON_glp(get_mat_col)(lp, ix, 0, 0); + + std::vector indices(length + 1); + std::vector values(length + 1); + + LEMON_glp(get_mat_col)(lp, ix, &indices[0], &values[0]); + + for (int i = 1; i <= length; ++i) { + *b = std::make_pair(indices[i], values[i]); + ++b; + } + } + + void LpGlpk::_setCoeff(int ix, int jx, Value value) + { + + if (LEMON_glp(get_num_cols)(lp) < LEMON_glp(get_num_rows)(lp)) { + + int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0); + + std::vector indices(length + 2); + std::vector values(length + 2); + + LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]); + + //The following code does not suppose that the elements of the + //array indices are sorted + bool found=false; + for (int i = 1; i <= length; ++i) { + if (indices[i]==jx){ + found=true; + values[i]=value; + break; + } + } + if (!found){ + ++length; + indices[length]=jx; + values[length]=value; + } + + LEMON_glp(set_mat_row)(lp, ix, length, &indices[0], &values[0]); + + } else { + + int length=LEMON_glp(get_mat_col)(lp, jx, 0, 0); + + std::vector indices(length + 2); + std::vector values(length + 2); + + LEMON_glp(get_mat_col)(lp, jx, &indices[0], &values[0]); + + //The following code does not suppose that the elements of the + //array indices are sorted + bool found=false; + for (int i = 1; i <= length; ++i) { + if (indices[i]==ix){ + found=true; + values[i]=value; + break; + } + } + if (!found){ + ++length; + indices[length]=ix; + values[length]=value; + } + + LEMON_glp(set_mat_col)(lp, jx, length, &indices[0], &values[0]); + } + + solved = false; + } + + LpGlpk::Value LpGlpk::_getCoeff(int ix, int jx) const + { + + int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0); + + std::vector indices(length + 1); + std::vector values(length + 1); + + LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]); + + //The following code does not suppose that the elements of the + //array indices are sorted + for (int i = 1; i <= length; ++i) { + if (indices[i]==jx){ + return values[i]; + } + } + return 0; + + } + + + void LpGlpk::_setColLowerBound(int i, Value lo) + { + if (lo==INF) { + //FIXME error + } + int b=LEMON_glp(get_col_type)(lp, i); + double up=LEMON_glp(get_col_ub)(lp, i); + if (lo==-INF) { + switch (b) { + case LEMON_GLP(FR): + case LEMON_GLP(LO): + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up); + break; + case LEMON_GLP(UP): + break; + case LEMON_GLP(DB): + case LEMON_GLP(FX): + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up); + break; + default: ; + //FIXME error + } + } else { + switch (b) { + case LEMON_GLP(FR): + case LEMON_GLP(LO): + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up); + break; + case LEMON_GLP(UP): + case LEMON_GLP(DB): + case LEMON_GLP(FX): + if (lo==up) + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up); + else + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up); + break; + default: ; + //FIXME error + } + } + + solved = false; + } + + LpGlpk::Value LpGlpk::_getColLowerBound(int i) const + { + int b=LEMON_glp(get_col_type)(lp, i); + switch (b) { + case LEMON_GLP(LO): + case LEMON_GLP(DB): + case LEMON_GLP(FX): + return LEMON_glp(get_col_lb)(lp, i); + default: ; + return -INF; + } + } + + void LpGlpk::_setColUpperBound(int i, Value up) + { + if (up==-INF) { + //FIXME error + } + int b=LEMON_glp(get_col_type)(lp, i); + double lo=LEMON_glp(get_col_lb)(lp, i); + if (up==INF) { + switch (b) { + case LEMON_GLP(FR): + case LEMON_GLP(LO): + break; + case LEMON_GLP(UP): + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up); + break; + case LEMON_GLP(DB): + case LEMON_GLP(FX): + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up); + break; + default: ; + //FIXME error + } + } else { + switch (b) { + case LEMON_GLP(FR): + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up); + break; + case LEMON_GLP(UP): + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up); + break; + case LEMON_GLP(LO): + case LEMON_GLP(DB): + case LEMON_GLP(FX): + if (lo==up) + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up); + else + LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up); + break; + default: ; + //FIXME error + } + } + + solved = false; + } + + LpGlpk::Value LpGlpk::_getColUpperBound(int i) const + { + int b=LEMON_glp(get_col_type)(lp, i); + switch (b) { + case LEMON_GLP(UP): + case LEMON_GLP(DB): + case LEMON_GLP(FX): + return LEMON_glp(get_col_ub)(lp, i); + default: ; + return INF; + } + } + + void LpGlpk::_setRowBounds(int i, Value lb, Value ub) + { + //Bad parameter + if (lb==INF || ub==-INF) { + //FIXME error + } + + if (lb == -INF){ + if (ub == INF){ + LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FR), lb, ub); + } + else{ + LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(UP), lb, ub); + } + } + else{ + if (ub==INF){ + LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(LO), lb, ub); + + } + else{ + if (lb == ub){ + LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FX), lb, ub); + } + else{ + LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(DB), lb, ub); + } + } + } + + solved = false; + } + + void LpGlpk::_getRowBounds(int i, Value &lb, Value &ub) const + { + + int b=LEMON_glp(get_row_type)(lp, i); + switch (b) { + case LEMON_GLP(FR): + case LEMON_GLP(UP): + lb = -INF; + break; + default: + lb=LEMON_glp(get_row_lb)(lp, i); + } + + switch (b) { + case LEMON_GLP(FR): + case LEMON_GLP(LO): + ub = INF; + break; + default: + ub=LEMON_glp(get_row_ub)(lp, i); + } + + } + + void LpGlpk::_setObjCoeff(int i, Value obj_coef) + { + //i=0 means the constant term (shift) + LEMON_glp(set_obj_coef)(lp, i, obj_coef); + + solved = false; + } + + LpGlpk::Value LpGlpk::_getObjCoeff(int i) const { + //i=0 means the constant term (shift) + return LEMON_glp(get_obj_coef)(lp, i); + } + + void LpGlpk::_clearObj() + { + for (int i=0;i<=LEMON_glp(get_num_cols)(lp);++i){ + LEMON_glp(set_obj_coef)(lp, i, 0); + } + + solved = false; + } + + LpGlpk::SolveExitStatus LpGlpk::_solve() + { + // A way to check the problem to be solved + //LEMON_glp(write_cpxlp(lp,"naittvan.cpx"); + + LEMON_lpx(std_basis)(lp); + int i = LEMON_lpx(simplex)(lp); + + switch (i) { + case LEMON_LPX(E_OK): + solved = true; + return SOLVED; + default: + return UNSOLVED; + } + } + + LpGlpk::Value LpGlpk::_getPrimal(int i) const + { + return LEMON_glp(get_col_prim)(lp,i); + } + + LpGlpk::Value LpGlpk::_getDual(int i) const + { + return LEMON_glp(get_row_dual)(lp,i); + } + + LpGlpk::Value LpGlpk::_getPrimalValue() const + { + return LEMON_glp(get_obj_val)(lp); + } + bool LpGlpk::_isBasicCol(int i) const + { + return (LEMON_glp(get_col_stat)(lp, i)==LEMON_GLP(BS)); + } + + + LpGlpk::SolutionStatus LpGlpk::_getPrimalStatus() const + { + if (!solved) return UNDEFINED; + int stat= LEMON_lpx(get_status)(lp); + switch (stat) { + case LEMON_LPX(UNDEF)://Undefined (no solve has been run yet) + return UNDEFINED; + case LEMON_LPX(NOFEAS)://There is no feasible solution (primal, I guess) + case LEMON_LPX(INFEAS)://Infeasible + return INFEASIBLE; + case LEMON_LPX(UNBND)://Unbounded + return INFINITE; + case LEMON_LPX(FEAS)://Feasible + return FEASIBLE; + case LEMON_LPX(OPT)://Feasible + return OPTIMAL; + default: + return UNDEFINED; //to avoid gcc warning + //FIXME error + } + } + + LpGlpk::SolutionStatus LpGlpk::_getDualStatus() const + { + if (!solved) return UNDEFINED; + switch (LEMON_lpx(get_dual_stat)(lp)) { + case LEMON_LPX(D_UNDEF)://Undefined (no solve has been run yet) + return UNDEFINED; + case LEMON_LPX(D_NOFEAS)://There is no dual feasible solution +// case LEMON_LPX(D_INFEAS://Infeasible + return INFEASIBLE; + case LEMON_LPX(D_FEAS)://Feasible + switch (LEMON_lpx(get_status)(lp)) { + case LEMON_LPX(NOFEAS): + return INFINITE; + case LEMON_LPX(OPT): + return OPTIMAL; + default: + return FEASIBLE; + } + default: + return UNDEFINED; //to avoid gcc warning + //FIXME error + } + } + + LpGlpk::ProblemTypes LpGlpk::_getProblemType() const + { + if (!solved) return UNKNOWN; + //int stat= LEMON_glp(get_status(lp); + int statp= LEMON_lpx(get_prim_stat)(lp); + int statd= LEMON_lpx(get_dual_stat)(lp); + if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_FEAS)) + return PRIMAL_DUAL_FEASIBLE; + if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_NOFEAS)) + return PRIMAL_FEASIBLE_DUAL_INFEASIBLE; + if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_FEAS)) + return PRIMAL_INFEASIBLE_DUAL_FEASIBLE; + if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_NOFEAS)) + return PRIMAL_DUAL_INFEASIBLE; + //In all other cases + return UNKNOWN; + } + + void LpGlpk::_setMax() + { + solved = false; + LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MAX)); + } + + void LpGlpk::_setMin() + { + solved = false; + LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MIN)); + } + + bool LpGlpk::_isMax() const + { + return (LEMON_glp(get_obj_dir)(lp)==LEMON_GLP(MAX)); + } + + + + void LpGlpk::messageLevel(int m) + { + LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_MSGLEV), m); + } + + void LpGlpk::presolver(bool b) + { + LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_PRESOL), b); + } + + +} //END OF NAMESPACE LEMON diff --git a/lemon/lp_glpk.h b/lemon/lp_glpk.h new file mode 100644 --- /dev/null +++ b/lemon/lp_glpk.h @@ -0,0 +1,139 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_LP_GLPK_H +#define LEMON_LP_GLPK_H + +///\file +///\brief Header of the LEMON-GLPK lp solver interface. +///\ingroup lp_group + +#include + +// forward declaration +#ifndef _GLP_PROB +#define _GLP_PROB +typedef struct { double _prob; } glp_prob; +/* LP/MIP problem object */ +#endif + +namespace lemon { + + + /// \brief Interface for the GLPK LP solver + /// + /// This class implements an interface for the GLPK LP solver. + ///\ingroup lp_group + class LpGlpk : virtual public LpSolverBase { + protected: + + typedef glp_prob LPX; + glp_prob* lp; + bool solved; + + public: + + typedef LpSolverBase Parent; + + LpGlpk(); + LpGlpk(const LpGlpk &); + ~LpGlpk(); + + protected: + virtual LpSolverBase* _newLp(); + virtual LpSolverBase* _copyLp(); + + virtual int _addCol(); + virtual int _addRow(); + virtual void _eraseCol(int i); + virtual void _eraseRow(int i); + virtual void _getColName(int col, std::string & name) const; + virtual void _setColName(int col, const std::string & name); + virtual int _colByName(const std::string& name) const; + virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e); + virtual void _getRowCoeffs(int i, RowIterator b) const; + virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e); + virtual void _getColCoeffs(int i, ColIterator b) const; + virtual void _setCoeff(int row, int col, Value value); + virtual Value _getCoeff(int row, int col) const; + + virtual void _setColLowerBound(int i, Value value); + virtual Value _getColLowerBound(int i) const; + virtual void _setColUpperBound(int i, Value value); + virtual Value _getColUpperBound(int i) const; + + virtual void _setRowBounds(int i, Value lower, Value upper); + virtual void _getRowBounds(int i, Value &lb, Value &ub) const; + virtual void _setObjCoeff(int i, Value obj_coef); + virtual Value _getObjCoeff(int i) const; + virtual void _clearObj(); + + ///\e + + ///\todo It should be clarified + /// + virtual SolveExitStatus _solve(); + virtual Value _getPrimal(int i) const; + virtual Value _getDual(int i) const; + virtual Value _getPrimalValue() const; + virtual bool _isBasicCol(int i) const; + ///\e + + ///\todo It should be clarified + /// + virtual SolutionStatus _getPrimalStatus() const; + virtual SolutionStatus _getDualStatus() const; + virtual ProblemTypes _getProblemType() const; + + virtual void _setMax(); + virtual void _setMin(); + + virtual bool _isMax() const; + + public: + ///Set the verbosity of the messages + + ///Set the verbosity of the messages + /// + ///\param m is the level of the messages output by the solver routines. + ///The possible values are: + ///- 0 --- no output (default value) + ///- 1 --- error messages only + ///- 2 --- normal output + ///- 3 --- full output (includes informational messages) + void messageLevel(int m); + ///Turns on or off the presolver + + ///Turns on (\c b is \c true) or off (\c b is \c false) the presolver + /// + ///The presolver is off by default. + void presolver(bool b); + + ///Pointer to the underlying GLPK data structure. + LPX *lpx() {return lp;} + + ///Returns the constraint identifier understood by GLPK. + int lpxRow(Row r) { return _lpId(r); } + + ///Returns the variable identifier understood by GLPK. + int lpxCol(Col c) { return _lpId(c); } + }; +} //END OF NAMESPACE LEMON + +#endif //LEMON_LP_GLPK_H + diff --git a/lemon/lp_skeleton.cc b/lemon/lp_skeleton.cc new file mode 100644 --- /dev/null +++ b/lemon/lp_skeleton.cc @@ -0,0 +1,187 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#include + +///\file +///\brief A skeleton file to implement LP solver interfaces +namespace lemon { + + LpSolverBase* LpSkeleton::_newLp() + { + LpSolverBase *tmp=0; + return tmp; + } + + LpSolverBase* LpSkeleton::_copyLp() + { + LpSolverBase *tmp=0; + return tmp; + } + + int LpSkeleton::_addCol() + { + return ++col_num; + } + + int LpSkeleton::_addRow() + { + return ++row_num; + } + + void LpSkeleton::_eraseCol(int ) { + } + + void LpSkeleton::_eraseRow(int) { + } + + void LpSkeleton::_getColName(int, std::string &) const { + } + + + void LpSkeleton::_setColName(int, const std::string &) { + } + + int LpSkeleton::_colByName(const std::string&) const { return -1; } + + + void LpSkeleton::_setRowCoeffs(int, ConstRowIterator, ConstRowIterator) { + } + + void LpSkeleton::_getRowCoeffs(int, RowIterator) const { + } + + void LpSkeleton::_setColCoeffs(int, ConstColIterator, ConstColIterator) { + } + + void LpSkeleton::_getColCoeffs(int, ColIterator) const { + } + + void LpSkeleton::_setCoeff(int, int, Value ) + { + } + + LpSkeleton::Value LpSkeleton::_getCoeff(int, int) const + { + return 0; + } + + + void LpSkeleton::_setColLowerBound(int, Value) + { + } + + LpSkeleton::Value LpSkeleton::_getColLowerBound(int) const + { + return 0; + } + + void LpSkeleton::_setColUpperBound(int, Value) + { + } + + LpSkeleton::Value LpSkeleton::_getColUpperBound(int) const + { + return 0; + } + +// void LpSkeleton::_setRowLowerBound(int, Value) +// { +// } + +// void LpSkeleton::_setRowUpperBound(int, Value) +// { +// } + + void LpSkeleton::_setRowBounds(int, Value, Value) + { + } + + void LpSkeleton::_getRowBounds(int, Value&, Value&) const + { + } + + void LpSkeleton::_setObjCoeff(int, Value) + { + } + + LpSkeleton::Value LpSkeleton::_getObjCoeff(int) const + { + return 0; + } + + void LpSkeleton::_setMax() + { + } + + void LpSkeleton::_setMin() + { + } + + bool LpSkeleton::_isMax() const + { + return true; + } + + + void LpSkeleton::_clearObj() + { + } + + LpSkeleton::SolveExitStatus LpSkeleton::_solve() + { + return SOLVED; + } + + LpSkeleton::Value LpSkeleton::_getPrimal(int) const + { + return 0; + } + + LpSkeleton::Value LpSkeleton::_getDual(int) const + { + return 0; + } + + LpSkeleton::Value LpSkeleton::_getPrimalValue() const + { + return 0; + } + + LpSkeleton::SolutionStatus LpSkeleton::_getPrimalStatus() const + { + return UNDEFINED; + } + + LpSkeleton::SolutionStatus LpSkeleton::_getDualStatus() const + { + return UNDEFINED; + } + + LpSkeleton::ProblemTypes LpSkeleton::_getProblemType() const + { + return UNKNOWN; + } + + bool LpSkeleton::_isBasicCol(int) const + { + return true; + } + +} //namespace lemon + diff --git a/lemon/lp_skeleton.h b/lemon/lp_skeleton.h new file mode 100644 --- /dev/null +++ b/lemon/lp_skeleton.h @@ -0,0 +1,183 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_LP_SKELETON +#define LEMON_LP_SKELETON + +#include + +///\file +///\brief A skeleton file to implement LP solver interfaces +namespace lemon { + + ///A skeleton class to implement LP solver interfaces + class LpSkeleton :public LpSolverBase { + int col_num,row_num; + + protected: + + ///\e + virtual LpSolverBase* _newLp(); + ///\e + virtual LpSolverBase* _copyLp(); + /// \e + virtual int _addCol(); + /// \e + virtual int _addRow(); + /// \e + virtual void _eraseCol(int i); + /// \e + virtual void _eraseRow(int i); + /// \e + virtual void _getColName(int col, std::string & name) const; + /// \e + virtual void _setColName(int col, const std::string & name); + /// \e + virtual int _colByName(const std::string& name) const; + + /// \e + virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e); + /// \e + virtual void _getRowCoeffs(int i, RowIterator b) const; + /// \e + virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e); + /// \e + virtual void _getColCoeffs(int i, ColIterator b) const; + + /// Set one element of the coefficient matrix + virtual void _setCoeff(int row, int col, Value value); + + /// Get one element of the coefficient matrix + virtual Value _getCoeff(int row, int col) const; + + /// The lower bound of a variable (column) have to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or -\ref INF. + virtual void _setColLowerBound(int i, Value value); + /// \e + + /// The lower bound of a variable (column) is an + /// extended number of type Value, i.e. a finite number of type + /// Value or -\ref INF. + virtual Value _getColLowerBound(int i) const; + + /// The upper bound of a variable (column) have to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or \ref INF. + virtual void _setColUpperBound(int i, Value value); + /// \e + + /// The upper bound of a variable (column) is an + /// extended number of type Value, i.e. a finite number of type + /// Value or \ref INF. + virtual Value _getColUpperBound(int i) const; + +// /// The lower bound of a linear expression (row) have to be given by an +// /// extended number of type Value, i.e. a finite number of type +// /// Value or -\ref INF. +// virtual void _setRowLowerBound(int i, Value value); +// /// \e + +// /// The upper bound of a linear expression (row) have to be given by an +// /// extended number of type Value, i.e. a finite number of type +// /// Value or \ref INF. +// virtual void _setRowUpperBound(int i, Value value); + + /// The lower and upper bound of a linear expression (row) have to be + /// given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or +/-\ref INF. + virtual void _setRowBounds(int i, Value lb, Value ub); + /// \e + + + /// The lower and the upper bound of + /// a constraint (row) are + /// extended numbers of type Value, i.e. finite numbers of type + /// Value, -\ref INF or \ref INF. + virtual void _getRowBounds(int i, Value &lb, Value &ub) const; + /// \e + + + /// \e + virtual void _clearObj(); + /// \e + virtual void _setObjCoeff(int i, Value obj_coef); + + /// \e + virtual Value _getObjCoeff(int i) const; + + ///\e + + ///\bug Wrong interface + /// + virtual SolveExitStatus _solve(); + + ///\e + + ///\bug Wrong interface + /// + virtual Value _getPrimal(int i) const; + + ///\e + + ///\bug Wrong interface + /// + virtual Value _getDual(int i) const; + + ///\e + + ///\bug Wrong interface + /// + virtual Value _getPrimalValue() const; + + ///\e + + ///\bug Wrong interface + /// + virtual SolutionStatus _getPrimalStatus() const; + + ////e + virtual SolutionStatus _getDualStatus() const; + + + ///\e + virtual ProblemTypes _getProblemType() const; + + ///\e + virtual void _setMax(); + ///\e + virtual void _setMin(); + + ///\e + virtual bool _isMax() const; + + + + ///\e + virtual bool _isBasicCol(int i) const; + + + + public: + LpSkeleton() : LpSolverBase(), col_num(0), row_num(0) {} + }; + +} //namespace lemon + +#endif // LEMON_LP_SKELETON diff --git a/lemon/lp_soplex.cc b/lemon/lp_soplex.cc new file mode 100644 --- /dev/null +++ b/lemon/lp_soplex.cc @@ -0,0 +1,316 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#include +#include + +#include + + +///\file +///\brief Implementation of the LEMON-SOPLEX lp solver interface. +namespace lemon { + + LpSoplex::LpSoplex() : LpSolverBase() { + rows.setIdHandler(relocateIdHandler); + cols.setIdHandler(relocateIdHandler); + soplex = new soplex::SoPlex; + solved = false; + } + + LpSoplex::~LpSoplex() { + delete soplex; + } + + LpSoplex::LpSoplex(const LpSoplex& lp) : LpSolverBase() { + rows = lp.rows; + rows.setIdHandler(relocateIdHandler); + + cols = lp.cols; + cols.setIdHandler(relocateIdHandler); + + soplex = new soplex::SoPlex; + (*static_cast(soplex)) = *(lp.soplex); + + colNames = lp.colNames; + invColNames = lp.invColNames; + + primal_value = lp.primal_value; + dual_value = lp.dual_value; + + } + + LpSolverBase* LpSoplex::_newLp() { + LpSoplex* newlp = new LpSoplex(); + return newlp; + } + + LpSolverBase* LpSoplex::_copyLp() { + LpSoplex* newlp = new LpSoplex(*this); + return newlp; + } + + int LpSoplex::_addCol() { + soplex::LPCol c; + c.setLower(-soplex::infinity); + c.setUpper(soplex::infinity); + soplex->addCol(c); + + colNames.push_back(std::string()); + primal_value.push_back(0.0); + solved = false; + + return soplex->nCols() - 1; + } + + int LpSoplex::_addRow() { + soplex::LPRow r; + r.setLhs(-soplex::infinity); + r.setRhs(soplex::infinity); + soplex->addRow(r); + + dual_value.push_back(0.0); + solved = false; + + return soplex->nRows() - 1; + } + + + void LpSoplex::_eraseCol(int i) { + soplex->removeCol(i); + invColNames.erase(colNames[i]); + colNames[i] = colNames.back(); + invColNames[colNames.back()] = i; + colNames.pop_back(); + primal_value[i] = primal_value.back(); + primal_value.pop_back(); + solved = false; + } + + void LpSoplex::_eraseRow(int i) { + soplex->removeRow(i); + dual_value[i] = dual_value.back(); + dual_value.pop_back(); + solved = false; + } + + void LpSoplex::_getColName(int c, std::string &name) const { + name = colNames[c]; + } + + void LpSoplex::_setColName(int c, const std::string &name) { + invColNames.erase(colNames[c]); + colNames[c] = name; + if (!name.empty()) { + invColNames.insert(std::make_pair(name, c)); + } + } + + int LpSoplex::_colByName(const std::string& name) const { + std::map::const_iterator it = + invColNames.find(name); + if (it != invColNames.end()) { + return it->second; + } else { + return -1; + } + } + + + void LpSoplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e) { + for (int j = 0; j < soplex->nCols(); ++j) { + soplex->changeElement(i, j, 0.0); + } + for(ConstRowIterator it = b; it != e; ++it) { + soplex->changeElement(i, it->first, it->second); + } + solved = false; + } + + void LpSoplex::_getRowCoeffs(int i, RowIterator b) const { + const soplex::SVector& vec = soplex->rowVector(i); + for (int k = 0; k < vec.size(); ++k) { + *b = std::make_pair(vec.index(k), vec.value(k)); + ++b; + } + } + + void LpSoplex::_setColCoeffs(int j, ConstColIterator b, ConstColIterator e) { + for (int i = 0; i < soplex->nRows(); ++i) { + soplex->changeElement(i, j, 0.0); + } + for(ConstColIterator it = b; it != e; ++it) { + soplex->changeElement(it->first, j, it->second); + } + solved = false; + } + + void LpSoplex::_getColCoeffs(int i, ColIterator b) const { + const soplex::SVector& vec = soplex->colVector(i); + for (int k = 0; k < vec.size(); ++k) { + *b = std::make_pair(vec.index(k), vec.value(k)); + ++b; + } + } + + void LpSoplex::_setCoeff(int i, int j, Value value) { + soplex->changeElement(i, j, value); + solved = false; + } + + LpSoplex::Value LpSoplex::_getCoeff(int i, int j) const { + return soplex->rowVector(i)[j]; + } + + void LpSoplex::_setColLowerBound(int i, Value value) { + soplex->changeLower(i, value != -INF ? value : -soplex::infinity); + solved = false; + } + + LpSoplex::Value LpSoplex::_getColLowerBound(int i) const { + double value = soplex->lower(i); + return value != -soplex::infinity ? value : -INF; + } + + void LpSoplex::_setColUpperBound(int i, Value value) { + soplex->changeUpper(i, value != INF ? value : soplex::infinity); + solved = false; + } + + LpSoplex::Value LpSoplex::_getColUpperBound(int i) const { + double value = soplex->upper(i); + return value != soplex::infinity ? value : INF; + } + + void LpSoplex::_setRowBounds(int i, Value lb, Value ub) { + soplex->changeRange(i, lb != -INF ? lb : -soplex::infinity, + ub != INF ? ub : soplex::infinity); + solved = false; + } + void LpSoplex::_getRowBounds(int i, Value &lower, Value &upper) const { + lower = soplex->lhs(i); + if (lower == -soplex::infinity) lower = -INF; + upper = soplex->rhs(i); + if (upper == -soplex::infinity) upper = INF; + } + + void LpSoplex::_setObjCoeff(int i, Value obj_coef) { + soplex->changeObj(i, obj_coef); + solved = false; + } + + LpSoplex::Value LpSoplex::_getObjCoeff(int i) const { + return soplex->obj(i); + } + + void LpSoplex::_clearObj() { + for (int i = 0; i < soplex->nCols(); ++i) { + soplex->changeObj(i, 0.0); + } + solved = false; + } + + LpSoplex::SolveExitStatus LpSoplex::_solve() { + soplex::SPxSolver::Status status = soplex->solve(); + + soplex::Vector pv(primal_value.size(), &primal_value[0]); + soplex->getPrimal(pv); + + soplex::Vector dv(dual_value.size(), &dual_value[0]); + soplex->getDual(dv); + + switch (status) { + case soplex::SPxSolver::OPTIMAL: + case soplex::SPxSolver::INFEASIBLE: + case soplex::SPxSolver::UNBOUNDED: + solved = true; + return SOLVED; + default: + return UNSOLVED; + } + } + + LpSoplex::Value LpSoplex::_getPrimal(int i) const { + return primal_value[i]; + } + + LpSoplex::Value LpSoplex::_getDual(int i) const { + return dual_value[i]; + } + + LpSoplex::Value LpSoplex::_getPrimalValue() const { + return soplex->objValue(); + } + + bool LpSoplex::_isBasicCol(int i) const { + return soplex->getBasisColStatus(i) == soplex::SPxSolver::BASIC; + } + + LpSoplex::SolutionStatus LpSoplex::_getPrimalStatus() const { + if (!solved) return UNDEFINED; + switch (soplex->status()) { + case soplex::SPxSolver::OPTIMAL: + return OPTIMAL; + case soplex::SPxSolver::UNBOUNDED: + return INFINITE; + case soplex::SPxSolver::INFEASIBLE: + return INFEASIBLE; + default: + return UNDEFINED; + } + } + + LpSoplex::SolutionStatus LpSoplex::_getDualStatus() const { + if (!solved) return UNDEFINED; + switch (soplex->status()) { + case soplex::SPxSolver::OPTIMAL: + return OPTIMAL; + case soplex::SPxSolver::UNBOUNDED: + return INFEASIBLE; + default: + return UNDEFINED; + } + } + + LpSoplex::ProblemTypes LpSoplex::_getProblemType() const { + if (!solved) return UNKNOWN; + switch (soplex->status()) { + case soplex::SPxSolver::OPTIMAL: + return PRIMAL_DUAL_FEASIBLE; + case soplex::SPxSolver::UNBOUNDED: + return PRIMAL_FEASIBLE_DUAL_INFEASIBLE; + default: + return UNKNOWN; + } + } + + void LpSoplex::_setMax() { + soplex->changeSense(soplex::SPxSolver::MAXIMIZE); + solved = false; + } + void LpSoplex::_setMin() { + soplex->changeSense(soplex::SPxSolver::MINIMIZE); + solved = false; + } + bool LpSoplex::_isMax() const { + return soplex->spxSense() == soplex::SPxSolver::MAXIMIZE; + } + + +} //namespace lemon + diff --git a/lemon/lp_soplex.h b/lemon/lp_soplex.h new file mode 100644 --- /dev/null +++ b/lemon/lp_soplex.h @@ -0,0 +1,120 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_LP_SOPLEX_H +#define LEMON_LP_SOPLEX_H + +///\file +///\brief Header of the LEMON-SOPLEX lp solver interface. + +#include +#include + +#include + +// Forward declaration +namespace soplex { + class SoPlex; +} + +namespace lemon { + + /// \ingroup lp_group + /// + /// \brief Interface for the SOPLEX solver + /// + /// This class implements an interface for the SoPlex LP solver. + /// The SoPlex library is an object oriented lp solver library + /// developed at the Konrad-Zuse-Zentrum für Informationstechnik + /// Berlin (ZIB). You can find detailed information about it at the + /// http://soplex.zib.de address. + class LpSoplex :virtual public LpSolverBase { + protected: + + _lp_bits::RelocateIdHandler relocateIdHandler; + + soplex::SoPlex* soplex; + bool solved; + + std::vector colNames; + std::map invColNames; + + std::vector primal_value; + std::vector dual_value; + + + public: + + typedef LpSolverBase Parent; + + + /// \e + LpSoplex(); + /// \e + LpSoplex(const LpSoplex&); + /// \e + ~LpSoplex(); + + protected: + + virtual LpSolverBase* _newLp(); + virtual LpSolverBase* _copyLp(); + + virtual int _addCol(); + virtual int _addRow(); + virtual void _eraseCol(int i); + virtual void _eraseRow(int i); + virtual void _getColName(int col, std::string & name) const; + virtual void _setColName(int col, const std::string & name); + virtual int _colByName(const std::string& name) const; + virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e); + virtual void _getRowCoeffs(int i, RowIterator b) const; + virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e); + virtual void _getColCoeffs(int i, ColIterator b) const; + virtual void _setCoeff(int row, int col, Value value); + virtual Value _getCoeff(int row, int col) const; + virtual void _setColLowerBound(int i, Value value); + virtual Value _getColLowerBound(int i) const; + virtual void _setColUpperBound(int i, Value value); + virtual Value _getColUpperBound(int i) const; + virtual void _setRowBounds(int i, Value lower, Value upper); + virtual void _getRowBounds(int i, Value &lower, Value &upper) const; + virtual void _setObjCoeff(int i, Value obj_coef); + virtual Value _getObjCoeff(int i) const; + virtual void _clearObj(); + + virtual SolveExitStatus _solve(); + virtual Value _getPrimal(int i) const; + virtual Value _getDual(int i) const; + virtual Value _getPrimalValue() const; + virtual bool _isBasicCol(int i) const; + + virtual SolutionStatus _getPrimalStatus() const; + virtual SolutionStatus _getDualStatus() const; + virtual ProblemTypes _getProblemType() const; + + + virtual void _setMax(); + virtual void _setMin(); + virtual bool _isMax() const; + + }; +} //END OF NAMESPACE LEMON + +#endif //LEMON_LP_SOPLEX_H + diff --git a/lemon/mip_cplex.cc b/lemon/mip_cplex.cc new file mode 100644 --- /dev/null +++ b/lemon/mip_cplex.cc @@ -0,0 +1,141 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +///\file +///\brief Implementation of the LEMON-CPLEX mip solver interface. + +#include + +extern "C" { +#include +} + +namespace lemon { + + MipCplex::MipCplex() { + //This is unnecessary: setting integrality constraints on + //variables will set this, too + + ///\todo The constant CPXPROB_MIP is + ///called CPXPROB_MILP in later versions +#if CPX_VERSION < 800 + CPXchgprobtype( env, lp, CPXPROB_MIP); +#else + CPXchgprobtype( env, lp, CPXPROB_MILP); +#endif + + } + + void MipCplex::_colType(int i, MipCplex::ColTypes col_type){ + + // Note If a variable is to be changed to binary, a call to CPXchgbds + // should also be made to change the bounds to 0 and 1. + + int indices[1]; + indices[0]=i; + char ctype[1]; + switch (col_type){ + case INT: + ctype[0]=CPX_INTEGER;//'I' + break; + case REAL: + ctype[0]=CPX_CONTINUOUS ;//'C' + break; + default:; + //FIXME problem + } + CPXchgctype (env, lp, 1, indices, ctype); + } + + MipCplex::ColTypes MipCplex::_colType(int i) const { + + char ctype[1]; + CPXgetctype (env, lp, ctype, i, i); + switch (ctype[0]){ + + case CPX_INTEGER: + return INT; + case CPX_CONTINUOUS: + return REAL; + default: + return REAL;//Error! + } + + } + + LpCplex::SolveExitStatus MipCplex::_solve(){ + + status = CPXmipopt (env, lp); + if (status==0) + return SOLVED; + else + return UNSOLVED; + + } + + + LpCplex::SolutionStatus MipCplex::_getMipStatus() const { + + int stat = CPXgetstat(env, lp); + + //Fortunately, MIP statuses did not change for cplex 8.0 + switch (stat) + { + case CPXMIP_OPTIMAL: + // Optimal integer solution has been found. + case CPXMIP_OPTIMAL_TOL: + // Optimal soluton with the tolerance defined by epgap or epagap has + // been found. + return OPTIMAL; + //This also exists in later issues + // case CPXMIP_UNBOUNDED: + //return INFINITE; + case CPXMIP_INFEASIBLE: + return INFEASIBLE; + default: + return UNDEFINED; + } + //Unboundedness not treated well: the following is from cplex 9.0 doc + // About Unboundedness + + // The treatment of models that are unbounded involves a few + // subtleties. Specifically, a declaration of unboundedness means that + // ILOG CPLEX has determined that the model has an unbounded + // ray. Given any feasible solution x with objective z, a multiple of + // the unbounded ray can be added to x to give a feasible solution + // with objective z-1 (or z+1 for maximization models). Thus, if a + // feasible solution exists, then the optimal objective is + // unbounded. Note that ILOG CPLEX has not necessarily concluded that + // a feasible solution exists. Users can call the routine CPXsolninfo + // to determine whether ILOG CPLEX has also concluded that the model + // has a feasible solution. + + } + + MipCplex::Value MipCplex::_getPrimal(int i) const { + Value x; + CPXgetmipx(env, lp, &x, i, i); + return x; + } + + MipCplex::Value MipCplex::_getPrimalValue() const { + Value objval; + CPXgetmipobjval(env, lp, &objval); + return objval; + } +} //END OF NAMESPACE LEMON diff --git a/lemon/mip_cplex.h b/lemon/mip_cplex.h new file mode 100644 --- /dev/null +++ b/lemon/mip_cplex.h @@ -0,0 +1,61 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_MIP_CPLEX_H +#define LEMON_MIP_CPLEX_H + +///\file +///\brief Header of the LEMON-CPLEX mip solver interface. +///\ingroup lp_group + + +#include + +namespace lemon { + + /// \brief Interface for the CPLEX MIP solver + /// + /// This class implements an interface for the CPLEX MIP solver. + ///\ingroup lp_group + class MipCplex : public MipSolverBase, public LpCplex{ + + public: + + typedef MipSolverBase ParentMip; + typedef LpCplex ParentLp; + + MipCplex(); + //~MipCplex(); + + + + + protected: + + virtual ColTypes _colType(int col) const; + virtual void _colType(int col, ColTypes col_type); + + virtual LpCplex::SolveExitStatus _solve(); + virtual LpCplex::SolutionStatus _getMipStatus() const; + virtual ParentLp::Value _getPrimal(int i) const; + virtual ParentLp::Value _getPrimalValue() const; + }; + +} //END OF NAMESPACE LEMON + +#endif // END OF LEMON_MIP_CPLEX_H diff --git a/lemon/mip_glpk.cc b/lemon/mip_glpk.cc new file mode 100644 --- /dev/null +++ b/lemon/mip_glpk.cc @@ -0,0 +1,154 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +///\file +///\brief Implementation of the LEMON-GLPK mip solver interface. + +#include + +extern "C" { +#include +} + +#if GLP_MAJOR_VERSION > 4 || (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION > 15) +#define LEMON_glp(func) (glp_##func) +#define LEMON_lpx(func) (lpx_##func) + +#define LEMON_GLP(def) (GLP_##def) +#define LEMON_LPX(def) (LPX_##def) + +#else + +#define LEMON_glp(func) (lpx_##func) +#define LEMON_lpx(func) (lpx_##func) + +#define LEMON_GLP(def) (LPX_##def) +#define LEMON_LPX(def) (LPX_##def) + +#endif + +namespace lemon { + + MipGlpk::MipGlpk() { +#if !(GLP_MAJOR_VERSION > 4 || \ + (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION > 15)) + LEMON_lpx(set_class)(lp,LEMON_GLP(MIP)); +#endif + } + + void MipGlpk::_colType(int i, MipGlpk::ColTypes col_type){ + switch (col_type){ + case INT: + LEMON_glp(set_col_kind)(lp,i,LEMON_GLP(IV)); + break; + case REAL: + LEMON_glp(set_col_kind)(lp,i,LEMON_GLP(CV)); + break; + default:; + //FIXME problem + } + } + + MipGlpk::ColTypes MipGlpk::_colType(int i) const { + switch (LEMON_glp(get_col_kind)(lp,i)){ + case LEMON_GLP(IV): + return INT;//Or binary + case LEMON_GLP(CV): + return REAL; + default: + return REAL;//Error! + } + + } + + LpGlpk::SolveExitStatus MipGlpk::_solve() { + int result = LEMON_lpx(simplex)(lp); + + // hack: mip does not contain integer variable +#if GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION == 16 + int tmp = -1; + if (LEMON_glp(get_num_int(lp)) == 0) { + tmp = LEMON_lpx(add_cols)(lp, 1); + LEMON_glp(set_col_bnds)(lp, tmp, LEMON_GLP(FX), 0.0, 0.0); + LEMON_glp(set_col_kind)(lp, tmp, LEMON_GLP(IV)); + } +#endif + + if (LEMON_lpx(get_status)(lp)==LEMON_LPX(OPT)) { + //Maybe we could try the routine lpx_intopt(lp), a revised + //version of lpx_integer + + result = LEMON_lpx(integer)(lp); + switch (result){ + case LEMON_LPX(E_OK): + solved = true; + break; + default: + solved = false; + } + } else { + solved = false; + } +#if GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION == 16 + if (tmp != -1) { + int tmpa[2]; + tmpa[1] = tmp; + LEMON_lpx(del_cols)(lp, 1, tmpa); + } +#endif + return solved ? SOLVED : UNSOLVED; + } + + + LpGlpk::SolutionStatus MipGlpk::_getMipStatus() const { + + if (LEMON_lpx(get_status)(lp)==LEMON_LPX(OPT)){ + //Meg kell nezni: ha az LP is infinite, akkor ez is, ha az is + //infeasible, akkor ez is, de ez lehet maskepp is infeasible. + int stat= LEMON_lpx(mip_status)(lp); + + switch (stat) { + case LEMON_LPX(I_UNDEF)://Undefined (no solve has been run yet) + return UNDEFINED; + case LEMON_LPX(I_NOFEAS)://There is no feasible integral solution + return INFEASIBLE; + // case LEMON_LPX(UNBND)://Unbounded + // return INFINITE; + case LEMON_LPX(I_FEAS)://Feasible + return FEASIBLE; + case LEMON_LPX(I_OPT)://Feasible + return OPTIMAL; + default: + return UNDEFINED; //to avoid gcc warning + //FIXME error + } + } + else + return UNDEFINED; //Maybe we could refine this: what does the LP + //relaxation look like + + } + + MipGlpk::Value MipGlpk::_getPrimal(int i) const { + return LEMON_glp(mip_col_val)(lp,i); + } + + MipGlpk::Value MipGlpk::_getPrimalValue() const { + return LEMON_glp(mip_obj_val)(lp); + } +} //END OF NAMESPACE LEMON diff --git a/lemon/mip_glpk.h b/lemon/mip_glpk.h new file mode 100644 --- /dev/null +++ b/lemon/mip_glpk.h @@ -0,0 +1,59 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_MIP_GLPK_H +#define LEMON_MIP_GLPK_H + +///\file +///\brief Header of the LEMON-GLPK mip solver interface. +///\ingroup lp_group + + +#include + +namespace lemon { + /// \brief Interface for the GLPK MIP solver + /// + /// This class implements an interface for the GLPK MIP solver. + ///\ingroup lp_group + class MipGlpk : public MipSolverBase, public LpGlpk{ + + public: + + typedef MipSolverBase ParentMip; + typedef LpGlpk ParentLp; + + MipGlpk(); + //~MipGlpk(); + + + + protected: + + virtual ColTypes _colType(int col) const; + virtual void _colType(int col, ColTypes col_type); + + virtual LpGlpk::SolveExitStatus _solve(); + virtual LpGlpk::SolutionStatus _getMipStatus() const; + virtual ParentLp::Value _getPrimal(int i) const; + virtual ParentLp::Value _getPrimalValue() const; + }; + +} //END OF NAMESPACE LEMON + +#endif // END OF LEMON_MIP_GLPK_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -18,6 +18,8 @@ hao_orlin_test heap_test kruskal_test + lp_test + mip_test maps_test max_matching_test radix_sort_test diff --git a/test/Makefile.am b/test/Makefile.am --- a/test/Makefile.am +++ b/test/Makefile.am @@ -33,6 +33,13 @@ test/time_measure_test \ test/unionfind_test +if HAVE_LP +check_PROGRAMS += test/lp_test +endif HAVE_LP +if HAVE_MIP +check_PROGRAMS += test/mip_test +endif HAVE_MIP + TESTS += $(check_PROGRAMS) XFAIL_TESTS += test/test_tools_fail$(EXEEXT) @@ -51,7 +58,9 @@ test_heap_test_SOURCES = test/heap_test.cc test_kruskal_test_SOURCES = test/kruskal_test.cc test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc +test_lp_test_SOURCES = test/lp_test.cc test_maps_test_SOURCES = test/maps_test.cc +test_mip_test_SOURCES = test/mip_test.cc test_max_matching_test_SOURCES = test/max_matching_test.cc test_path_test_SOURCES = test/path_test.cc test_preflow_test_SOURCES = test/preflow_test.cc diff --git a/test/lp_test.cc b/test/lp_test.cc new file mode 100644 --- /dev/null +++ b/test/lp_test.cc @@ -0,0 +1,423 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#include +#include +#include "test_tools.h" +#include + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifdef HAVE_GLPK +#include +#endif + +#ifdef HAVE_CPLEX +#include +#endif + +#ifdef HAVE_SOPLEX +#include +#endif + +using namespace lemon; + +void lpTest(LpSolverBase & lp) +{ + + + + typedef LpSolverBase LP; + + std::vector x(10); + // for(int i=0;i<10;i++) x.push_back(lp.addCol()); + lp.addColSet(x); + lp.colLowerBound(x,1); + lp.colUpperBound(x,1); + lp.colBounds(x,1,2); +#ifndef GYORSITAS + + std::vector y(10); + lp.addColSet(y); + + lp.colLowerBound(y,1); + lp.colUpperBound(y,1); + lp.colBounds(y,1,2); + + std::map z; + + z.insert(std::make_pair(12,INVALID)); + z.insert(std::make_pair(2,INVALID)); + z.insert(std::make_pair(7,INVALID)); + z.insert(std::make_pair(5,INVALID)); + + lp.addColSet(z); + + lp.colLowerBound(z,1); + lp.colUpperBound(z,1); + lp.colBounds(z,1,2); + + { + LP::Expr e,f,g; + LP::Col p1,p2,p3,p4,p5; + LP::Constr c; + + p1=lp.addCol(); + p2=lp.addCol(); + p3=lp.addCol(); + p4=lp.addCol(); + p5=lp.addCol(); + + e[p1]=2; + e.constComp()=12; + e[p1]+=2; + e.constComp()+=12; + e[p1]-=2; + e.constComp()-=12; + + e=2; + e=2.2; + e=p1; + e=f; + + e+=2; + e+=2.2; + e+=p1; + e+=f; + + e-=2; + e-=2.2; + e-=p1; + e-=f; + + e*=2; + e*=2.2; + e/=2; + e/=2.2; + + e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+ + (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+ + (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+ + 2.2*f+f*2.2+f/2.2+ + 2*f+f*2+f/2+ + 2.2*p1+p1*2.2+p1/2.2+ + 2*p1+p1*2+p1/2 + ); + + + c = (e <= f ); + c = (e <= 2.2); + c = (e <= 2 ); + c = (e <= p1 ); + c = (2.2<= f ); + c = (2 <= f ); + c = (p1 <= f ); + c = (p1 <= p2 ); + c = (p1 <= 2.2); + c = (p1 <= 2 ); + c = (2.2<= p2 ); + c = (2 <= p2 ); + + c = (e >= f ); + c = (e >= 2.2); + c = (e >= 2 ); + c = (e >= p1 ); + c = (2.2>= f ); + c = (2 >= f ); + c = (p1 >= f ); + c = (p1 >= p2 ); + c = (p1 >= 2.2); + c = (p1 >= 2 ); + c = (2.2>= p2 ); + c = (2 >= p2 ); + + c = (e == f ); + c = (e == 2.2); + c = (e == 2 ); + c = (e == p1 ); + c = (2.2== f ); + c = (2 == f ); + c = (p1 == f ); + //c = (p1 == p2 ); + c = (p1 == 2.2); + c = (p1 == 2 ); + c = (2.2== p2 ); + c = (2 == p2 ); + + c = (2 <= e <= 3); + c = (2 <= p1<= 3); + + c = (2 >= e >= 3); + c = (2 >= p1>= 3); + + e[x[3]]=2; + e[x[3]]=4; + e[x[3]]=1; + e.constComp()=12; + + lp.addRow(LP::INF,e,23); + lp.addRow(LP::INF,3.0*(x[1]+x[2]/2)-x[3],23); + lp.addRow(LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23); + + lp.addRow(x[1]+x[3]<=x[5]-3); + lp.addRow(-7<=x[1]+x[3]-12<=3); + lp.addRow(x[1]<=x[5]); + + std::ostringstream buf; + + + //Checking the simplify function + +// //How to check the simplify function? A map gives no information +// //on the question whether a given key is or is not stored in it, or +// //it does? +// Yes, it does, using the find() function. + e=((p1+p2)+(p1-p2)); + e.simplify(); + buf << "Coeff. of p2 should be 0"; + // std::cout<0, buf.str()); + + tolerance=0.02; + e.simplify(tolerance); + buf << "Coeff. of p2 should be 0"; + check(e.find(p2)==e.end(), buf.str()); + + + } + + { + LP::DualExpr e,f,g; + LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID, + p4 = INVALID, p5 = INVALID; + + e[p1]=2; + e[p1]+=2; + e[p1]-=2; + + e=p1; + e=f; + + e+=p1; + e+=f; + + e-=p1; + e-=f; + + e*=2; + e*=2.2; + e/=2; + e/=2.2; + + e=((p1+p2)+(p1-p2)+ + (p1+f)+(f+p1)+(f+g)+ + (p1-f)+(f-p1)+(f-g)+ + 2.2*f+f*2.2+f/2.2+ + 2*f+f*2+f/2+ + 2.2*p1+p1*2.2+p1/2.2+ + 2*p1+p1*2+p1/2 + ); + } + +#endif +} + +void solveAndCheck(LpSolverBase& lp, LpSolverBase::SolutionStatus stat, + double exp_opt) { + using std::string; + lp.solve(); + //int decimal,sign; + std::ostringstream buf; + buf << "Primalstatus should be: " << int(stat); + + // itoa(stat,buf1, 10); + check(lp.primalStatus()==stat, buf.str()); + + if (stat == LpSolverBase::OPTIMAL) { + std::ostringstream sbuf; + sbuf << "Wrong optimal value: the right optimum is " << exp_opt; + check(std::abs(lp.primalValue()-exp_opt) < 1e-3, sbuf.str()); + //+ecvt(exp_opt,2) + } +} + +void aTest(LpSolverBase & lp) +{ + typedef LpSolverBase LP; + + //The following example is very simple + + typedef LpSolverBase::Row Row; + typedef LpSolverBase::Col Col; + + + Col x1 = lp.addCol(); + Col x2 = lp.addCol(); + + + //Constraints + Row upright=lp.addRow(x1+x2 <=1); + lp.addRow(x1+x2 >=-1); + lp.addRow(x1-x2 <=1); + lp.addRow(x1-x2 >=-1); + //Nonnegativity of the variables + lp.colLowerBound(x1, 0); + lp.colLowerBound(x2, 0); + //Objective function + lp.obj(x1+x2); + + lp.max(); + + //Testing the problem retrieving routines + check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!"); + check(lp.isMax(),"This is a maximization!"); + check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!"); + // std::cout<objCoeff(x1)==1,"First term should be 1 in the obj function!"); + check(clp->isMax(),"This is a maximization!"); + check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!"); + // std::cout<colLowerBound(x1)==0, + "The lower bound for variable x1 should be 0."); + check( clp->colUpperBound(x1)==LpSolverBase::INF, + "The upper bound for variable x1 should be infty."); + + clp->getRowBounds(upright,lb,ub); + check( lb==-LpSolverBase::INF, + "The lower bound for the first row should be -infty."); + check( ub==1,"The upper bound for the first row should be 1."); + e = clp->row(upright); + check( e.size() == 2, "The row retrieval gives back wrong expression."); + check( e[x1] == 1, "The first coefficient should 1."); + check( e[x2] == 1, "The second coefficient should 1."); + + de = clp->col(x1); + check( de.size() == 4, "The col retrieval gives back wrong expression."); + check( de[upright] == 1, "The first coefficient should 1."); + + delete clp; + + //Maximization of x1+x2 + //over the triangle with vertices (0,0) (0,1) (1,0) + double expected_opt=1; + solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt); + + //Minimization + lp.min(); + expected_opt=0; + solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt); + + //Vertex (-1,0) instead of (0,0) + lp.colLowerBound(x1, -LpSolverBase::INF); + expected_opt=-1; + solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt); + + //Erase one constraint and return to maximization + lp.eraseRow(upright); + lp.max(); + expected_opt=LpSolverBase::INF; + solveAndCheck(lp, LpSolverBase::INFINITE, expected_opt); + + //Infeasibilty + lp.addRow(x1+x2 <=-2); + solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt); + + //Change problem and forget to solve + lp.min(); + check(lp.primalStatus()==LpSolverBase::UNDEFINED, + "Primalstatus should be UNDEFINED"); + + +// lp.solve(); +// if (lp.primalStatus()==LpSolverBase::OPTIMAL){ +// std::cout<< "Z = "< +#endif + +#ifdef HAVE_CPLEX +#include +#endif + +#ifdef HAVE_GLPK +#include +#endif + + +using namespace lemon; + +void solveAndCheck(MipSolverBase& lp, MipSolverBase::SolutionStatus stat, + double exp_opt) { + using std::string; + + lp.solve(); + //int decimal,sign; + std::ostringstream buf; + buf << "Primalstatus should be: " << int(stat) + <<" and it is "<