# HG changeset patch # User Balazs Dezso # Date 2008-12-02 22:48:28 # Node ID ed54c0d13df089ac0e2a36df2361f33989b98680 # Parent 7afc121e068927b0503e6c0cdcc58eb381ec50b1 Thorough redesign of the LP/MIP interface (#44) - Redesigned class structure - Redesigned iterators - Some functions in the basic interface redesigned - More complete setting functions - Ray retrieving functions - Lot of improvements - Cplex common env - CLP macro definition to config.h.in - Update lp.h to also use soplex and clp - Remove default_solver_name - New solverName() function in solvers - Handle exceptions for MipCplex test - Rename tolerance parameter to epsilon - Rename MapIt to CoeffIt - Lot of documentation improvements - Various bugfixes diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -53,6 +53,7 @@ LX_CHECK_GLPK LX_CHECK_CPLEX LX_CHECK_SOPLEX +LX_CHECK_CLP AM_CONDITIONAL([HAVE_LP], [test x"$lx_lp_found" = x"yes"]) AM_CONDITIONAL([HAVE_MIP], [test x"$lx_mip_found" = x"yes"]) @@ -120,6 +121,7 @@ echo GLPK support.................. : $lx_glpk_found echo CPLEX support................. : $lx_cplex_found echo SOPLEX support................ : $lx_soplex_found +echo CLP support................... : $lx_clp_found echo echo Build demo programs........... : $enable_demo echo Build additional tools........ : $enable_tools diff --git a/lemon/CMakeLists.txt b/lemon/CMakeLists.txt --- a/lemon/CMakeLists.txt +++ b/lemon/CMakeLists.txt @@ -4,9 +4,6 @@ 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 @@ -18,25 +18,31 @@ lemon_libemon_la_CXXFLAGS = \ $(GLPK_CFLAGS) \ $(CPLEX_CFLAGS) \ - $(SOPLEX_CXXFLAGS) + $(SOPLEX_CXXFLAGS) \ + $(CLP_CXXFLAGS) lemon_libemon_la_LDFLAGS = \ $(GLPK_LIBS) \ $(CPLEX_LIBS) \ - $(SOPLEX_LIBS) + $(SOPLEX_LIBS) \ + $(CLP_LIBS) if HAVE_GLPK -lemon_libemon_la_SOURCES += lemon/lp_glpk.cc lemon/mip_glpk.cc +lemon_libemon_la_SOURCES += lemon/lp_glpk.cc endif if HAVE_CPLEX -lemon_libemon_la_SOURCES += lemon/lp_cplex.cc lemon/mip_cplex.cc +lemon_libemon_la_SOURCES += lemon/lp_cplex.cc endif if HAVE_SOPLEX lemon_libemon_la_SOURCES += lemon/lp_soplex.cc endif +if HAVE_CLP +lemon_libemon_la_SOURCES += lemon/lp_clp.cc +endif + lemon_HEADERS += \ lemon/adaptors.h \ lemon/arg_parser.h \ @@ -65,12 +71,12 @@ lemon/list_graph.h \ lemon/lp.h \ lemon/lp_base.h \ + lemon/lp_clp.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/list_graph.h \ lemon/maps.h \ lemon/math.h \ lemon/max_matching.h \ @@ -94,9 +100,9 @@ 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/solver_bits.h \ lemon/bits/traits.h \ lemon/bits/variant.h \ lemon/bits/vector_map.h diff --git a/lemon/bits/lp_id.h b/lemon/bits/lp_id.h deleted file mode 100644 --- a/lemon/bits/lp_id.h +++ /dev/null @@ -1,157 +0,0 @@ -/* -*- 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/bits/solver_bits.h b/lemon/bits/solver_bits.h new file mode 100644 --- /dev/null +++ b/lemon/bits/solver_bits.h @@ -0,0 +1,191 @@ +/* -*- 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_SOLVER_BITS_H +#define LEMON_BITS_SOLVER_BITS_H + +namespace lemon { + + namespace _solver_bits { + + class VarIndex { + private: + struct ItemT { + int prev, next; + int index; + }; + std::vector items; + int first_item, last_item, first_free_item; + + std::vector cross; + + public: + + VarIndex() + : first_item(-1), last_item(-1), first_free_item(-1) { + } + + void clear() { + first_item = -1; + first_free_item = -1; + items.clear(); + cross.clear(); + } + + int addIndex(int idx) { + int n; + if (first_free_item == -1) { + n = items.size(); + items.push_back(ItemT()); + } else { + n = first_free_item; + first_free_item = items[n].next; + if (first_free_item != -1) { + items[first_free_item].prev = -1; + } + } + items[n].index = idx; + if (static_cast(cross.size()) <= idx) { + cross.resize(idx + 1, -1); + } + cross[idx] = n; + + items[n].prev = last_item; + items[n].next = -1; + if (last_item != -1) { + items[last_item].next = n; + } else { + first_item = n; + } + last_item = n; + + return n; + } + + int addIndex(int idx, int n) { + while (n >= static_cast(items.size())) { + items.push_back(ItemT()); + items.back().prev = -1; + items.back().next = first_free_item; + if (first_free_item != -1) { + items[first_free_item].prev = items.size() - 1; + } + first_free_item = items.size() - 1; + } + if (items[n].next != -1) { + items[items[n].next].prev = items[n].prev; + } + if (items[n].prev != -1) { + items[items[n].prev].next = items[n].next; + } else { + first_free_item = items[n].next; + } + + items[n].index = idx; + if (static_cast(cross.size()) <= idx) { + cross.resize(idx + 1, -1); + } + cross[idx] = n; + + items[n].prev = last_item; + items[n].next = -1; + if (last_item != -1) { + items[last_item].next = n; + } else { + first_item = n; + } + last_item = n; + + return n; + } + + void eraseIndex(int idx) { + int n = cross[idx]; + + if (items[n].prev != -1) { + items[items[n].prev].next = items[n].next; + } else { + first_item = items[n].next; + } + if (items[n].next != -1) { + items[items[n].next].prev = items[n].prev; + } else { + last_item = items[n].prev; + } + + if (first_free_item != -1) { + items[first_free_item].prev = n; + } + items[n].next = first_free_item; + items[n].prev = -1; + first_free_item = n; + + while (!cross.empty() && cross.back() == -1) { + cross.pop_back(); + } + } + + int maxIndex() const { + return cross.size() - 1; + } + + void shiftIndices(int idx) { + for (int i = idx + 1; i < static_cast(cross.size()); ++i) { + cross[i - 1] = cross[i]; + if (cross[i] != -1) { + --items[cross[i]].index; + } + } + cross.back() = -1; + cross.pop_back(); + while (!cross.empty() && cross.back() == -1) { + cross.pop_back(); + } + } + + void relocateIndex(int idx, int jdx) { + cross[idx] = cross[jdx]; + items[cross[jdx]].index = idx; + cross[jdx] = -1; + + while (!cross.empty() && cross.back() == -1) { + cross.pop_back(); + } + } + + int operator[](int idx) const { + return cross[idx]; + } + + int operator()(int fdx) const { + return items[fdx].index; + } + + void firstItem(int& fdx) const { + fdx = first_item; + } + + void nextItem(int& fdx) const { + fdx = items[fdx].next; + } + + }; + } +} + +#endif diff --git a/lemon/config.h.in b/lemon/config.h.in --- a/lemon/config.h.in +++ b/lemon/config.h.in @@ -11,4 +11,7 @@ #undef HAVE_GLPK /* Define to 1 if you have SOPLEX */ -#undef HAVE_SOPLEX \ No newline at end of file +#undef HAVE_SOPLEX + +/* Define to 1 if you have CLP */ +#undef HAVE_CLP diff --git a/lemon/lp.h b/lemon/lp.h --- a/lemon/lp.h +++ b/lemon/lp.h @@ -24,12 +24,12 @@ #ifdef HAVE_GLPK #include -#include #elif HAVE_CPLEX #include -#include #elif HAVE_SOPLEX #include +#elif HAVE_CLP +#include #endif ///\file @@ -43,45 +43,48 @@ ///The default LP solver identifier. ///\ingroup lp_group /// - ///Currently, the possible values are \c GLPK or \c CPLEX -#define DEFAULT_LP SOLVER + ///Currently, the possible values are \c LP_GLPK, \c LP_CPLEX, \c + ///LP_SOPLEX or \c LP_CLP +#define LEMON_DEFAULT_LP SOLVER ///The default LP solver ///The default LP solver. ///\ingroup lp_group /// - ///Currently, it is either \c LpGlpk or \c LpCplex + ///Currently, it is either \c LpGlpk, \c LpCplex, \c LpSoplex or \c LpClp typedef LpGlpk Lp; - ///The default LP solver identifier string - ///The default LP solver identifier string. + ///The default MIP solver identifier + + ///The default MIP solver identifier. ///\ingroup lp_group /// - ///Currently, the possible values are "GLPK" or "CPLEX" - const char default_solver_name[]="SOLVER"; + ///Currently, the possible values are \c MIP_GLPK or \c MIP_CPLEX +#define LEMON_DEFAULT_MIP SOLVER + ///The default MIP solver. - ///The default ILP solver. - - ///The default ILP solver. + ///The default MIP solver. ///\ingroup lp_group /// - ///Currently, it is either \c LpGlpk or \c LpCplex + ///Currently, it is either \c MipGlpk or \c MipCplex typedef MipGlpk Mip; #else #ifdef HAVE_GLPK -#define DEFAULT_LP GLPK +# define LEMON_DEFAULT_LP LP_GLPK typedef LpGlpk Lp; +# define LEMON_DEFAULT_MIP MIP_GLPK typedef MipGlpk Mip; - const char default_solver_name[]="GLPK"; #elif HAVE_CPLEX -#define DEFAULT_LP CPLEX +# define LEMON_DEFAULT_LP LP_CPLEX typedef LpCplex Lp; +# define LEMON_DEFAULT_MIP MIP_CPLEX typedef MipCplex Mip; - const char default_solver_name[]="CPLEX"; #elif HAVE_SOPLEX -#define DEFAULT_LP SOPLEX +# define DEFAULT_LP LP_SOPLEX typedef LpSoplex Lp; - const char default_solver_name[]="SOPLEX"; +#elif HAVE_CLP +# define DEFAULT_LP LP_CLP + typedef LpClp Lp; #endif #endif diff --git a/lemon/lp_base.cc b/lemon/lp_base.cc --- a/lemon/lp_base.cc +++ b/lemon/lp_base.cc @@ -22,14 +22,7 @@ #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(); + const LpBase::Value LpBase::INF = std::numeric_limits::infinity(); + const LpBase::Value LpBase::NaN = std::numeric_limits::quiet_NaN(); } //namespace lemon diff --git a/lemon/lp_base.h b/lemon/lp_base.h --- a/lemon/lp_base.h +++ b/lemon/lp_base.h @@ -25,41 +25,28 @@ #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. + ///Common base class for LP and MIP solvers - /// 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 + ///Usually this class is not used directly, please use one of the concrete + ///implementations of the solver interface. ///\ingroup lp_group - class LpSolverBase { + class LpBase { protected: - _lp_bits::LpId rows; - _lp_bits::LpId cols; + _solver_bits::VarIndex rows; + _solver_bits::VarIndex cols; public: @@ -74,38 +61,12 @@ 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 + ///Direction of the optimization + enum Sense { + /// Minimization + MIN, + /// Maximization + MAX }; ///The floating point type used by the solver @@ -115,11 +76,10 @@ ///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; + friend class RowIt; ///Refer to a column of the LP. @@ -128,43 +88,93 @@ ///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) + ///\note This class is similar to other Item types in LEMON, like + ///Node and Arc types in digraph. class Col { + friend class LpBase; protected: - int id; - friend class LpSolverBase; - friend class MipSolverBase; - explicit Col(int _id) : id(_id) {} + int _id; + explicit Col(int id) : _id(id) {} public: typedef Value ExprValue; - typedef True LpSolverCol; + typedef True LpCol; + /// Default constructor + + /// \warning The default constructor sets the Col to an + /// undefined value. 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;} + /// Invalid constructor \& conversion. + + /// This constructor initializes the Col to be invalid. + /// \sa Invalid for more details. + Col(const Invalid&) : _id(-1) {} + /// Equality operator + + /// Two \ref Col "Col"s are equal if and only if they point to + /// the same LP column or both are invalid. + bool operator==(Col c) const {return _id == c._id;} + /// Inequality operator + + /// \sa operator==(Col c) + /// + bool operator!=(Col c) const {return _id != c._id;} + /// Artificial ordering operator. + + /// To allow the use of this object in std::map or similar + /// associative container we require this. + /// + /// \note This operator only have to define some strict ordering of + /// the items; this order has nothing to do with the iteration + /// ordering of the items. + bool operator<(Col c) const {return _id < c._id;} }; + ///Iterator for iterate over the columns of an LP problem + + /// Its usage is quite simple, for example you can count the number + /// of columns in an LP \c lp: + ///\code + /// int count=0; + /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count; + ///\endcode class ColIt : public Col { - const LpSolverBase *_lp; + const LpBase *_solver; public: + /// Default constructor + + /// \warning The default constructor sets the iterator + /// to an undefined value. ColIt() {} - ColIt(const LpSolverBase &lp) : _lp(&lp) + /// Sets the iterator to the first Col + + /// Sets the iterator to the first Col. + /// + ColIt(const LpBase &solver) : _solver(&solver) { - _lp->cols.firstFix(id); + _solver->cols.firstItem(_id); } + /// Invalid constructor \& conversion + + /// Initialize the iterator to be invalid. + /// \sa Invalid for more details. ColIt(const Invalid&) : Col(INVALID) {} + /// Next column + + /// Assign the iterator to the next column. + /// ColIt &operator++() { - _lp->cols.nextFix(id); + _solver->cols.nextItem(_id); return *this; } }; - static int id(const Col& col) { return col.id; } - + /// \brief Returns the ID of the column. + static int id(const Col& col) { return col._id; } + /// \brief Returns the column with the given ID. + /// + /// \pre The argument should be a valid column ID in the LP problem. + static Col colFromId(int id) { return Col(id); } ///Refer to a row of the LP. @@ -173,61 +183,93 @@ ///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) + ///\note This class is similar to other Item types in LEMON, like + ///Node and Arc types in digraph. class Row { + friend class LpBase; protected: - int id; - friend class LpSolverBase; - explicit Row(int _id) : id(_id) {} + int _id; + explicit Row(int id) : _id(id) {} public: typedef Value ExprValue; - typedef True LpSolverRow; + typedef True LpRow; + /// Default constructor + + /// \warning The default constructor sets the Row to an + /// undefined value. Row() {} - Row(const Invalid&) : id(-1) {} + /// Invalid constructor \& conversion. + + /// This constructor initializes the Row to be invalid. + /// \sa Invalid for more details. + Row(const Invalid&) : _id(-1) {} + /// Equality operator - 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;} + /// Two \ref Row "Row"s are equal if and only if they point to + /// the same LP row or both are invalid. + bool operator==(Row r) const {return _id == r._id;} + /// Inequality operator + + /// \sa operator==(Row r) + /// + bool operator!=(Row r) const {return _id != r._id;} + /// Artificial ordering operator. + + /// To allow the use of this object in std::map or similar + /// associative container we require this. + /// + /// \note This operator only have to define some strict ordering of + /// the items; this order has nothing to do with the iteration + /// ordering of the items. + bool operator<(Row r) const {return _id < r._id;} }; + ///Iterator for iterate over the rows of an LP problem + + /// Its usage is quite simple, for example you can count the number + /// of rows in an LP \c lp: + ///\code + /// int count=0; + /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count; + ///\endcode class RowIt : public Row { - const LpSolverBase *_lp; + const LpBase *_solver; public: + /// Default constructor + + /// \warning The default constructor sets the iterator + /// to an undefined value. RowIt() {} - RowIt(const LpSolverBase &lp) : _lp(&lp) + /// Sets the iterator to the first Row + + /// Sets the iterator to the first Row. + /// + RowIt(const LpBase &solver) : _solver(&solver) { - _lp->rows.firstFix(id); + _solver->rows.firstItem(_id); } + /// Invalid constructor \& conversion + + /// Initialize the iterator to be invalid. + /// \sa Invalid for more details. RowIt(const Invalid&) : Row(INVALID) {} + /// Next row + + /// Assign the iterator to the next row. + /// RowIt &operator++() { - _lp->rows.nextFix(id); + _solver->rows.nextItem(_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)); - } - + /// \brief Returns the ID of the row. + static int id(const Row& row) { return row._id; } + /// \brief Returns the row with the given ID. + /// + /// \pre The argument should be a valid row ID in the LP problem. + static Row rowFromId(int id) { return Row(id); } public: @@ -238,10 +280,6 @@ /// ///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; @@ -250,10 +288,10 @@ ///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; + ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i) + /// s+=*i * primal(i); ///\endcode - ///(This code computes the sum of all coefficients). + ///(This code computes the primal value of the expression). ///- Numbers (double's) ///and variables (\ref Col "Col"s) directly convert to an ///\ref Expr and the usual linear operations are defined, so @@ -262,7 +300,7 @@ ///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. + ///are valid expressions. ///The usual assignment operations are also defined. ///\code ///e=v+w; @@ -270,105 +308,218 @@ ///e*=3.4; ///e/=5; ///\endcode - ///- The constant member can be set and read by \ref constComp() + ///- The constant member can be set and read by dereference + /// operator (unary *) + /// ///\code - ///e.constComp()=12; - ///double c=e.constComp(); + ///*e=12; + ///double c=*e; ///\endcode /// - ///\note \ref clear() not only sets all coefficients to 0 but also - ///clears the constant components. - /// ///\sa Constr - /// - class Expr : public std::map - { + class Expr { + friend class LpBase; public: - typedef LpSolverBase::Col Key; - typedef LpSolverBase::Value Value; + /// The key type of the expression + typedef LpBase::Col Key; + /// The value type of the expression + typedef LpBase::Value Value; protected: - typedef std::map Base; + Value const_comp; + std::map comps; - 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)); + typedef True SolverExpr; + /// Default constructor + + /// Construct an empty expression, the coefficients and + /// the constant component are initialized to zero. + Expr() : const_comp(0) {} + /// Construct an expression from a column + + /// Construct an expression, which has a term with \c c variable + /// and 1.0 coefficient. + Expr(const Col &c) : const_comp(0) { + typedef std::map::value_type pair_type; + comps.insert(pair_type(id(c), 1)); } - ///\e + /// Construct an expression from a constant + + /// Construct an expression, which's constant component is \c v. + /// 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; + /// Returns the coefficient of the column + Value operator[](const Col& c) const { + std::map::const_iterator it=comps.find(id(c)); + if (it != comps.end()) { + return it->second; + } else { + return 0; } } - - void simplify() const { - const_cast(this)->simplify(); + /// Returns the coefficient of the column + Value& operator[](const Col& c) { + return comps[id(c)]; + } + /// Sets the coefficient of the column + void set(const Col &c, const Value &v) { + if (v != 0.0) { + typedef std::map::value_type pair_type; + comps.insert(pair_type(id(c), v)); + } else { + comps.erase(id(c)); + } + } + /// Returns the constant component of the expression + Value& operator*() { return const_comp; } + /// Returns the constant component of the expression + const Value& operator*() const { return const_comp; } + /// \brief Removes the coefficients which's absolute value does + /// not exceed \c epsilon. It also sets to zero the constant + /// component, if it does not exceed epsilon in absolute value. + void simplify(Value epsilon = 0.0) { + std::map::iterator it=comps.begin(); + while (it != comps.end()) { + std::map::iterator jt=it; + ++jt; + if (std::fabs((*it).second) <= epsilon) comps.erase(it); + it=jt; + } + if (std::fabs(const_comp) <= epsilon) const_comp = 0; } - ///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)(this)->simplify(epsilon); } ///Sets all coefficients and the constant component to 0. void clear() { - Base::clear(); + comps.clear(); const_comp=0; } - ///\e + ///Compound assignment Expr &operator+=(const Expr &e) { - for (Base::const_iterator j=e.begin(); j!=e.end(); ++j) - (*this)[j->first]+=j->second; + for (std::map::const_iterator it=e.comps.begin(); + it!=e.comps.end(); ++it) + comps[it->first]+=it->second; const_comp+=e.const_comp; return *this; } - ///\e + ///Compound assignment Expr &operator-=(const Expr &e) { - for (Base::const_iterator j=e.begin(); j!=e.end(); ++j) - (*this)[j->first]-=j->second; + for (std::map::const_iterator it=e.comps.begin(); + it!=e.comps.end(); ++it) + comps[it->first]-=it->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; + ///Multiply with a constant + Expr &operator*=(const Value &v) { + for (std::map::iterator it=comps.begin(); + it!=comps.end(); ++it) + it->second*=v; + const_comp*=v; return *this; } - ///\e + ///Division with a constant Expr &operator/=(const Value &c) { - for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) - j->second/=c; + for (std::map::iterator it=comps.begin(); + it!=comps.end(); ++it) + it->second/=c; const_comp/=c; return *this; } + ///Iterator over the expression + + ///The iterator iterates over the terms of the expression. + /// + ///\code + ///double s=0; + ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i) + /// s+= *i * primal(i); + ///\endcode + class CoeffIt { + private: + + std::map::iterator _it, _end; + + public: + + /// Sets the iterator to the first term + + /// Sets the iterator to the first term of the expression. + /// + CoeffIt(Expr& e) + : _it(e.comps.begin()), _end(e.comps.end()){} + + /// Convert the iterator to the column of the term + operator Col() const { + return colFromId(_it->first); + } + + /// Returns the coefficient of the term + Value& operator*() { return _it->second; } + + /// Returns the coefficient of the term + const Value& operator*() const { return _it->second; } + /// Next term + + /// Assign the iterator to the next term. + /// + CoeffIt& operator++() { ++_it; return *this; } + + /// Equality operator + bool operator==(Invalid) const { return _it == _end; } + /// Inequality operator + bool operator!=(Invalid) const { return _it != _end; } + }; + + /// Const iterator over the expression + + ///The iterator iterates over the terms of the expression. + /// + ///\code + ///double s=0; + ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i) + /// s+=*i * primal(i); + ///\endcode + class ConstCoeffIt { + private: + + std::map::const_iterator _it, _end; + + public: + + /// Sets the iterator to the first term + + /// Sets the iterator to the first term of the expression. + /// + ConstCoeffIt(const Expr& e) + : _it(e.comps.begin()), _end(e.comps.end()){} + + /// Convert the iterator to the column of the term + operator Col() const { + return colFromId(_it->first); + } + + /// Returns the coefficient of the term + const Value& operator*() const { return _it->second; } + + /// Next term + + /// Assign the iterator to the next term. + /// + ConstCoeffIt& operator++() { ++_it; return *this; } + + /// Equality operator + bool operator==(Invalid) const { return _it == _end; } + /// Inequality operator + bool operator!=(Invalid) const { return _it != _end; } + }; + }; ///Linear constraint @@ -394,13 +545,13 @@ /// 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. + ///\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 fail an assertion. class Constr { public: - typedef LpSolverBase::Expr Expr; + typedef LpBase::Expr Expr; typedef Expr::Key Key; typedef Expr::Value Value; @@ -411,15 +562,8 @@ ///\e Constr() : _expr(), _lb(NaN), _ub(NaN) {} ///\e - Constr(Value lb,const Expr &e,Value ub) : + 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 @@ -453,11 +597,11 @@ const Value &upperBound() const { return _ub; } ///Is the constraint lower bounded? bool lowerBounded() const { - return isFinite(_lb); + return _lb != -INF && !std::isnan(_lb); } ///Is the constraint upper bounded? bool upperBounded() const { - return isFinite(_ub); + return _ub != INF && !std::isnan(_ub); } }; @@ -470,11 +614,6 @@ /// ///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; @@ -483,8 +622,8 @@ ///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; + ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i) + /// s+=*i; ///\endcode ///(This code computes the sum of all coefficients). ///- Numbers (double's) @@ -495,7 +634,7 @@ ///2*v-3.12*(v-w/2) ///v*2.1+(3*v+(v*12+w)*3)/2 ///\endcode - ///are valid \ref DualExpr "DualExpr"essions. + ///are valid \ref DualExpr dual expressions. ///The usual assignment operations are also defined. ///\code ///e=v+w; @@ -505,146 +644,237 @@ ///\endcode /// ///\sa Expr - /// - class DualExpr : public std::map - { + class DualExpr { + friend class LpBase; public: - typedef LpSolverBase::Row Key; - typedef LpSolverBase::Value Value; + /// The key type of the expression + typedef LpBase::Row Key; + /// The value type of the expression + typedef LpBase::Value Value; protected: - typedef std::map Base; + std::map comps; public: - typedef True IsLinExpression; - ///\e - DualExpr() : Base() { } - ///\e - DualExpr(const Key &v) { - Base::insert(std::make_pair(v, 1)); + typedef True SolverExpr; + /// Default constructor + + /// Construct an empty expression, the coefficients are + /// initialized to zero. + DualExpr() {} + /// Construct an expression from a row + + /// Construct an expression, which has a term with \c r dual + /// variable and 1.0 coefficient. + DualExpr(const Row &r) { + typedef std::map::value_type pair_type; + comps.insert(pair_type(id(r), 1)); } - ///\e - void set(const Key &v,const Value &c) { - Base::insert(std::make_pair(v, c)); + /// Returns the coefficient of the row + Value operator[](const Row& r) const { + std::map::const_iterator it = comps.find(id(r)); + if (it != comps.end()) { + return it->second; + } else { + return 0; + } } - - ///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; + /// Returns the coefficient of the row + Value& operator[](const Row& r) { + return comps[id(r)]; + } + /// Sets the coefficient of the row + void set(const Row &r, const Value &v) { + if (v != 0.0) { + typedef std::map::value_type pair_type; + comps.insert(pair_type(id(r), v)); + } else { + comps.erase(id(r)); + } + } + /// \brief Removes the coefficients which's absolute value does + /// not exceed \c epsilon. + void simplify(Value epsilon = 0.0) { + std::map::iterator it=comps.begin(); + while (it != comps.end()) { + std::map::iterator jt=it; + ++jt; + if (std::fabs((*it).second) <= epsilon) comps.erase(it); + it=jt; } } - 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)(this)->simplify(epsilon); } ///Sets all coefficients to 0. void clear() { - Base::clear(); + comps.clear(); + } + ///Compound assignment + DualExpr &operator+=(const DualExpr &e) { + for (std::map::const_iterator it=e.comps.begin(); + it!=e.comps.end(); ++it) + comps[it->first]+=it->second; + return *this; + } + ///Compound assignment + DualExpr &operator-=(const DualExpr &e) { + for (std::map::const_iterator it=e.comps.begin(); + it!=e.comps.end(); ++it) + comps[it->first]-=it->second; + return *this; + } + ///Multiply with a constant + DualExpr &operator*=(const Value &v) { + for (std::map::iterator it=comps.begin(); + it!=comps.end(); ++it) + it->second*=v; + return *this; + } + ///Division with a constant + DualExpr &operator/=(const Value &v) { + for (std::map::iterator it=comps.begin(); + it!=comps.end(); ++it) + it->second/=v; + 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 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; - } + ///Iterator over the expression + + ///The iterator iterates over the terms of the expression. + /// + ///\code + ///double s=0; + ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i) + /// s+= *i * dual(i); + ///\endcode + class CoeffIt { + private: + + std::map::iterator _it, _end; + + public: + + /// Sets the iterator to the first term + + /// Sets the iterator to the first term of the expression. + /// + CoeffIt(DualExpr& e) + : _it(e.comps.begin()), _end(e.comps.end()){} + + /// Convert the iterator to the row of the term + operator Row() const { + return rowFromId(_it->first); + } + + /// Returns the coefficient of the term + Value& operator*() { return _it->second; } + + /// Returns the coefficient of the term + const Value& operator*() const { return _it->second; } + + /// Next term + + /// Assign the iterator to the next term. + /// + CoeffIt& operator++() { ++_it; return *this; } + + /// Equality operator + bool operator==(Invalid) const { return _it == _end; } + /// Inequality operator + bool operator!=(Invalid) const { return _it != _end; } + }; + + ///Iterator over the expression + + ///The iterator iterates over the terms of the expression. + /// + ///\code + ///double s=0; + ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i) + /// s+= *i * dual(i); + ///\endcode + class ConstCoeffIt { + private: + + std::map::const_iterator _it, _end; + + public: + + /// Sets the iterator to the first term + + /// Sets the iterator to the first term of the expression. + /// + ConstCoeffIt(const DualExpr& e) + : _it(e.comps.begin()), _end(e.comps.end()){} + + /// Convert the iterator to the row of the term + operator Row() const { + return rowFromId(_it->first); + } + + /// Returns the coefficient of the term + const Value& operator*() const { return _it->second; } + + /// Next term + + /// Assign the iterator to the next term. + /// + ConstCoeffIt& operator++() { ++_it; return *this; } + + /// Equality operator + bool operator==(Invalid) const { return _it == _end; } + /// Inequality operator + bool operator!=(Invalid) const { return _it != _end; } + }; }; - private: + protected: - template - class MappedOutputIterator { + class InsertIterator { + private: + + std::map& _host; + const _solver_bits::VarIndex& _index; + 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) {} + InsertIterator(std::map& host, + const _solver_bits::VarIndex& index) + : _host(host), _index(index) {} - MappedOutputIterator& operator*() { + InsertIterator& operator=(const std::pair& value) { + typedef std::map::value_type pair_type; + _host.insert(pair_type(_index[value.first], value.second)); return *this; } - MappedOutputIterator& operator=(const std::pair& value) { - *base = std::make_pair(lp._item(value.first, typename _Expr::Key()), - value.second); - return *this; - } + InsertIterator& operator*() { return *this; } + InsertIterator& operator++() { return *this; } + InsertIterator operator++(int) { 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 { + class ExprIterator { + private: + std::map::const_iterator _host_it; + const _solver_bits::VarIndex& _index; public: - typedef typename Expr::const_iterator Base; - - typedef typename Base::iterator_category iterator_category; - typedef typename Base::difference_type difference_type; + typedef std::bidirectional_iterator_tag iterator_category; + typedef std::ptrdiff_t difference_type; typedef const std::pair value_type; typedef value_type reference; + class pointer { public: pointer(value_type& _value) : value(_value) {} @@ -653,82 +883,49 @@ value_type value; }; - MappedInputIterator(const Base& _base, const LpSolverBase& _lp) - : base(_base), lp(_lp) {} + ExprIterator(const std::map::const_iterator& host_it, + const _solver_bits::VarIndex& index) + : _host_it(host_it), _index(index) {} reference operator*() { - return std::make_pair(lp._lpId(base->first), base->second); + return std::make_pair(_index(_host_it->first), _host_it->second); } pointer operator->() { return pointer(operator*()); } - MappedInputIterator& operator++() { - ++base; - return *this; + ExprIterator& operator++() { ++_host_it; return *this; } + ExprIterator operator++(int) { + ExprIterator tmp(*this); ++_host_it; return tmp; } - MappedInputIterator operator++(int) { - MappedInputIterator tmp(*this); - ++base; - return tmp; + ExprIterator& operator--() { --_host_it; return *this; } + ExprIterator operator--(int) { + ExprIterator tmp(*this); --_host_it; return tmp; } - bool operator==(const MappedInputIterator& it) const { - return base == it.base; + bool operator==(const ExprIterator& it) const { + return _host_it == it._host_it; } - bool operator!=(const MappedInputIterator& it) const { - return base != it.base; + bool operator!=(const ExprIterator& it) const { + return _host_it != it._host_it; } - 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; + //Abstract virtual functions + virtual LpBase* _newSolver() const = 0; + virtual LpBase* _cloneSolver() const = 0; - /// STL compatible iterator for lp col - typedef MappedOutputIterator RowIterator; - /// STL compatible iterator for lp row - typedef MappedOutputIterator ColIterator; + virtual int _addColId(int col) { return cols.addIndex(col); } + virtual int _addRowId(int row) { return rows.addIndex(row); } - //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 void _eraseColId(int col) { cols.eraseIndex(col); } + virtual void _eraseRowId(int row) { rows.eraseIndex(row); } virtual int _addCol() = 0; virtual int _addRow() = 0; @@ -736,93 +933,95 @@ 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 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 _getRowName(int row, std::string& name) const = 0; + virtual void _setRowName(int row, const std::string& name) = 0; + virtual int _rowByName(const std::string& name) const = 0; + + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0; + virtual void _getRowCoeffs(int i, InsertIterator b) const = 0; + + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0; + virtual void _getColCoeffs(int i, InsertIterator 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 _setRowLowerBound(int i, Value value) = 0; + virtual Value _getRowLowerBound(int i) const = 0; + + virtual void _setRowUpperBound(int i, Value value) = 0; + virtual Value _getRowUpperBound(int i) const = 0; + + virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0; + virtual void _getObjCoeffs(InsertIterator b) 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 _setSense(Sense) = 0; + virtual Sense _getSense() const = 0; - virtual void _setMax() = 0; - virtual void _setMin() = 0; + virtual void _clear() = 0; - - virtual bool _isMax() const = 0; + virtual const char* _solverName() const = 0; //Own protected stuff //Constant component of the objective function Value obj_const_comp; + LpBase() : rows(), cols(), obj_const_comp(0) {} + public: - ///\e - LpSolverBase() : obj_const_comp(0) {} - - ///\e - virtual ~LpSolverBase() {} + /// Virtual destructor + virtual ~LpBase() {} ///Creates a new LP problem - LpSolverBase* newLp() {return _newLp();} + LpBase* newSolver() {return _newSolver();} ///Makes a copy of the LP problem - LpSolverBase* copyLp() {return _copyLp();} + LpBase* cloneSolver() {return _cloneSolver();} + + ///Gives back the name of the solver. + const char* solverName() const {return _solverName();} ///\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;} + Col addCol() { Col c; c._id = _addColId(_addCol()); return c;} - ///\brief Adds several new columns - ///(i.e a variables) at once + ///\brief Adds several new columns (i.e variables) at once /// - ///This magic function takes a container as its argument - ///and fills its elements - ///with new columns (i.e. variables) + ///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 + ///\ref Col as its \c values_type like ///\code - ///std::vector - ///std::list + ///std::vector + ///std::list ///\endcode ///- a standard STL compatible iterable container with - ///\ref Col as its \c mapped_type - ///like + ///\ref Col as its \c mapped_type like ///\code - ///std::map + ///std::map ///\endcode ///- an iterable lemon \ref concepts::WriteMap "write map" like ///\code - ///ListGraph::NodeMap - ///ListGraph::EdgeMap + ///ListGraph::NodeMap + ///ListGraph::ArcMap ///\endcode ///\return The number of the created column. #ifdef DOXYGEN @@ -830,14 +1029,14 @@ int addColSet(T &t) { return 0;} #else template - typename enable_if::type + 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; @@ -848,7 +1047,7 @@ return s; } template - typename enable_if::type addColSet(T &t,dummy<2> = 2) { int s=0; @@ -866,26 +1065,26 @@ ///\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) { + void col(Col c, const DualExpr &e) { e.simplify(); - _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this), - ConstColIterator(e.end(), *this)); + _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), cols), + ExprIterator(e.comps.end(), cols)); } ///Get a column (i.e a dual constraint) of the LP - ///\param r is the column to get + ///\param c 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)); + _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows)); 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 + ///\param o 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) { @@ -899,32 +1098,28 @@ ///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;} + Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;} - ///\brief Add several new rows - ///(i.e a constraints) at once + ///\brief Add several new rows (i.e constraints) at once /// - ///This magic function takes a container as its argument - ///and fills its elements - ///with new row (i.e. variables) + ///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 + ///\ref Row as its \c values_type like ///\code - ///std::vector - ///std::list + ///std::vector + ///std::list ///\endcode ///- a standard STL compatible iterable container with - ///\ref Row as its \c mapped_type - ///like + ///\ref Row as its \c mapped_type like ///\code - ///std::map + ///std::map ///\endcode ///- an iterable lemon \ref concepts::WriteMap "write map" like ///\code - ///ListGraph::NodeMap - ///ListGraph::EdgeMap + ///ListGraph::NodeMap + ///ListGraph::ArcMap ///\endcode ///\return The number of rows created. #ifdef DOXYGEN @@ -932,16 +1127,15 @@ int addRowSet(T &t) { return 0;} #else template - typename enable_if::type - addRowSet(T &t,dummy<0> = 0) { + 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) { + 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(); @@ -950,9 +1144,8 @@ return s; } template - typename enable_if::type - addRowSet(T &t,dummy<2> = 2) { + typename enable_if::type + addRowSet(T &t, dummy<2> = 2) { int s=0; for(typename T::MapIt i(t); i!=INVALID; ++i) { @@ -969,15 +1162,12 @@ ///\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()); + _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols), + ExprIterator(e.comps.end(), cols)); + _setRowLowerBound(rows(id(r)),l - *e); + _setRowUpperBound(rows(id(r)),u - *e); } ///Set a row (i.e a constraint) of the LP @@ -996,7 +1186,7 @@ ///\return the expression associated to the row Expr row(Row r) const { Expr e; - _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this)); + _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols)); return e; } @@ -1006,8 +1196,6 @@ ///\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); @@ -1023,39 +1211,37 @@ row(r,c); return r; } - ///Erase a coloumn (i.e a variable) from the LP + ///Erase a column (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); + ///\param c is the column to be deleted + void erase(Col c) { + _eraseCol(cols(id(c))); + _eraseColId(cols(id(c))); } - ///Erase a row (i.e a constraint) from the LP + ///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); + void erase(Row r) { + _eraseRow(rows(id(r))); + _eraseRowId(rows(id(r))); } /// Get the name of a column - ///\param c is the coresponding coloumn + ///\param c is the coresponding column ///\return The name of the colunm std::string colName(Col c) const { std::string name; - _getColName(_lpId(c), name); + _getColName(cols(id(c)), name); return name; } /// Set the name of a column - ///\param c is the coresponding coloumn + ///\param c is the coresponding column ///\param name The name to be given void colName(Col c, const std::string& name) { - _setColName(_lpId(c), name); + _setColName(cols(id(c)), name); } /// Get the column by its name @@ -1064,27 +1250,52 @@ ///\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); + return k != -1 ? Col(cols[k]) : Col(INVALID); + } + + /// Get the name of a row + + ///\param r is the coresponding row + ///\return The name of the row + std::string rowName(Row r) const { + std::string name; + _getRowName(rows(id(r)), name); + return name; + } + + /// Set the name of a row + + ///\param r is the coresponding row + ///\param name The name to be given + void rowName(Row r, const std::string& name) { + _setRowName(rows(id(r)), name); + } + + /// Get the row by its name + + ///\param name The name of the row + ///\return the proper row or \c INVALID + Row rowByName(const std::string& name) const { + int k = _rowByName(name); + return k != -1 ? Row(rows[k]) : Row(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 c is the column 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); + _setCoeff(rows(id(r)),cols(id(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 + ///\param r is the row of the element + ///\param c is the column of the element ///\return the corresponding coefficient - Value coeff(Row r, Col c) const { - return _getCoeff(_lpId(r),_lpId(c)); + return _getCoeff(rows(id(r)),cols(id(c))); } /// Set the lower bound of a column (i.e a variable) @@ -1093,39 +1304,39 @@ /// 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); + _setColLowerBound(cols(id(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 function returns the lower bound for column (variable) \c c /// (this might be -\ref INF as well). - ///\return The lower bound for coloumn \t c + ///\return The lower bound for column \c c Value colLowerBound(Col c) const { - return _getColLowerBound(_lpId(c)); + return _getColLowerBound(cols(id(c))); } ///\brief Set the lower bound of several columns - ///(i.e a variables) at once + ///(i.e 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. + ///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 + 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) { @@ -1133,7 +1344,7 @@ } } template - typename enable_if::type colLowerBound(T &t, Value value,dummy<2> = 2) { for(typename T::MapIt i(t); i!=INVALID; ++i){ @@ -1148,39 +1359,39 @@ /// 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); + _setColUpperBound(cols(id(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 function returns the upper bound for column (variable) \c c /// (this might be \ref INF as well). - ///\return The upper bound for coloumn \t c + /// \return The upper bound for column \c c Value colUpperBound(Col c) const { - return _getColUpperBound(_lpId(c)); + return _getColUpperBound(cols(id(c))); } ///\brief Set the upper bound of several columns - ///(i.e a variables) at once + ///(i.e 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. + ///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 + 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) { @@ -1188,7 +1399,7 @@ } } template - typename enable_if::type colUpperBound(T &t, Value value,dummy<2> = 2) { for(typename T::MapIt i(t); i!=INVALID; ++i){ @@ -1204,12 +1415,12 @@ /// 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); + _setColLowerBound(cols(id(c)),lower); + _setColUpperBound(cols(id(c)),upper); } ///\brief Set the lower and the upper bound of several columns - ///(i.e a variables) at once + ///(i.e variables) at once /// ///This magic function takes a container as its argument ///and applies the function on all of its elements. @@ -1222,23 +1433,21 @@ void colBounds(T &t, Value lower, Value upper) { return 0;} #else template - typename enable_if::type + 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 + 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 + 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); @@ -1246,76 +1455,371 @@ } #endif + /// Set the lower bound of a row (i.e a constraint) - /// 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); + /// The lower bound of a constraint (row) has to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or -\ref INF. + void rowLowerBound(Row r, Value value) { + _setRowLowerBound(rows(id(r)),value); } - /// Get the lower and the upper bounds of a row (i.e a constraint) + /// Get the lower bound 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); + /// This function returns the lower bound for row (constraint) \c c + /// (this might be -\ref INF as well). + ///\return The lower bound for row \c r + Value rowLowerBound(Row r) const { + return _getRowLowerBound(rows(id(r))); + } + + /// Set the upper bound of a row (i.e a constraint) + + /// The upper bound of a constraint (row) has to be given by an + /// extended number of type Value, i.e. a finite number of type + /// Value or -\ref INF. + void rowUpperBound(Row r, Value value) { + _setRowUpperBound(rows(id(r)),value); + } + + /// Get the upper bound of a row (i.e a constraint) + + /// This function returns the upper bound for row (constraint) \c c + /// (this might be -\ref INF as well). + ///\return The upper bound for row \c r + Value rowUpperBound(Row r) const { + return _getRowUpperBound(rows(id(r))); } ///Set an element of the objective function - void objCoeff(Col c, Value v) {_setObjCoeff(_lpId(c),v); }; + void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); }; ///Get an element of the objective function - Value objCoeff(Col c) const { return _getObjCoeff(_lpId(c)); }; + Value objCoeff(Col c) const { return _getObjCoeff(cols(id(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(); + /// + void obj(const Expr& e) { + _setObjCoeffs(ExprIterator(e.comps.begin(), cols), + ExprIterator(e.comps.end(), cols)); + obj_const_comp = *e; } ///Get the objective function - ///\return the objective function as a linear expression of type \ref Expr. + ///\return the objective function as a linear expression of type + ///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)); - } - } + _getObjCoeffs(InsertIterator(e.comps, cols)); + *e = obj_const_comp; return e; } - ///Maximize - void max() { _setMax(); } - ///Minimize - void min() { _setMin(); } + ///Set the direction of optimization + void sense(Sense sense) { _setSense(sense); } - ///Query function: is this a maximization problem? - bool isMax() const {return _isMax(); } + ///Query the direction of the optimization + Sense sense() const {return _getSense(); } - ///Query function: is this a minimization problem? - bool isMin() const {return !isMax(); } + ///Set the sense to maximization + void max() { _setSense(MAX); } + + ///Set the sense to maximization + void min() { _setSense(MIN); } + + ///Clears the problem + void clear() { _clear(); } ///@} + }; + + /// Addition + + ///\relates LpBase::Expr + /// + inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) { + LpBase::Expr tmp(a); + tmp+=b; + return tmp; + } + ///Substraction + + ///\relates LpBase::Expr + /// + inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) { + LpBase::Expr tmp(a); + tmp-=b; + return tmp; + } + ///Multiply with constant + + ///\relates LpBase::Expr + /// + inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) { + LpBase::Expr tmp(a); + tmp*=b; + return tmp; + } + + ///Multiply with constant + + ///\relates LpBase::Expr + /// + inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) { + LpBase::Expr tmp(b); + tmp*=a; + return tmp; + } + ///Divide with constant + + ///\relates LpBase::Expr + /// + inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) { + LpBase::Expr tmp(a); + tmp/=b; + return tmp; + } + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator<=(const LpBase::Expr &e, + const LpBase::Expr &f) { + return LpBase::Constr(0, f - e, LpBase::INF); + } + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator<=(const LpBase::Value &e, + const LpBase::Expr &f) { + return LpBase::Constr(e, f, LpBase::NaN); + } + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator<=(const LpBase::Expr &e, + const LpBase::Value &f) { + return LpBase::Constr(- LpBase::INF, e, f); + } + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator>=(const LpBase::Expr &e, + const LpBase::Expr &f) { + return LpBase::Constr(0, e - f, LpBase::INF); + } + + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator>=(const LpBase::Value &e, + const LpBase::Expr &f) { + return LpBase::Constr(LpBase::NaN, f, e); + } + + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator>=(const LpBase::Expr &e, + const LpBase::Value &f) { + return LpBase::Constr(f, e, LpBase::INF); + } + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator==(const LpBase::Expr &e, + const LpBase::Value &f) { + return LpBase::Constr(f, e, f); + } + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator==(const LpBase::Expr &e, + const LpBase::Expr &f) { + return LpBase::Constr(0, f - e, 0); + } + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator<=(const LpBase::Value &n, + const LpBase::Constr &c) { + LpBase::Constr tmp(c); + LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint"); + tmp.lowerBound()=n; + return tmp; + } + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator<=(const LpBase::Constr &c, + const LpBase::Value &n) + { + LpBase::Constr tmp(c); + LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint"); + tmp.upperBound()=n; + return tmp; + } + + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator>=(const LpBase::Value &n, + const LpBase::Constr &c) { + LpBase::Constr tmp(c); + LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint"); + tmp.upperBound()=n; + return tmp; + } + ///Create constraint + + ///\relates LpBase::Constr + /// + inline LpBase::Constr operator>=(const LpBase::Constr &c, + const LpBase::Value &n) + { + LpBase::Constr tmp(c); + LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint"); + tmp.lowerBound()=n; + return tmp; + } + + ///Addition + + ///\relates LpBase::DualExpr + /// + inline LpBase::DualExpr operator+(const LpBase::DualExpr &a, + const LpBase::DualExpr &b) { + LpBase::DualExpr tmp(a); + tmp+=b; + return tmp; + } + ///Substraction + + ///\relates LpBase::DualExpr + /// + inline LpBase::DualExpr operator-(const LpBase::DualExpr &a, + const LpBase::DualExpr &b) { + LpBase::DualExpr tmp(a); + tmp-=b; + return tmp; + } + ///Multiply with constant + + ///\relates LpBase::DualExpr + /// + inline LpBase::DualExpr operator*(const LpBase::DualExpr &a, + const LpBase::Value &b) { + LpBase::DualExpr tmp(a); + tmp*=b; + return tmp; + } + + ///Multiply with constant + + ///\relates LpBase::DualExpr + /// + inline LpBase::DualExpr operator*(const LpBase::Value &a, + const LpBase::DualExpr &b) { + LpBase::DualExpr tmp(b); + tmp*=a; + return tmp; + } + ///Divide with constant + + ///\relates LpBase::DualExpr + /// + inline LpBase::DualExpr operator/(const LpBase::DualExpr &a, + const LpBase::Value &b) { + LpBase::DualExpr tmp(a); + tmp/=b; + return tmp; + } + + /// \ingroup lp_group + /// + /// \brief Common base class for LP solvers + /// + /// This class is an abstract base class for LP solvers. This class + /// provides a full interface for set and modify an LP problem, + /// solve it and retrieve the solution. You can use one of the + /// descendants as a concrete implementation, or the \c Lp + /// default LP solver. However, if you would like to handle LP + /// solvers as reference or pointer in a generic way, you can use + /// this class directly. + class LpSolver : virtual public LpBase { + public: + + /// The problem types for primal and dual problems + enum ProblemType { + ///Feasible solution hasn't been found (but may exist). + 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 + UNBOUNDED = 4 + }; + + ///The basis status of variables + enum VarStatus { + /// The variable is in the basis + BASIC, + /// The variable is free, but not basic + FREE, + /// The variable has active lower bound + LOWER, + /// The variable has active upper bound + UPPER, + /// The variable is non-basic and fixed + FIXED + }; + + protected: + + virtual SolveExitStatus _solve() = 0; + + virtual Value _getPrimal(int i) const = 0; + virtual Value _getDual(int i) const = 0; + + virtual Value _getPrimalRay(int i) const = 0; + virtual Value _getDualRay(int i) const = 0; + + virtual Value _getPrimalValue() const = 0; + + virtual VarStatus _getColStatus(int i) const = 0; + virtual VarStatus _getRowStatus(int i) const = 0; + + virtual ProblemType _getPrimalType() const = 0; + virtual ProblemType _getDualType() const = 0; + + public: ///\name Solve the LP @@ -1326,8 +1830,6 @@ ///\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(); } ///@} @@ -1336,368 +1838,241 @@ ///@{ - /// The status of the primal problem (the original LP problem) - SolutionStatus primalStatus() const { - return _getPrimalStatus(); + /// The type of the primal problem + ProblemType primalType() const { + return _getPrimalType(); } - /// The status of the dual (of the original LP) problem - SolutionStatus dualStatus() const { - return _getDualStatus(); + /// The type of the dual problem + ProblemType dualType() const { + return _getDualType(); } - ///The type of the original LP problem - ProblemTypes problemType() const { - return _getProblemType(); + /// Return the primal value of the column + + /// Return the primal value of the column. + /// \pre The problem is solved. + Value primal(Col c) const { return _getPrimal(cols(id(c))); } + + /// Return the primal value of the expression + + /// Return the primal value of the expression, i.e. the dot + /// product of the primal solution and the expression. + /// \pre The problem is solved. + Value primal(const Expr& e) const { + double res = *e; + for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) { + res += *c * primal(c); + } + return res; } + /// Returns a component of the primal ray + + /// The primal ray is solution of the modified primal problem, + /// where we change each finite bound to 0, and we looking for a + /// negative objective value in case of minimization, and positive + /// objective value for maximization. If there is such solution, + /// that proofs the unsolvability of the dual problem, and if a + /// feasible primal solution exists, then the unboundness of + /// primal problem. + /// + /// \pre The problem is solved and the dual problem is infeasible. + /// \note Some solvers does not provide primal ray calculation + /// functions. + Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); } - ///\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 the dual value of the row + + /// Return the dual value of the row. + /// \pre The problem is solved. + Value dual(Row r) const { return _getDual(rows(id(r))); } + + /// Return the dual value of the dual expression + + /// Return the dual value of the dual expression, i.e. the dot + /// product of the dual solution and the dual expression. + /// \pre The problem is solved. + Value dual(const DualExpr& e) const { + double res = 0.0; + for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) { + res += *r * dual(r); } 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; - } + /// Returns a component of the dual ray + + /// The dual ray is solution of the modified primal problem, where + /// we change each finite bound to 0 (i.e. the objective function + /// coefficients in the primal problem), and we looking for a + /// ositive objective value. If there is such solution, that + /// proofs the unsolvability of the primal problem, and if a + /// feasible dual solution exists, then the unboundness of + /// dual problem. + /// + /// \pre The problem is solved and the primal problem is infeasible. + /// \note Some solvers does not provide dual ray calculation + /// functions. + Value dualRay(Row r) const { return _getDualRay(rows(id(r))); } - ///\e - bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); } + /// Return the basis status of the column - ///\e + /// \see VarStatus + VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); } + + /// Return the basis status of the row + + /// \see VarStatus + VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); } + + ///The value of the objective function ///\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;} + Value primal() const { return _getPrimalValue()+obj_const_comp;} ///@} + LpSolver* newSolver() {return _newSolver();} + LpSolver* cloneSolver() {return _cloneSolver();} + + protected: + + virtual LpSolver* _newSolver() const = 0; + virtual LpSolver* _cloneSolver() const = 0; }; /// \ingroup lp_group /// /// \brief Common base class for MIP solvers - /// \todo Much more docs - class MipSolverBase : virtual public LpSolverBase{ + /// + /// This class is an abstract base class for MIP solvers. This class + /// provides a full interface for set and modify an MIP problem, + /// solve it and retrieve the solution. You can use one of the + /// descendants as a concrete implementation, or the \c Lp + /// default MIP solver. However, if you would like to handle MIP + /// solvers as reference or pointer in a generic way, you can use + /// this class directly. + class MipSolver : virtual public LpBase { public: - ///Possible variable (coloumn) types (e.g. real, integer, binary etc.) + /// The problem types for MIP problems + enum ProblemType { + ///Feasible solution hasn't been found (but may exist). + 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 + /// + ///The Mip or at least the relaxed problem is unbounded + UNBOUNDED = 4 + }; + + ///\name Solve the MIP + + ///@{ + + /// Solve the MIP problem at hand + /// + ///\return The result of the optimization procedure. Possible + ///values and their meanings can be found in the documentation of + ///\ref SolveExitStatus. + SolveExitStatus solve() { return _solve(); } + + ///@} + + ///\name Setting column type + ///@{ + + ///Possible variable (column) types (e.g. real, integer, binary etc.) enum ColTypes { - ///Continuous variable + ///Continuous variable (default) REAL = 0, ///Integer variable - - ///Unfortunately, cplex 7.5 somewhere writes something like - ///#define INTEGER 'I' - INT = 1 - ///\todo No support for other types yet. + INTEGER = 1 }; - ///Sets the type of the given coloumn to the given type + ///Sets the type of the given column to the given type + + ///Sets the type of the given column 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); + _setColType(cols(id(c)),col_type); } ///Gives back the type of the column. + + ///Gives back the type of the column. /// - ///Gives back the type of the column. ColTypes colType(Col c) const { - return _colType(_lpId(c)); + return _getColType(cols(id(c))); + } + ///@} + + ///\name Obtain the solution + + ///@{ + + /// The type of the MIP problem + ProblemType type() const { + return _getType(); } - ///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); + /// Return the value of the row in the solution + + /// Return the value of the row in the solution. + /// \pre The problem is solved. + Value sol(Col c) const { return _getSol(cols(id(c))); } + + /// Return the value of the expression in the solution + + /// Return the value of the expression in the solution, i.e. the + /// dot product of the solution and the expression. + /// \pre The problem is solved. + Value sol(const Expr& e) const { + double res = *e; + for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) { + res += *c * sol(c); + } + return res; } - - ///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(); - } + ///The value of the objective function + + ///\return + ///- \ref INF or -\ref INF means either infeasibility or unboundedness + /// of the 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 solValue() const { return _getSolValue()+obj_const_comp;} + ///@} protected: - virtual ColTypes _colType(int col) const = 0; - virtual void _colType(int col, ColTypes col_type) = 0; - virtual SolutionStatus _getMipStatus() const = 0; + virtual SolveExitStatus _solve() = 0; + virtual ColTypes _getColType(int col) const = 0; + virtual void _setColType(int col, ColTypes col_type) = 0; + virtual ProblemType _getType() const = 0; + virtual Value _getSol(int i) const = 0; + virtual Value _getSolValue() const = 0; + public: + + MipSolver* newSolver() {return _newSolver();} + MipSolver* cloneSolver() {return _cloneSolver();} + + protected: + + virtual MipSolver* _newSolver() const = 0; + virtual MipSolver* _cloneSolver() 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 diff --git a/lemon/lp_clp.cc b/lemon/lp_clp.cc new file mode 100644 --- /dev/null +++ b/lemon/lp_clp.cc @@ -0,0 +1,437 @@ +/* -*- 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 + +namespace lemon { + + LpClp::LpClp() { + _prob = new ClpSimplex(); + _init_temporals(); + messageLevel(MESSAGE_NO_OUTPUT); + } + + LpClp::LpClp(const LpClp& other) { + _prob = new ClpSimplex(*other._prob); + rows = other.rows; + cols = other.cols; + _init_temporals(); + messageLevel(MESSAGE_NO_OUTPUT); + } + + LpClp::~LpClp() { + delete _prob; + _clear_temporals(); + } + + void LpClp::_init_temporals() { + _primal_ray = 0; + _dual_ray = 0; + } + + void LpClp::_clear_temporals() { + if (_primal_ray) { + delete[] _primal_ray; + _primal_ray = 0; + } + if (_dual_ray) { + delete[] _dual_ray; + _dual_ray = 0; + } + } + + LpClp* LpClp::_newSolver() const { + LpClp* newlp = new LpClp; + return newlp; + } + + LpClp* LpClp::_cloneSolver() const { + LpClp* copylp = new LpClp(*this); + return copylp; + } + + const char* LpClp::_solverName() const { return "LpClp"; } + + int LpClp::_addCol() { + _prob->addColumn(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX, 0.0); + return _prob->numberColumns() - 1; + } + + int LpClp::_addRow() { + _prob->addRow(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX); + return _prob->numberRows() - 1; + } + + + void LpClp::_eraseCol(int c) { + _col_names_ref.erase(_prob->getColumnName(c)); + _prob->deleteColumns(1, &c); + } + + void LpClp::_eraseRow(int r) { + _row_names_ref.erase(_prob->getRowName(r)); + _prob->deleteRows(1, &r); + } + + void LpClp::_eraseColId(int i) { + cols.eraseIndex(i); + cols.shiftIndices(i); + } + + void LpClp::_eraseRowId(int i) { + rows.eraseIndex(i); + rows.shiftIndices(i); + } + + void LpClp::_getColName(int c, std::string& name) const { + name = _prob->getColumnName(c); + } + + void LpClp::_setColName(int c, const std::string& name) { + _prob->setColumnName(c, const_cast(name)); + _col_names_ref[name] = c; + } + + int LpClp::_colByName(const std::string& name) const { + std::map::const_iterator it = _col_names_ref.find(name); + return it != _col_names_ref.end() ? it->second : -1; + } + + void LpClp::_getRowName(int r, std::string& name) const { + name = _prob->getRowName(r); + } + + void LpClp::_setRowName(int r, const std::string& name) { + _prob->setRowName(r, const_cast(name)); + _row_names_ref[name] = r; + } + + int LpClp::_rowByName(const std::string& name) const { + std::map::const_iterator it = _row_names_ref.find(name); + return it != _row_names_ref.end() ? it->second : -1; + } + + + void LpClp::_setRowCoeffs(int ix, ExprIterator b, ExprIterator e) { + std::map coeffs; + + int n = _prob->clpMatrix()->getNumCols(); + + const int* indices = _prob->clpMatrix()->getIndices(); + const double* elements = _prob->clpMatrix()->getElements(); + + for (int i = 0; i < n; ++i) { + CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i]; + CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i]; + + const int* it = std::lower_bound(indices + begin, indices + end, ix); + if (it != indices + end && *it == ix && elements[it - indices] != 0.0) { + coeffs[i] = 0.0; + } + } + + for (ExprIterator it = b; it != e; ++it) { + coeffs[it->first] = it->second; + } + + for (std::map::iterator it = coeffs.begin(); + it != coeffs.end(); ++it) { + _prob->modifyCoefficient(ix, it->first, it->second); + } + } + + void LpClp::_getRowCoeffs(int ix, InsertIterator b) const { + int n = _prob->clpMatrix()->getNumCols(); + + const int* indices = _prob->clpMatrix()->getIndices(); + const double* elements = _prob->clpMatrix()->getElements(); + + for (int i = 0; i < n; ++i) { + CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i]; + CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i]; + + const int* it = std::lower_bound(indices + begin, indices + end, ix); + if (it != indices + end && *it == ix) { + *b = std::make_pair(i, elements[it - indices]); + } + } + } + + void LpClp::_setColCoeffs(int ix, ExprIterator b, ExprIterator e) { + std::map coeffs; + + CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix]; + CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix]; + + const int* indices = _prob->clpMatrix()->getIndices(); + const double* elements = _prob->clpMatrix()->getElements(); + + for (CoinBigIndex i = begin; i != end; ++i) { + if (elements[i] != 0.0) { + coeffs[indices[i]] = 0.0; + } + } + for (ExprIterator it = b; it != e; ++it) { + coeffs[it->first] = it->second; + } + for (std::map::iterator it = coeffs.begin(); + it != coeffs.end(); ++it) { + _prob->modifyCoefficient(it->first, ix, it->second); + } + } + + void LpClp::_getColCoeffs(int ix, InsertIterator b) const { + CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix]; + CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix]; + + const int* indices = _prob->clpMatrix()->getIndices(); + const double* elements = _prob->clpMatrix()->getElements(); + + for (CoinBigIndex i = begin; i != end; ++i) { + *b = std::make_pair(indices[i], elements[i]); + ++b; + } + } + + void LpClp::_setCoeff(int ix, int jx, Value value) { + _prob->modifyCoefficient(ix, jx, value); + } + + LpClp::Value LpClp::_getCoeff(int ix, int jx) const { + CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix]; + CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix]; + + const int* indices = _prob->clpMatrix()->getIndices(); + const double* elements = _prob->clpMatrix()->getElements(); + + const int* it = std::lower_bound(indices + begin, indices + end, jx); + if (it != indices + end && *it == jx) { + return elements[it - indices]; + } else { + return 0.0; + } + } + + void LpClp::_setColLowerBound(int i, Value lo) { + _prob->setColumnLower(i, lo == - INF ? - COIN_DBL_MAX : lo); + } + + LpClp::Value LpClp::_getColLowerBound(int i) const { + double val = _prob->getColLower()[i]; + return val == - COIN_DBL_MAX ? - INF : val; + } + + void LpClp::_setColUpperBound(int i, Value up) { + _prob->setColumnUpper(i, up == INF ? COIN_DBL_MAX : up); + } + + LpClp::Value LpClp::_getColUpperBound(int i) const { + double val = _prob->getColUpper()[i]; + return val == COIN_DBL_MAX ? INF : val; + } + + void LpClp::_setRowLowerBound(int i, Value lo) { + _prob->setRowLower(i, lo == - INF ? - COIN_DBL_MAX : lo); + } + + LpClp::Value LpClp::_getRowLowerBound(int i) const { + double val = _prob->getRowLower()[i]; + return val == - COIN_DBL_MAX ? - INF : val; + } + + void LpClp::_setRowUpperBound(int i, Value up) { + _prob->setRowUpper(i, up == INF ? COIN_DBL_MAX : up); + } + + LpClp::Value LpClp::_getRowUpperBound(int i) const { + double val = _prob->getRowUpper()[i]; + return val == COIN_DBL_MAX ? INF : val; + } + + void LpClp::_setObjCoeffs(ExprIterator b, ExprIterator e) { + int num = _prob->clpMatrix()->getNumCols(); + for (int i = 0; i < num; ++i) { + _prob->setObjectiveCoefficient(i, 0.0); + } + for (ExprIterator it = b; it != e; ++it) { + _prob->setObjectiveCoefficient(it->first, it->second); + } + } + + void LpClp::_getObjCoeffs(InsertIterator b) const { + int num = _prob->clpMatrix()->getNumCols(); + for (int i = 0; i < num; ++i) { + Value coef = _prob->getObjCoefficients()[i]; + if (coef != 0.0) { + *b = std::make_pair(i, coef); + ++b; + } + } + } + + void LpClp::_setObjCoeff(int i, Value obj_coef) { + _prob->setObjectiveCoefficient(i, obj_coef); + } + + LpClp::Value LpClp::_getObjCoeff(int i) const { + return _prob->getObjCoefficients()[i]; + } + + LpClp::SolveExitStatus LpClp::_solve() { + return _prob->primal() >= 0 ? SOLVED : UNSOLVED; + } + + LpClp::SolveExitStatus LpClp::solvePrimal() { + return _prob->primal() >= 0 ? SOLVED : UNSOLVED; + } + + LpClp::SolveExitStatus LpClp::solveDual() { + return _prob->dual() >= 0 ? SOLVED : UNSOLVED; + } + + LpClp::SolveExitStatus LpClp::solveBarrier() { + return _prob->barrier() >= 0 ? SOLVED : UNSOLVED; + } + + LpClp::Value LpClp::_getPrimal(int i) const { + return _prob->primalColumnSolution()[i]; + } + LpClp::Value LpClp::_getPrimalValue() const { + return _prob->objectiveValue(); + } + + LpClp::Value LpClp::_getDual(int i) const { + return _prob->dualRowSolution()[i]; + } + + LpClp::Value LpClp::_getPrimalRay(int i) const { + if (!_primal_ray) { + _primal_ray = _prob->unboundedRay(); + LEMON_ASSERT(_primal_ray != 0, "Primal ray is not provided"); + } + return _primal_ray[i]; + } + + LpClp::Value LpClp::_getDualRay(int i) const { + if (!_dual_ray) { + _dual_ray = _prob->infeasibilityRay(); + LEMON_ASSERT(_dual_ray != 0, "Dual ray is not provided"); + } + return _dual_ray[i]; + } + + LpClp::VarStatus LpClp::_getColStatus(int i) const { + switch (_prob->getColumnStatus(i)) { + case ClpSimplex::basic: + return BASIC; + case ClpSimplex::isFree: + return FREE; + case ClpSimplex::atUpperBound: + return UPPER; + case ClpSimplex::atLowerBound: + return LOWER; + case ClpSimplex::isFixed: + return FIXED; + case ClpSimplex::superBasic: + return FREE; + default: + LEMON_ASSERT(false, "Wrong column status"); + return VarStatus(); + } + } + + LpClp::VarStatus LpClp::_getRowStatus(int i) const { + switch (_prob->getColumnStatus(i)) { + case ClpSimplex::basic: + return BASIC; + case ClpSimplex::isFree: + return FREE; + case ClpSimplex::atUpperBound: + return UPPER; + case ClpSimplex::atLowerBound: + return LOWER; + case ClpSimplex::isFixed: + return FIXED; + case ClpSimplex::superBasic: + return FREE; + default: + LEMON_ASSERT(false, "Wrong row status"); + return VarStatus(); + } + } + + + LpClp::ProblemType LpClp::_getPrimalType() const { + if (_prob->isProvenOptimal()) { + return OPTIMAL; + } else if (_prob->isProvenPrimalInfeasible()) { + return INFEASIBLE; + } else if (_prob->isProvenDualInfeasible()) { + return UNBOUNDED; + } else { + return UNDEFINED; + } + } + + LpClp::ProblemType LpClp::_getDualType() const { + if (_prob->isProvenOptimal()) { + return OPTIMAL; + } else if (_prob->isProvenDualInfeasible()) { + return INFEASIBLE; + } else if (_prob->isProvenPrimalInfeasible()) { + return INFEASIBLE; + } else { + return UNDEFINED; + } + } + + void LpClp::_setSense(LpClp::Sense sense) { + switch (sense) { + case MIN: + _prob->setOptimizationDirection(1); + break; + case MAX: + _prob->setOptimizationDirection(-1); + break; + } + } + + LpClp::Sense LpClp::_getSense() const { + double dir = _prob->optimizationDirection(); + if (dir > 0.0) { + return MIN; + } else { + return MAX; + } + } + + void LpClp::_clear() { + delete _prob; + _prob = new ClpSimplex(); + rows.clear(); + cols.clear(); + _col_names_ref.clear(); + _clear_temporals(); + } + + void LpClp::messageLevel(MessageLevel m) { + _prob->setLogLevel(static_cast(m)); + } + +} //END OF NAMESPACE LEMON diff --git a/lemon/lp_clp.h b/lemon/lp_clp.h new file mode 100644 --- /dev/null +++ b/lemon/lp_clp.h @@ -0,0 +1,179 @@ +/* -*- 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_CLP_H +#define LEMON_LP_CLP_H + +///\file +///\brief Header of the LEMON-CLP lp solver interface. + +#include +#include + +#include + +class ClpSimplex; + +namespace lemon { + + /// \ingroup lp_group + /// + /// \brief Interface for the CLP solver + /// + /// This class implements an interface for the Clp LP solver. The + /// Clp library is an object oriented lp solver library developed at + /// the IBM. The CLP is part of the COIN-OR package and it can be + /// used with Common Public License. + class LpClp : public LpSolver { + protected: + + ClpSimplex* _prob; + + std::map _col_names_ref; + std::map _row_names_ref; + + public: + + /// \e + LpClp(); + /// \e + LpClp(const LpClp&); + /// \e + ~LpClp(); + + protected: + + mutable double* _primal_ray; + mutable double* _dual_ray; + + void _init_temporals(); + void _clear_temporals(); + + protected: + + virtual LpClp* _newSolver() const; + virtual LpClp* _cloneSolver() const; + + virtual const char* _solverName() const; + + virtual int _addCol(); + virtual int _addRow(); + + virtual void _eraseCol(int i); + virtual void _eraseRow(int i); + + virtual void _eraseColId(int i); + virtual void _eraseRowId(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 _getRowName(int row, std::string& name) const; + virtual void _setRowName(int row, const std::string& name); + virtual int _rowByName(const std::string& name) const; + + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e); + virtual void _getRowCoeffs(int i, InsertIterator b) const; + + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e); + virtual void _getColCoeffs(int i, InsertIterator 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 Value _getRowLowerBound(int i) const; + virtual void _setRowUpperBound(int i, Value value); + virtual Value _getRowUpperBound(int i) const; + + virtual void _setObjCoeffs(ExprIterator, ExprIterator); + virtual void _getObjCoeffs(InsertIterator) const; + + virtual void _setObjCoeff(int i, Value obj_coef); + virtual Value _getObjCoeff(int i) const; + + virtual void _setSense(Sense sense); + virtual Sense _getSense() const; + + virtual SolveExitStatus _solve(); + + virtual Value _getPrimal(int i) const; + virtual Value _getDual(int i) const; + + virtual Value _getPrimalValue() const; + + virtual Value _getPrimalRay(int i) const; + virtual Value _getDualRay(int i) const; + + virtual VarStatus _getColStatus(int i) const; + virtual VarStatus _getRowStatus(int i) const; + + virtual ProblemType _getPrimalType() const; + virtual ProblemType _getDualType() const; + + virtual void _clear(); + + public: + + ///Solves LP with primal simplex method. + SolveExitStatus solvePrimal(); + + ///Solves LP with dual simplex method. + SolveExitStatus solveDual(); + + ///Solves LP with barrier method. + SolveExitStatus solveBarrier(); + + ///Returns the constraint identifier understood by CLP. + int clpRow(Row r) const { return rows(id(r)); } + + ///Returns the variable identifier understood by CLP. + int clpCol(Col c) const { return cols(id(c)); } + + ///Enum for \c messageLevel() parameter + enum MessageLevel { + /// no output (default value) + MESSAGE_NO_OUTPUT = 0, + /// print final solution + MESSAGE_FINAL_SOLUTION = 1, + /// print factorization + MESSAGE_FACTORIZATION = 2, + /// normal output + MESSAGE_NORMAL_OUTPUT = 3, + /// verbose output + MESSAGE_VERBOSE_OUTPUT = 4 + }; + ///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. + void messageLevel(MessageLevel m); + + }; + +} //END OF NAMESPACE LEMON + +#endif //LEMON_LP_CLP_H + diff --git a/lemon/lp_cplex.cc b/lemon/lp_cplex.cc --- a/lemon/lp_cplex.cc +++ b/lemon/lp_cplex.cc @@ -18,6 +18,8 @@ #include #include +#include + #include extern "C" { @@ -29,167 +31,226 @@ ///\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"); + CplexEnv::LicenseError::LicenseError(int status) { + if (!CPXgeterrorstring(0, status, _message)) { + std::strcpy(_message, "Cplex unknown error"); + } } - LpCplex::LpCplex(const LpCplex& cplex) : LpSolverBase() { - env = CPXopenCPLEX(&status); - lp = CPXcloneprob(env, cplex.lp, &status); + CplexEnv::CplexEnv() { + int status; + _cnt = new int; + _env = CPXopenCPLEX(&status); + if (_env == 0) { + delete _cnt; + _cnt = 0; + throw LicenseError(status); + } + } + + CplexEnv::CplexEnv(const CplexEnv& other) { + _env = other._env; + _cnt = other._cnt; + ++(*_cnt); + } + + CplexEnv& CplexEnv::operator=(const CplexEnv& other) { + _env = other._env; + _cnt = other._cnt; + ++(*_cnt); + return *this; + } + + CplexEnv::~CplexEnv() { + --(*_cnt); + if (*_cnt == 0) { + delete _cnt; + CPXcloseCPLEX(&_env); + } + } + + CplexBase::CplexBase() : LpBase() { + int status; + _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem"); + } + + CplexBase::CplexBase(const CplexEnv& env) + : LpBase(), _env(env) { + int status; + _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem"); + } + + CplexBase::CplexBase(const CplexBase& cplex) + : LpBase() { + int status; + _prob = CPXcloneprob(cplexEnv(), cplex._prob, &status); rows = cplex.rows; cols = cplex.cols; } - LpCplex::~LpCplex() { - CPXfreeprob(env,&lp); - CPXcloseCPLEX(&env); + CplexBase::~CplexBase() { + CPXfreeprob(cplexEnv(),&_prob); } - 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); + int CplexBase::_addCol() { + int i = CPXgetnumcols(cplexEnv(), _prob); + double lb = -INF, ub = INF; + CPXnewcols(cplexEnv(), _prob, 1, 0, &lb, &ub, 0, 0); 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); + int CplexBase::_addRow() { + int i = CPXgetnumrows(cplexEnv(), _prob); + const double ub = INF; + const char s = 'L'; + CPXnewrows(cplexEnv(), _prob, 1, &ub, &s, 0, 0); return i; } - void LpCplex::_eraseCol(int i) { - CPXdelcols(env, lp, i, i); + void CplexBase::_eraseCol(int i) { + CPXdelcols(cplexEnv(), _prob, i, i); } - void LpCplex::_eraseRow(int i) { - CPXdelrows(env, lp, i, i); + void CplexBase::_eraseRow(int i) { + CPXdelrows(cplexEnv(), _prob, 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) { + void CplexBase::_eraseColId(int i) { + cols.eraseIndex(i); + cols.shiftIndices(i); + } + void CplexBase::_eraseRowId(int i) { + rows.eraseIndex(i); + rows.shiftIndices(i); + } + + void CplexBase::_getColName(int col, std::string &name) const { + int size; + CPXgetcolname(cplexEnv(), _prob, 0, 0, 0, &size, col, col); + if (size == 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]; + size *= -1; + std::vector buf(size); + char *cname; + int tmp; + CPXgetcolname(cplexEnv(), _prob, &cname, &buf.front(), size, + &tmp, col, col); + name = cname; } - 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); + void CplexBase::_setColName(int col, const std::string &name) { + char *cname; + cname = const_cast(name.c_str()); + CPXchgcolname(cplexEnv(), _prob, 1, &col, &cname); } - int LpCplex::_colByName(const std::string& name) const - { + int CplexBase::_colByName(const std::string& name) const { int index; - if (CPXgetcolindex(env, lp, + if (CPXgetcolindex(cplexEnv(), _prob, 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) + void CplexBase::_getRowName(int row, std::string &name) const { + int size; + CPXgetrowname(cplexEnv(), _prob, 0, 0, 0, &size, row, row); + if (size == 0) { + name.clear(); + return; + } + + size *= -1; + std::vector buf(size); + char *cname; + int tmp; + CPXgetrowname(cplexEnv(), _prob, &cname, &buf.front(), size, + &tmp, row, row); + name = cname; + } + + void CplexBase::_setRowName(int row, const std::string &name) { + char *cname; + cname = const_cast(name.c_str()); + CPXchgrowname(cplexEnv(), _prob, 1, &row, &cname); + } + + int CplexBase::_rowByName(const std::string& name) const { + int index; + if (CPXgetrowindex(cplexEnv(), _prob, + const_cast(name.c_str()), &index) == 0) { + return index; + } + return -1; + } + + void CplexBase::_setRowCoeffs(int i, ExprIterator b, + ExprIterator e) { std::vector indices; std::vector rowlist; std::vector values; - for(ConstRowIterator it=b; it!=e; ++it) { + for(ExprIterator 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]); + CPXchgcoeflist(cplexEnv(), _prob, values.size(), + &rowlist.front(), &indices.front(), &values.front()); } - void LpCplex::_getRowCoeffs(int i, RowIterator b) const { + void CplexBase::_getRowCoeffs(int i, InsertIterator b) const { int tmp1, tmp2, tmp3, length; - CPXgetrows(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i); + CPXgetrows(cplexEnv(), _prob, &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], + CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2, + &indices.front(), &values.front(), 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) - { + void CplexBase::_setColCoeffs(int i, ExprIterator b, ExprIterator e) { std::vector indices; std::vector collist; std::vector values; - for(ConstColIterator it=b; it!=e; ++it) { + for(ExprIterator 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]); + CPXchgcoeflist(cplexEnv(), _prob, values.size(), + &indices.front(), &collist.front(), &values.front()); } - void LpCplex::_getColCoeffs(int i, ColIterator b) const { + void CplexBase::_getColCoeffs(int i, InsertIterator b) const { int tmp1, tmp2, tmp3, length; - CPXgetcols(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i); + CPXgetcols(cplexEnv(), _prob, &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], + CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2, + &indices.front(), &values.front(), length, &tmp3, i, i); for (int i = 0; i < length; ++i) { @@ -199,175 +260,209 @@ } - void LpCplex::_setCoeff(int row, int col, Value value) - { - CPXchgcoef(env, lp, row, col, value); + void CplexBase::_setCoeff(int row, int col, Value value) { + CPXchgcoef(cplexEnv(), _prob, row, col, value); } - LpCplex::Value LpCplex::_getCoeff(int row, int col) const - { - LpCplex::Value value; - CPXgetcoef(env, lp, row, col, &value); + CplexBase::Value CplexBase::_getCoeff(int row, int col) const { + CplexBase::Value value; + CPXgetcoef(cplexEnv(), _prob, row, col, &value); return value; } - void LpCplex::_setColLowerBound(int i, Value value) + void CplexBase::_setColLowerBound(int i, Value value) { + const char s = 'L'; + CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value); + } + + CplexBase::Value CplexBase::_getColLowerBound(int i) const { + CplexBase::Value res; + CPXgetlb(cplexEnv(), _prob, &res, i, i); + return res <= -CPX_INFBOUND ? -INF : res; + } + + void CplexBase::_setColUpperBound(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); + const char s = 'U'; + CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value); + } + + CplexBase::Value CplexBase::_getColUpperBound(int i) const { + CplexBase::Value res; + CPXgetub(cplexEnv(), _prob, &res, i, i); + return res >= CPX_INFBOUND ? INF : res; + } + + CplexBase::Value CplexBase::_getRowLowerBound(int i) const { + char s; + CPXgetsense(cplexEnv(), _prob, &s, i, i); + CplexBase::Value res; + + switch (s) { + case 'G': + case 'R': + case 'E': + CPXgetrhs(cplexEnv(), _prob, &res, i, i); + return res <= -CPX_INFBOUND ? -INF : res; + default: + return -INF; + } + } + + CplexBase::Value CplexBase::_getRowUpperBound(int i) const { + char s; + CPXgetsense(cplexEnv(), _prob, &s, i, i); + CplexBase::Value res; + + switch (s) { + case 'L': + case 'E': + CPXgetrhs(cplexEnv(), _prob, &res, i, i); + return res >= CPX_INFBOUND ? INF : res; + case 'R': + CPXgetrhs(cplexEnv(), _prob, &res, i, i); + { + double rng; + CPXgetrngval(cplexEnv(), _prob, &rng, i, i); + res += rng; + } + return res >= CPX_INFBOUND ? INF : res; + default: + return INF; + } + } + + //This is easier to implement + void CplexBase::_set_row_bounds(int i, Value lb, Value ub) { + if (lb == -INF) { + const char s = 'L'; + CPXchgsense(cplexEnv(), _prob, 1, &i, &s); + CPXchgrhs(cplexEnv(), _prob, 1, &i, &ub); + } else if (ub == INF) { + const char s = 'G'; + CPXchgsense(cplexEnv(), _prob, 1, &i, &s); + CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb); + } else if (lb == ub){ + const char s = 'E'; + CPXchgsense(cplexEnv(), _prob, 1, &i, &s); + CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb); + } else { + const char s = 'R'; + CPXchgsense(cplexEnv(), _prob, 1, &i, &s); + CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb); + double len = ub - lb; + CPXchgrngval(cplexEnv(), _prob, 1, &i, &len); + } + } + + void CplexBase::_setRowLowerBound(int i, Value lb) + { + LEMON_ASSERT(lb != INF, "Invalid bound"); + _set_row_bounds(i, lb, CplexBase::_getRowUpperBound(i)); + } + + void CplexBase::_setRowUpperBound(int i, Value ub) + { + + LEMON_ASSERT(ub != -INF, "Invalid bound"); + _set_row_bounds(i, CplexBase::_getRowLowerBound(i), ub); + } + + void CplexBase::_setObjCoeffs(ExprIterator b, ExprIterator e) + { + std::vector indices; + std::vector values; + for(ExprIterator it=b; it!=e; ++it) { + indices.push_back(it->first); + values.push_back(it->second); + } + CPXchgobj(cplexEnv(), _prob, values.size(), + &indices.front(), &values.front()); } - LpCplex::Value LpCplex::_getColLowerBound(int i) const + void CplexBase::_getObjCoeffs(InsertIterator b) const { - LpCplex::Value x; - CPXgetlb (env, lp, &x, i, i); - if (x <= -CPX_INFBOUND) x = -INF; - return x; - } + int num = CPXgetnumcols(cplexEnv(), _prob); + std::vector x(num); - 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); - } + CPXgetobj(cplexEnv(), _prob, &x.front(), 0, num - 1); + for (int i = 0; i < num; ++i) { + if (x[i] != 0.0) { + *b = std::make_pair(i, x[i]); + ++b; } } } -// 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 + void CplexBase::_setObjCoeff(int i, Value obj_coef) { - 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; - } + CPXchgobj(cplexEnv(), _prob, 1, &i, &obj_coef); } - void LpCplex::_setObjCoeff(int i, Value obj_coef) - { - CPXchgcoef(env, lp, -1, i, obj_coef); - } - - LpCplex::Value LpCplex::_getObjCoeff(int i) const + CplexBase::Value CplexBase::_getObjCoeff(int i) const { Value x; - CPXgetcoef(env, lp, -1, i, &x); + CPXgetobj(cplexEnv(), _prob, &x, i, i); return x; } - void LpCplex::_clearObj() - { - for (int i=0;i< CPXgetnumcols(env, lp);++i){ - CPXchgcoef(env, lp, -1, i, 0); + void CplexBase::_setSense(CplexBase::Sense sense) { + switch (sense) { + case MIN: + CPXchgobjsen(cplexEnv(), _prob, CPX_MIN); + break; + case MAX: + CPXchgobjsen(cplexEnv(), _prob, CPX_MAX); + break; } + } + CplexBase::Sense CplexBase::_getSense() const { + switch (CPXgetobjsen(cplexEnv(), _prob)) { + case CPX_MIN: + return MIN; + case CPX_MAX: + return MAX; + default: + LEMON_ASSERT(false, "Invalid sense"); + return CplexBase::Sense(); + } } + + void CplexBase::_clear() { + CPXfreeprob(cplexEnv(),&_prob); + int status; + _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem"); + rows.clear(); + cols.clear(); + } + + // LpCplex members + + LpCplex::LpCplex() + : LpBase(), CplexBase(), LpSolver() {} + + LpCplex::LpCplex(const CplexEnv& env) + : LpBase(), CplexBase(env), LpSolver() {} + + LpCplex::LpCplex(const LpCplex& other) + : LpBase(), CplexBase(other), LpSolver() {} + + LpCplex::~LpCplex() {} + + LpCplex* LpCplex::_newSolver() const { return new LpCplex; } + LpCplex* LpCplex::_cloneSolver() const {return new LpCplex(*this); } + + const char* LpCplex::_solverName() const { return "LpCplex"; } + + void LpCplex::_clear_temporals() { + _col_status.clear(); + _row_status.clear(); + _primal_ray.clear(); + _dual_ray.clear(); + } + // 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 @@ -377,32 +472,24 @@ // 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); + LpCplex::SolveExitStatus LpCplex::convertStatus(int status) { #if CPX_VERSION >= 800 - if (status) - { + if (status == 0) { + switch (CPXgetstat(cplexEnv(), _prob)) { + case CPX_STAT_OPTIMAL: + case CPX_STAT_INFEASIBLE: + case CPX_STAT_UNBOUNDED: + return SOLVED; + default: + return UNSOLVED; + } + } else { 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){ + if (status == 0) { //We want to exclude some cases - switch (CPXgetstat(env, lp)){ + switch (CPXgetstat(cplexEnv(), _prob)) { case CPX_OBJ_LIM: case CPX_IT_LIM_FEAS: case CPX_IT_LIM_INFEAS: @@ -412,115 +499,179 @@ default: return SOLVED; } - } - else{ + } else { return UNSOLVED; } #endif } - LpCplex::Value LpCplex::_getPrimal(int i) const - { + LpCplex::SolveExitStatus LpCplex::_solve() { + _clear_temporals(); + return convertStatus(CPXlpopt(cplexEnv(), _prob)); + } + + LpCplex::SolveExitStatus LpCplex::solvePrimal() { + _clear_temporals(); + return convertStatus(CPXprimopt(cplexEnv(), _prob)); + } + + LpCplex::SolveExitStatus LpCplex::solveDual() { + _clear_temporals(); + return convertStatus(CPXdualopt(cplexEnv(), _prob)); + } + + LpCplex::SolveExitStatus LpCplex::solveBarrier() { + _clear_temporals(); + return convertStatus(CPXbaropt(cplexEnv(), _prob)); + } + + LpCplex::Value LpCplex::_getPrimal(int i) const { Value x; - CPXgetx(env, lp, &x, i, i); + CPXgetx(cplexEnv(), _prob, &x, i, i); return x; } - LpCplex::Value LpCplex::_getDual(int i) const - { + LpCplex::Value LpCplex::_getDual(int i) const { Value y; - CPXgetpi(env, lp, &y, i, i); + CPXgetpi(cplexEnv(), _prob, &y, i, i); return y; } - LpCplex::Value LpCplex::_getPrimalValue() const - { + 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); + CPXgetobjval(cplexEnv(), _prob, &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); + + LpCplex::VarStatus LpCplex::_getColStatus(int i) const { + if (_col_status.empty()) { + _col_status.resize(CPXgetnumcols(cplexEnv(), _prob)); + CPXgetbase(cplexEnv(), _prob, &_col_status.front(), 0); + } + switch (_col_status[i]) { + case CPX_BASIC: + return BASIC; + case CPX_FREE_SUPER: + return FREE; + case CPX_AT_LOWER: + return LOWER; + case CPX_AT_UPPER: + return UPPER; + default: + LEMON_ASSERT(false, "Wrong column status"); + return LpCplex::VarStatus(); + } } -//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. + LpCplex::VarStatus LpCplex::_getRowStatus(int i) const { + if (_row_status.empty()) { + _row_status.resize(CPXgetnumrows(cplexEnv(), _prob)); + CPXgetbase(cplexEnv(), _prob, 0, &_row_status.front()); + } + switch (_row_status[i]) { + case CPX_BASIC: + return BASIC; + case CPX_AT_LOWER: + { + char s; + CPXgetsense(cplexEnv(), _prob, &s, i, i); + return s != 'L' ? LOWER : UPPER; + } + case CPX_AT_UPPER: + return UPPER; + default: + LEMON_ASSERT(false, "Wrong row status"); + return LpCplex::VarStatus(); + } + } -// 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 + LpCplex::Value LpCplex::_getPrimalRay(int i) const { + if (_primal_ray.empty()) { + _primal_ray.resize(CPXgetnumcols(cplexEnv(), _prob)); + CPXgetray(cplexEnv(), _prob, &_primal_ray.front()); + } + return _primal_ray[i]; + } -//Some more interesting stuff: + LpCplex::Value LpCplex::_getDualRay(int i) const { + if (_dual_ray.empty()) { -// 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 + } + return _dual_ray[i]; + } + + //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_PROBMETHOD 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){ + void statusSwitch(CPXENVptr cplexEnv(),int& stat){ int lpmethod; - CPXgetintparam (env,CPX_PARAM_LPMETHOD,&lpmethod); + CPXgetintparam (cplexEnv(),CPX_PARAM_PROBMETHOD,&lpmethod); if (lpmethod==2){ if (stat==CPX_UNBOUNDED){ stat=CPX_INFEASIBLE; @@ -535,8 +686,213 @@ void statusSwitch(CPXENVptr,int&){} #endif - LpCplex::SolutionStatus LpCplex::_getPrimalStatus() const - { + LpCplex::ProblemType LpCplex::_getPrimalType() 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(cplexEnv(), _prob); +#if CPX_VERSION >= 800 + switch (stat) + { + case CPX_STAT_OPTIMAL: + return OPTIMAL; + case CPX_STAT_UNBOUNDED: + return UNBOUNDED; + case CPX_STAT_INFEASIBLE: + return INFEASIBLE; + default: + return UNDEFINED; + } +#else + statusSwitch(cplexEnv(),stat); + //CPXgetstat(cplexEnv(), _prob); + //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 UNBOUNDED; + 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 UNBOUNDED;//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::ProblemType LpCplex::_getDualType() const { + int stat = CPXgetstat(cplexEnv(), _prob); +#if CPX_VERSION >= 800 + switch (stat) { + case CPX_STAT_OPTIMAL: + return OPTIMAL; + case CPX_STAT_UNBOUNDED: + return INFEASIBLE; + default: + return UNDEFINED; + } +#else + statusSwitch(cplexEnv(),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 + } + + // MipCplex members + + MipCplex::MipCplex() + : LpBase(), CplexBase(), MipSolver() { + +#if CPX_VERSION < 800 + CPXchgprobtype(cplexEnv(), _prob, CPXPROB_MIP); +#else + CPXchgprobtype(cplexEnv(), _prob, CPXPROB_MILP); +#endif + } + + MipCplex::MipCplex(const CplexEnv& env) + : LpBase(), CplexBase(env), MipSolver() { + +#if CPX_VERSION < 800 + CPXchgprobtype(cplexEnv(), _prob, CPXPROB_MIP); +#else + CPXchgprobtype(cplexEnv(), _prob, CPXPROB_MILP); +#endif + + } + + MipCplex::MipCplex(const MipCplex& other) + : LpBase(), CplexBase(other), MipSolver() {} + + MipCplex::~MipCplex() {} + + MipCplex* MipCplex::_newSolver() const { return new MipCplex; } + MipCplex* MipCplex::_cloneSolver() const {return new MipCplex(*this); } + + const char* MipCplex::_solverName() const { return "MipCplex"; } + + void MipCplex::_setColType(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. + + switch (col_type){ + case INTEGER: { + const char t = 'I'; + CPXchgctype (cplexEnv(), _prob, 1, &i, &t); + } break; + case REAL: { + const char t = 'C'; + CPXchgctype (cplexEnv(), _prob, 1, &i, &t); + } break; + default: + break; + } + } + + MipCplex::ColTypes MipCplex::_getColType(int i) const { + char t; + CPXgetctype (cplexEnv(), _prob, &t, i, i); + switch (t) { + case 'I': + return INTEGER; + case 'C': + return REAL; + default: + LEMON_ASSERT(false, "Invalid column type"); + return ColTypes(); + } + + } + + MipCplex::SolveExitStatus MipCplex::_solve() { + int status; + status = CPXmipopt (cplexEnv(), _prob); + if (status==0) + return SOLVED; + else + return UNSOLVED; + + } + + + MipCplex::ProblemType MipCplex::_getType() const { + + int stat = CPXgetstat(cplexEnv(), _prob); + + //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 UNBOUNDED; + case CPXMIP_INFEASIBLE: + return INFEASIBLE; + default: + return UNDEFINED; + } //Unboundedness not treated well: the following is from cplex 9.0 doc // About Unboundedness @@ -551,148 +907,18 @@ // 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 + MipCplex::Value MipCplex::_getSol(int i) const { + Value x; + CPXgetmipx(cplexEnv(), _prob, &x, i, i); + return x; } - 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; + MipCplex::Value MipCplex::_getSolValue() const { + Value objval; + CPXgetmipobjval(cplexEnv(), _prob, &objval); + return objval; } } //namespace lemon diff --git a/lemon/lp_cplex.h b/lemon/lp_cplex.h --- a/lemon/lp_cplex.h +++ b/lemon/lp_cplex.h @@ -29,84 +29,227 @@ namespace lemon { - - /// \brief Interface for the CPLEX solver + /// \brief Reference counted wrapper around cpxenv pointer /// - /// This class implements an interface for the CPLEX LP solver. - class LpCplex :virtual public LpSolverBase { + /// The cplex uses environment object which is responsible for + /// checking the proper license usage. This class provides a simple + /// interface for share the environment object between different + /// problems. + class CplexEnv { + friend class CplexBase; + private: + cpxenv* _env; + mutable int* _cnt; public: - typedef LpSolverBase Parent; + /// \brief This exception is thrown when the license check is not + /// sufficient + class LicenseError : public Exception { + friend class CplexEnv; + private: - /// \e - int status; - cpxenv* env; - cpxlp* lp; + LicenseError(int status); + char _message[510]; + public: - /// \e - LpCplex(); - /// \e - LpCplex(const LpCplex&); - /// \e - ~LpCplex(); + /// The short error message + virtual const char* what() const throw() { + return _message; + } + }; + + /// Constructor + CplexEnv(); + /// Shallow copy constructor + CplexEnv(const CplexEnv&); + /// Shallow assignement + CplexEnv& operator=(const CplexEnv&); + /// Destructor + virtual ~CplexEnv(); protected: - virtual LpSolverBase* _newLp(); - virtual LpSolverBase* _copyLp(); + cpxenv* cplexEnv() { return _env; } + const cpxenv* cplexEnv() const { return _env; } + }; + + /// \brief Base interface for the CPLEX LP and MIP solver + /// + /// This class implements the common interface of the CPLEX LP and + /// MIP solvers. + /// \ingroup lp_group + class CplexBase : virtual public LpBase { + protected: + + CplexEnv _env; + cpxlp* _prob; + + CplexBase(); + CplexBase(const CplexEnv&); + CplexBase(const CplexBase &); + virtual ~CplexBase(); 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 void _eraseColId(int i); + virtual void _eraseRowId(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 _getRowName(int row, std::string& name) const; + virtual void _setRowName(int row, const std::string& name); + virtual int _rowByName(const std::string& name) const; + + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e); + virtual void _getRowCoeffs(int i, InsertIterator b) const; + + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e); + virtual void _getColCoeffs(int i, InsertIterator 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; + private: + void _set_row_bounds(int i, Value lb, Value ub); + protected: + + virtual void _setRowLowerBound(int i, Value value); + virtual Value _getRowLowerBound(int i) const; + + virtual void _setRowUpperBound(int i, Value value); + virtual Value _getRowUpperBound(int i) const; + + virtual void _setObjCoeffs(ExprIterator b, ExprIterator e); + virtual void _getObjCoeffs(InsertIterator b) const; + virtual void _setObjCoeff(int i, Value obj_coef); virtual Value _getObjCoeff(int i) const; - virtual void _clearObj(); + virtual void _setSense(Sense sense); + virtual Sense _getSense() const; + + virtual void _clear(); + + public: + + /// Returns the used \c CplexEnv instance + const CplexEnv& env() const { return _env; } + /// + const cpxenv* cplexEnv() const { return _env.cplexEnv(); } + + cpxlp* cplexLp() { return _prob; } + const cpxlp* cplexLp() const { return _prob; } + + }; + + /// \brief Interface for the CPLEX LP solver + /// + /// This class implements an interface for the CPLEX LP solver. + ///\ingroup lp_group + class LpCplex : public CplexBase, public LpSolver { + public: + /// \e + LpCplex(); + /// \e + LpCplex(const CplexEnv&); + /// \e + LpCplex(const LpCplex&); + /// \e + virtual ~LpCplex(); + + private: + + // these values cannot retrieved element by element + mutable std::vector _col_status; + mutable std::vector _row_status; + + mutable std::vector _primal_ray; + mutable std::vector _dual_ray; + + void _clear_temporals(); + + SolveExitStatus convertStatus(int status); + + protected: + + virtual LpCplex* _cloneSolver() const; + virtual LpCplex* _newSolver() const; + + virtual const char* _solverName() const; 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 VarStatus _getColStatus(int i) const; + virtual VarStatus _getRowStatus(int i) const; + virtual Value _getPrimalRay(int i) const; + virtual Value _getDualRay(int i) const; - virtual void _setMax(); - virtual void _setMin(); - - virtual bool _isMax() const; + virtual ProblemType _getPrimalType() const; + virtual ProblemType _getDualType() const; public: - cpxenv* cplexEnv() { return env; } - cpxlp* cplexLp() { return lp; } + /// Solve with primal simplex method + SolveExitStatus solvePrimal(); + + /// Solve with dual simplex method + SolveExitStatus solveDual(); + + /// Solve with barrier method + SolveExitStatus solveBarrier(); }; + + /// \brief Interface for the CPLEX MIP solver + /// + /// This class implements an interface for the CPLEX MIP solver. + ///\ingroup lp_group + class MipCplex : public CplexBase, public MipSolver { + public: + /// \e + MipCplex(); + /// \e + MipCplex(const CplexEnv&); + /// \e + MipCplex(const MipCplex&); + /// \e + virtual ~MipCplex(); + + protected: + + virtual MipCplex* _cloneSolver() const; + virtual MipCplex* _newSolver() const; + + virtual const char* _solverName() const; + + virtual ColTypes _getColType(int col) const; + virtual void _setColType(int col, ColTypes col_type); + + virtual SolveExitStatus _solve(); + virtual ProblemType _getType() const; + virtual Value _getSol(int i) const; + virtual Value _getSolValue() const; + + }; + } //END OF NAMESPACE LEMON #endif //LEMON_LP_CPLEX_H diff --git a/lemon/lp_glpk.cc b/lemon/lp_glpk.cc --- a/lemon/lp_glpk.cc +++ b/lemon/lp_glpk.cc @@ -17,628 +17,936 @@ */ ///\file -///\brief Implementation of the LEMON-GLPK lp solver interface. +///\brief Implementation of the LEMON GLPK LP and MIP solver interface. #include -//#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 +#include 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); + // GlpkBase members + + GlpkBase::GlpkBase() : LpBase() { + lp = glp_create_prob(); + glp_create_index(lp); } - 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; + GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() { + lp = glp_create_prob(); + glp_copy_prob(lp, other.lp, GLP_ON); + glp_create_index(lp); + rows = other.rows; + cols = other.cols; } - LpGlpk::~LpGlpk() { - LEMON_glp(delete_prob)(lp); + GlpkBase::~GlpkBase() { + 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; + int GlpkBase::_addCol() { + int i = glp_add_cols(lp, 1); + glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0); 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; + int GlpkBase::_addRow() { + int i = glp_add_rows(lp, 1); + glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0); return i; } - - void LpGlpk::_eraseCol(int i) { + void GlpkBase::_eraseCol(int i) { int ca[2]; - ca[1]=i; - LEMON_glp(del_cols)(lp, 1, ca); - solved = false; + ca[1] = i; + glp_del_cols(lp, 1, ca); } - void LpGlpk::_eraseRow(int i) { + void GlpkBase::_eraseRow(int i) { int ra[2]; - ra[1]=i; - LEMON_glp(del_rows)(lp, 1, ra); - solved = false; + ra[1] = i; + glp_del_rows(lp, 1, ra); } - void LpGlpk::_getColName(int c, std::string & name) const - { - - const char *n = LEMON_glp(get_col_name)(lp,c); - name = n?n:""; + void GlpkBase::_eraseColId(int i) { + cols.eraseIndex(i); + cols.shiftIndices(i); } + void GlpkBase::_eraseRowId(int i) { + rows.eraseIndex(i); + rows.shiftIndices(i); + } - void LpGlpk::_setColName(int c, const std::string & name) - { - LEMON_glp(set_col_name)(lp,c,const_cast(name.c_str())); + void GlpkBase::_getColName(int c, std::string& name) const { + const char *str = glp_get_col_name(lp, c); + if (str) name = str; + else name.clear(); + } + + void GlpkBase::_setColName(int c, const std::string & name) { + 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())); + int GlpkBase::_colByName(const std::string& name) const { + int k = glp_find_col(lp, const_cast(name.c_str())); return k > 0 ? k : -1; } + void GlpkBase::_getRowName(int r, std::string& name) const { + const char *str = glp_get_row_name(lp, r); + if (str) name = str; + else name.clear(); + } - void LpGlpk::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e) - { - std::vector indices; + void GlpkBase::_setRowName(int r, const std::string & name) { + glp_set_row_name(lp, r, const_cast(name.c_str())); + + } + + int GlpkBase::_rowByName(const std::string& name) const { + int k = glp_find_row(lp, const_cast(name.c_str())); + return k > 0 ? k : -1; + } + + void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) { + std::vector indexes; std::vector values; - indices.push_back(0); + indexes.push_back(0); values.push_back(0); - for(ConstRowIterator it=b; it!=e; ++it) { - indices.push_back(it->first); + for(ExprIterator it = b; it != e; ++it) { + indexes.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; + glp_set_mat_row(lp, i, values.size() - 1, + &indexes.front(), &values.front()); } - void LpGlpk::_getRowCoeffs(int ix, RowIterator b) const - { - int length = LEMON_glp(get_mat_row)(lp, ix, 0, 0); + void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const { + int length = glp_get_mat_row(lp, ix, 0, 0); - std::vector indices(length + 1); + std::vector indexes(length + 1); std::vector values(length + 1); - LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]); + glp_get_mat_row(lp, ix, &indexes.front(), &values.front()); for (int i = 1; i <= length; ++i) { - *b = std::make_pair(indices[i], values[i]); + *b = std::make_pair(indexes[i], values[i]); ++b; } } - void LpGlpk::_setColCoeffs(int ix, ConstColIterator b, ConstColIterator e) { + void GlpkBase::_setColCoeffs(int ix, ExprIterator b, + ExprIterator e) { - std::vector indices; + std::vector indexes; std::vector values; - indices.push_back(0); + indexes.push_back(0); values.push_back(0); - for(ConstColIterator it=b; it!=e; ++it) { - indices.push_back(it->first); + for(ExprIterator it = b; it != e; ++it) { + indexes.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; + glp_set_mat_col(lp, ix, values.size() - 1, + &indexes.front(), &values.front()); } - void LpGlpk::_getColCoeffs(int ix, ColIterator b) const - { - int length = LEMON_glp(get_mat_col)(lp, ix, 0, 0); + void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const { + int length = glp_get_mat_col(lp, ix, 0, 0); - std::vector indices(length + 1); + std::vector indexes(length + 1); std::vector values(length + 1); - LEMON_glp(get_mat_col)(lp, ix, &indices[0], &values[0]); + glp_get_mat_col(lp, ix, &indexes.front(), &values.front()); - for (int i = 1; i <= length; ++i) { - *b = std::make_pair(indices[i], values[i]); + for (int i = 1; i <= length; ++i) { + *b = std::make_pair(indexes[i], values[i]); ++b; } } - void LpGlpk::_setCoeff(int ix, int jx, Value value) - { + void GlpkBase::_setCoeff(int ix, int jx, Value value) { - if (LEMON_glp(get_num_cols)(lp) < LEMON_glp(get_num_rows)(lp)) { + if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) { - int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0); + int length = glp_get_mat_row(lp, ix, 0, 0); - std::vector indices(length + 2); + std::vector indexes(length + 2); std::vector values(length + 2); - LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]); + glp_get_mat_row(lp, ix, &indexes.front(), &values.front()); //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; + //array indexes are sorted + bool found = false; + for (int i = 1; i <= length; ++i) { + if (indexes[i] == jx) { + found = true; + values[i] = value; break; } } - if (!found){ + if (!found) { ++length; - indices[length]=jx; - values[length]=value; + indexes[length] = jx; + values[length] = value; } - LEMON_glp(set_mat_row)(lp, ix, length, &indices[0], &values[0]); + glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front()); } else { - int length=LEMON_glp(get_mat_col)(lp, jx, 0, 0); + int length = glp_get_mat_col(lp, jx, 0, 0); - std::vector indices(length + 2); + std::vector indexes(length + 2); std::vector values(length + 2); - LEMON_glp(get_mat_col)(lp, jx, &indices[0], &values[0]); + glp_get_mat_col(lp, jx, &indexes.front(), &values.front()); //The following code does not suppose that the elements of the - //array indices are sorted - bool found=false; + //array indexes are sorted + bool found = false; for (int i = 1; i <= length; ++i) { - if (indices[i]==ix){ - found=true; - values[i]=value; + if (indexes[i] == ix) { + found = true; + values[i] = value; break; } } - if (!found){ + if (!found) { ++length; - indices[length]=ix; - values[length]=value; + indexes[length] = ix; + values[length] = value; } - LEMON_glp(set_mat_col)(lp, jx, length, &indices[0], &values[0]); + glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front()); } - solved = false; } - LpGlpk::Value LpGlpk::_getCoeff(int ix, int jx) const - { + GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const { - int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0); + int length = glp_get_mat_row(lp, ix, 0, 0); - std::vector indices(length + 1); + std::vector indexes(length + 1); std::vector values(length + 1); - LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]); + glp_get_mat_row(lp, ix, &indexes.front(), &values.front()); - //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){ + for (int i = 1; i <= length; ++i) { + if (indexes[i] == jx) { return values[i]; } } + return 0; + } + + void GlpkBase::_setColLowerBound(int i, Value lo) { + LEMON_ASSERT(lo != INF, "Invalid bound"); + + int b = glp_get_col_type(lp, i); + double up = glp_get_col_ub(lp, i); + if (lo == -INF) { + switch (b) { + case GLP_FR: + case GLP_LO: + glp_set_col_bnds(lp, i, GLP_FR, lo, up); + break; + case GLP_UP: + break; + case GLP_DB: + case GLP_FX: + glp_set_col_bnds(lp, i, GLP_UP, lo, up); + break; + default: + break; + } + } else { + switch (b) { + case GLP_FR: + case GLP_LO: + glp_set_col_bnds(lp, i, GLP_LO, lo, up); + break; + case GLP_UP: + case GLP_DB: + case GLP_FX: + if (lo == up) + glp_set_col_bnds(lp, i, GLP_FX, lo, up); + else + glp_set_col_bnds(lp, i, GLP_DB, lo, up); + break; + default: + break; + } + } + } + + GlpkBase::Value GlpkBase::_getColLowerBound(int i) const { + int b = glp_get_col_type(lp, i); + switch (b) { + case GLP_LO: + case GLP_DB: + case GLP_FX: + return glp_get_col_lb(lp, i); + default: + return -INF; + } + } + + void GlpkBase::_setColUpperBound(int i, Value up) { + LEMON_ASSERT(up != -INF, "Invalid bound"); + + int b = glp_get_col_type(lp, i); + double lo = glp_get_col_lb(lp, i); + if (up == INF) { + switch (b) { + case GLP_FR: + case GLP_LO: + break; + case GLP_UP: + glp_set_col_bnds(lp, i, GLP_FR, lo, up); + break; + case GLP_DB: + case GLP_FX: + glp_set_col_bnds(lp, i, GLP_LO, lo, up); + break; + default: + break; + } + } else { + switch (b) { + case GLP_FR: + glp_set_col_bnds(lp, i, GLP_UP, lo, up); + break; + case GLP_UP: + glp_set_col_bnds(lp, i, GLP_UP, lo, up); + break; + case GLP_LO: + case GLP_DB: + case GLP_FX: + if (lo == up) + glp_set_col_bnds(lp, i, GLP_FX, lo, up); + else + glp_set_col_bnds(lp, i, GLP_DB, lo, up); + break; + default: + break; + } + } } - - 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) { + GlpkBase::Value GlpkBase::_getColUpperBound(int i) const { + int b = glp_get_col_type(lp, i); 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: ; + case GLP_UP: + case GLP_DB: + case GLP_FX: + return 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 - } + void GlpkBase::_setRowLowerBound(int i, Value lo) { + LEMON_ASSERT(lo != INF, "Invalid bound"); - if (lb == -INF){ - if (ub == INF){ - LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FR), lb, ub); + int b = glp_get_row_type(lp, i); + double up = glp_get_row_ub(lp, i); + if (lo == -INF) { + switch (b) { + case GLP_FR: + case GLP_LO: + glp_set_row_bnds(lp, i, GLP_FR, lo, up); + break; + case GLP_UP: + break; + case GLP_DB: + case GLP_FX: + glp_set_row_bnds(lp, i, GLP_UP, lo, up); + break; + default: + break; } - else{ - LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(UP), lb, ub); + } else { + switch (b) { + case GLP_FR: + case GLP_LO: + glp_set_row_bnds(lp, i, GLP_LO, lo, up); + break; + case GLP_UP: + case GLP_DB: + case GLP_FX: + if (lo == up) + glp_set_row_bnds(lp, i, GLP_FX, lo, up); + else + glp_set_row_bnds(lp, i, GLP_DB, lo, up); + break; + default: + break; } } - 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; + GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const { + int b = glp_get_row_type(lp, i); + switch (b) { + case GLP_LO: + case GLP_DB: + case GLP_FX: + return glp_get_row_lb(lp, i); default: - return UNSOLVED; + return -INF; } } - LpGlpk::Value LpGlpk::_getPrimal(int i) const - { - return LEMON_glp(get_col_prim)(lp,i); - } + void GlpkBase::_setRowUpperBound(int i, Value up) { + LEMON_ASSERT(up != -INF, "Invalid bound"); - 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 + int b = glp_get_row_type(lp, i); + double lo = glp_get_row_lb(lp, i); + if (up == INF) { + switch (b) { + case GLP_FR: + case GLP_LO: + break; + case GLP_UP: + glp_set_row_bnds(lp, i, GLP_FR, lo, up); + break; + case GLP_DB: + case GLP_FX: + glp_set_row_bnds(lp, i, GLP_LO, lo, up); + break; + default: + break; + } + } else { + switch (b) { + case GLP_FR: + glp_set_row_bnds(lp, i, GLP_UP, lo, up); + break; + case GLP_UP: + glp_set_row_bnds(lp, i, GLP_UP, lo, up); + break; + case GLP_LO: + case GLP_DB: + case GLP_FX: + if (lo == up) + glp_set_row_bnds(lp, i, GLP_FX, lo, up); + else + glp_set_row_bnds(lp, i, GLP_DB, lo, up); + break; + default: + break; + } } } - 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; - } + GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const { + int b = glp_get_row_type(lp, i); + switch (b) { + case GLP_UP: + case GLP_DB: + case GLP_FX: + return glp_get_row_ub(lp, i); default: - return UNDEFINED; //to avoid gcc warning - //FIXME error + return INF; } } - 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 GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) { + for (int i = 1; i <= glp_get_num_cols(lp); ++i) { + glp_set_obj_coef(lp, i, 0.0); + } + for (ExprIterator it = b; it != e; ++it) { + glp_set_obj_coef(lp, it->first, it->second); + } } - void LpGlpk::_setMax() - { - solved = false; - LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MAX)); + void GlpkBase::_getObjCoeffs(InsertIterator b) const { + for (int i = 1; i <= glp_get_num_cols(lp); ++i) { + Value val = glp_get_obj_coef(lp, i); + if (val != 0.0) { + *b = std::make_pair(i, val); + ++b; + } + } } - void LpGlpk::_setMin() - { - solved = false; - LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MIN)); + void GlpkBase::_setObjCoeff(int i, Value obj_coef) { + //i = 0 means the constant term (shift) + glp_set_obj_coef(lp, i, obj_coef); } - bool LpGlpk::_isMax() const - { - return (LEMON_glp(get_obj_dir)(lp)==LEMON_GLP(MAX)); + GlpkBase::Value GlpkBase::_getObjCoeff(int i) const { + //i = 0 means the constant term (shift) + return glp_get_obj_coef(lp, i); } - - - void LpGlpk::messageLevel(int m) - { - LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_MSGLEV), m); + void GlpkBase::_setSense(GlpkBase::Sense sense) { + switch (sense) { + case MIN: + glp_set_obj_dir(lp, GLP_MIN); + break; + case MAX: + glp_set_obj_dir(lp, GLP_MAX); + break; + } } - void LpGlpk::presolver(bool b) - { - LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_PRESOL), b); + GlpkBase::Sense GlpkBase::_getSense() const { + switch(glp_get_obj_dir(lp)) { + case GLP_MIN: + return MIN; + case GLP_MAX: + return MAX; + default: + LEMON_ASSERT(false, "Wrong sense"); + return GlpkBase::Sense(); + } } + void GlpkBase::_clear() { + glp_erase_prob(lp); + rows.clear(); + cols.clear(); + } + + // LpGlpk members + + LpGlpk::LpGlpk() + : LpBase(), GlpkBase(), LpSolver() { + messageLevel(MESSAGE_NO_OUTPUT); + } + + LpGlpk::LpGlpk(const LpGlpk& other) + : LpBase(other), GlpkBase(other), LpSolver(other) { + messageLevel(MESSAGE_NO_OUTPUT); + } + + LpGlpk* LpGlpk::_newSolver() const { return new LpGlpk; } + LpGlpk* LpGlpk::_cloneSolver() const { return new LpGlpk(*this); } + + const char* LpGlpk::_solverName() const { return "LpGlpk"; } + + void LpGlpk::_clear_temporals() { + _primal_ray.clear(); + _dual_ray.clear(); + } + + LpGlpk::SolveExitStatus LpGlpk::_solve() { + return solvePrimal(); + } + + LpGlpk::SolveExitStatus LpGlpk::solvePrimal() { + _clear_temporals(); + + glp_smcp smcp; + glp_init_smcp(&smcp); + + switch (_message_level) { + case MESSAGE_NO_OUTPUT: + smcp.msg_lev = GLP_MSG_OFF; + break; + case MESSAGE_ERROR_MESSAGE: + smcp.msg_lev = GLP_MSG_ERR; + break; + case MESSAGE_NORMAL_OUTPUT: + smcp.msg_lev = GLP_MSG_ON; + break; + case MESSAGE_FULL_OUTPUT: + smcp.msg_lev = GLP_MSG_ALL; + break; + } + + if (glp_simplex(lp, &smcp) != 0) return UNSOLVED; + return SOLVED; + } + + LpGlpk::SolveExitStatus LpGlpk::solveDual() { + _clear_temporals(); + + glp_smcp smcp; + glp_init_smcp(&smcp); + + switch (_message_level) { + case MESSAGE_NO_OUTPUT: + smcp.msg_lev = GLP_MSG_OFF; + break; + case MESSAGE_ERROR_MESSAGE: + smcp.msg_lev = GLP_MSG_ERR; + break; + case MESSAGE_NORMAL_OUTPUT: + smcp.msg_lev = GLP_MSG_ON; + break; + case MESSAGE_FULL_OUTPUT: + smcp.msg_lev = GLP_MSG_ALL; + break; + } + smcp.meth = GLP_DUAL; + + if (glp_simplex(lp, &smcp) != 0) return UNSOLVED; + return SOLVED; + } + + LpGlpk::Value LpGlpk::_getPrimal(int i) const { + return glp_get_col_prim(lp, i); + } + + LpGlpk::Value LpGlpk::_getDual(int i) const { + return glp_get_row_dual(lp, i); + } + + LpGlpk::Value LpGlpk::_getPrimalValue() const { + return glp_get_obj_val(lp); + } + + LpGlpk::VarStatus LpGlpk::_getColStatus(int i) const { + switch (glp_get_col_stat(lp, i)) { + case GLP_BS: + return BASIC; + case GLP_UP: + return UPPER; + case GLP_LO: + return LOWER; + case GLP_NF: + return FREE; + case GLP_NS: + return FIXED; + default: + LEMON_ASSERT(false, "Wrong column status"); + return LpGlpk::VarStatus(); + } + } + + LpGlpk::VarStatus LpGlpk::_getRowStatus(int i) const { + switch (glp_get_row_stat(lp, i)) { + case GLP_BS: + return BASIC; + case GLP_UP: + return UPPER; + case GLP_LO: + return LOWER; + case GLP_NF: + return FREE; + case GLP_NS: + return FIXED; + default: + LEMON_ASSERT(false, "Wrong row status"); + return LpGlpk::VarStatus(); + } + } + + LpGlpk::Value LpGlpk::_getPrimalRay(int i) const { + if (_primal_ray.empty()) { + int row_num = glp_get_num_rows(lp); + int col_num = glp_get_num_cols(lp); + + _primal_ray.resize(col_num + 1, 0.0); + + int index = glp_get_unbnd_ray(lp); + if (index != 0) { + // The primal ray is found in primal simplex second phase + LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) : + glp_get_col_stat(lp, index - row_num)) != GLP_BS, + "Wrong primal ray"); + + bool negate = glp_get_obj_dir(lp) == GLP_MAX; + + if (index > row_num) { + _primal_ray[index - row_num] = 1.0; + if (glp_get_col_dual(lp, index - row_num) > 0) { + negate = !negate; + } + } else { + if (glp_get_row_dual(lp, index) > 0) { + negate = !negate; + } + } + + std::vector ray_indexes(row_num + 1); + std::vector ray_values(row_num + 1); + int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(), + &ray_values.front()); + + for (int i = 1; i <= ray_length; ++i) { + if (ray_indexes[i] > row_num) { + _primal_ray[ray_indexes[i] - row_num] = ray_values[i]; + } + } + + if (negate) { + for (int i = 1; i <= col_num; ++i) { + _primal_ray[i] = - _primal_ray[i]; + } + } + } else { + for (int i = 1; i <= col_num; ++i) { + _primal_ray[i] = glp_get_col_prim(lp, i); + } + } + } + return _primal_ray[i]; + } + + LpGlpk::Value LpGlpk::_getDualRay(int i) const { + if (_dual_ray.empty()) { + int row_num = glp_get_num_rows(lp); + + _dual_ray.resize(row_num + 1, 0.0); + + int index = glp_get_unbnd_ray(lp); + if (index != 0) { + // The dual ray is found in dual simplex second phase + LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) : + glp_get_col_stat(lp, index - row_num)) == GLP_BS, + + "Wrong dual ray"); + + int idx; + bool negate = false; + + if (index > row_num) { + idx = glp_get_col_bind(lp, index - row_num); + if (glp_get_col_prim(lp, index - row_num) > + glp_get_col_ub(lp, index - row_num)) { + negate = true; + } + } else { + idx = glp_get_row_bind(lp, index); + if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) { + negate = true; + } + } + + _dual_ray[idx] = negate ? - 1.0 : 1.0; + + glp_btran(lp, &_dual_ray.front()); + } else { + double eps = 1e-7; + // The dual ray is found in primal simplex first phase + // We assume that the glpk minimizes the slack to get feasible solution + for (int i = 1; i <= row_num; ++i) { + int index = glp_get_bhead(lp, i); + if (index <= row_num) { + double res = glp_get_row_prim(lp, index); + if (res > glp_get_row_ub(lp, index) + eps) { + _dual_ray[i] = -1; + } else if (res < glp_get_row_lb(lp, index) - eps) { + _dual_ray[i] = 1; + } else { + _dual_ray[i] = 0; + } + _dual_ray[i] *= glp_get_rii(lp, index); + } else { + double res = glp_get_col_prim(lp, index - row_num); + if (res > glp_get_col_ub(lp, index - row_num) + eps) { + _dual_ray[i] = -1; + } else if (res < glp_get_col_lb(lp, index - row_num) - eps) { + _dual_ray[i] = 1; + } else { + _dual_ray[i] = 0; + } + _dual_ray[i] /= glp_get_sjj(lp, index - row_num); + } + } + + glp_btran(lp, &_dual_ray.front()); + + for (int i = 1; i <= row_num; ++i) { + _dual_ray[i] /= glp_get_rii(lp, i); + } + } + } + return _dual_ray[i]; + } + + LpGlpk::ProblemType LpGlpk::_getPrimalType() const { + if (glp_get_status(lp) == GLP_OPT) + return OPTIMAL; + switch (glp_get_prim_stat(lp)) { + case GLP_UNDEF: + return UNDEFINED; + case GLP_FEAS: + case GLP_INFEAS: + if (glp_get_dual_stat(lp) == GLP_NOFEAS) { + return UNBOUNDED; + } else { + return UNDEFINED; + } + case GLP_NOFEAS: + return INFEASIBLE; + default: + LEMON_ASSERT(false, "Wrong primal type"); + return LpGlpk::ProblemType(); + } + } + + LpGlpk::ProblemType LpGlpk::_getDualType() const { + if (glp_get_status(lp) == GLP_OPT) + return OPTIMAL; + switch (glp_get_dual_stat(lp)) { + case GLP_UNDEF: + return UNDEFINED; + case GLP_FEAS: + case GLP_INFEAS: + if (glp_get_prim_stat(lp) == GLP_NOFEAS) { + return UNBOUNDED; + } else { + return UNDEFINED; + } + case GLP_NOFEAS: + return INFEASIBLE; + default: + LEMON_ASSERT(false, "Wrong primal type"); + return LpGlpk::ProblemType(); + } + } + + void LpGlpk::presolver(bool b) { + lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0); + } + + void LpGlpk::messageLevel(MessageLevel m) { + _message_level = m; + } + + // MipGlpk members + + MipGlpk::MipGlpk() + : LpBase(), GlpkBase(), MipSolver() { + messageLevel(MESSAGE_NO_OUTPUT); + } + + MipGlpk::MipGlpk(const MipGlpk& other) + : LpBase(), GlpkBase(other), MipSolver() { + messageLevel(MESSAGE_NO_OUTPUT); + } + + void MipGlpk::_setColType(int i, MipGlpk::ColTypes col_type) { + switch (col_type) { + case INTEGER: + glp_set_col_kind(lp, i, GLP_IV); + break; + case REAL: + glp_set_col_kind(lp, i, GLP_CV); + break; + } + } + + MipGlpk::ColTypes MipGlpk::_getColType(int i) const { + switch (glp_get_col_kind(lp, i)) { + case GLP_IV: + case GLP_BV: + return INTEGER; + default: + return REAL; + } + + } + + MipGlpk::SolveExitStatus MipGlpk::_solve() { + glp_smcp smcp; + glp_init_smcp(&smcp); + + switch (_message_level) { + case MESSAGE_NO_OUTPUT: + smcp.msg_lev = GLP_MSG_OFF; + break; + case MESSAGE_ERROR_MESSAGE: + smcp.msg_lev = GLP_MSG_ERR; + break; + case MESSAGE_NORMAL_OUTPUT: + smcp.msg_lev = GLP_MSG_ON; + break; + case MESSAGE_FULL_OUTPUT: + smcp.msg_lev = GLP_MSG_ALL; + break; + } + smcp.meth = GLP_DUAL; + + if (glp_simplex(lp, &smcp) != 0) return UNSOLVED; + if (glp_get_status(lp) != GLP_OPT) return SOLVED; + + glp_iocp iocp; + glp_init_iocp(&iocp); + + switch (_message_level) { + case MESSAGE_NO_OUTPUT: + iocp.msg_lev = GLP_MSG_OFF; + break; + case MESSAGE_ERROR_MESSAGE: + iocp.msg_lev = GLP_MSG_ERR; + break; + case MESSAGE_NORMAL_OUTPUT: + iocp.msg_lev = GLP_MSG_ON; + break; + case MESSAGE_FULL_OUTPUT: + iocp.msg_lev = GLP_MSG_ALL; + break; + } + + if (glp_intopt(lp, &iocp) != 0) return UNSOLVED; + return SOLVED; + } + + + MipGlpk::ProblemType MipGlpk::_getType() const { + switch (glp_get_status(lp)) { + case GLP_OPT: + switch (glp_mip_status(lp)) { + case GLP_UNDEF: + return UNDEFINED; + case GLP_NOFEAS: + return INFEASIBLE; + case GLP_FEAS: + return FEASIBLE; + case GLP_OPT: + return OPTIMAL; + default: + LEMON_ASSERT(false, "Wrong problem type."); + return MipGlpk::ProblemType(); + } + case GLP_NOFEAS: + return INFEASIBLE; + case GLP_INFEAS: + case GLP_FEAS: + if (glp_get_dual_stat(lp) == GLP_NOFEAS) { + return UNBOUNDED; + } else { + return UNDEFINED; + } + default: + LEMON_ASSERT(false, "Wrong problem type."); + return MipGlpk::ProblemType(); + } + } + + MipGlpk::Value MipGlpk::_getSol(int i) const { + return glp_mip_col_val(lp, i); + } + + MipGlpk::Value MipGlpk::_getSolValue() const { + return glp_mip_obj_val(lp); + } + + MipGlpk* MipGlpk::_newSolver() const { return new MipGlpk; } + MipGlpk* MipGlpk::_cloneSolver() const {return new MipGlpk(*this); } + + const char* MipGlpk::_solverName() const { return "MipGlpk"; } + + void MipGlpk::messageLevel(MessageLevel m) { + _message_level = m; + } } //END OF NAMESPACE LEMON diff --git a/lemon/lp_glpk.h b/lemon/lp_glpk.h --- a/lemon/lp_glpk.h +++ b/lemon/lp_glpk.h @@ -35,88 +35,137 @@ namespace lemon { - /// \brief Interface for the GLPK LP solver + /// \brief Base interface for the GLPK LP and MIP solver /// - /// This class implements an interface for the GLPK LP solver. - ///\ingroup lp_group - class LpGlpk : virtual public LpSolverBase { + /// This class implements the common interface of the GLPK LP and MIP solver. + /// \ingroup lp_group + class GlpkBase : virtual public LpBase { protected: typedef glp_prob LPX; glp_prob* lp; - bool solved; - public: - - typedef LpSolverBase Parent; - - LpGlpk(); - LpGlpk(const LpGlpk &); - ~LpGlpk(); + GlpkBase(); + GlpkBase(const GlpkBase&); + virtual ~GlpkBase(); 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 void _eraseColId(int i); + virtual void _eraseRowId(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 _getRowName(int row, std::string& name) const; + virtual void _setRowName(int row, const std::string& name); + virtual int _rowByName(const std::string& name) const; + + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e); + virtual void _getRowCoeffs(int i, InsertIterator b) const; + + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e); + virtual void _getColCoeffs(int i, InsertIterator 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 _setRowLowerBound(int i, Value value); + virtual Value _getRowLowerBound(int i) const; + + virtual void _setRowUpperBound(int i, Value value); + virtual Value _getRowUpperBound(int i) const; + + virtual void _setObjCoeffs(ExprIterator b, ExprIterator e); + virtual void _getObjCoeffs(InsertIterator b) const; + virtual void _setObjCoeff(int i, Value obj_coef); virtual Value _getObjCoeff(int i) const; - virtual void _clearObj(); + + virtual void _setSense(Sense); + virtual Sense _getSense() const; + + virtual void _clear(); + + public: + + ///Pointer to the underlying GLPK data structure. + LPX *lpx() {return lp;} + ///Const pointer to the underlying GLPK data structure. + const LPX *lpx() const {return lp;} + + ///Returns the constraint identifier understood by GLPK. + int lpxRow(Row r) const { return rows(id(r)); } + + ///Returns the variable identifier understood by GLPK. + int lpxCol(Col c) const { return cols(id(c)); } + + }; + + /// \brief Interface for the GLPK LP solver + /// + /// This class implements an interface for the GLPK LP solver. + ///\ingroup lp_group + class LpGlpk : public GlpkBase, public LpSolver { + public: ///\e + LpGlpk(); + ///\e + LpGlpk(const LpGlpk&); + + private: + + mutable std::vector _primal_ray; + mutable std::vector _dual_ray; + + void _clear_temporals(); + + protected: + + virtual LpGlpk* _cloneSolver() const; + virtual LpGlpk* _newSolver() const; + + virtual const char* _solverName() const; + + virtual SolveExitStatus _solve(); + virtual Value _getPrimal(int i) const; + virtual Value _getDual(int i) const; + + virtual Value _getPrimalValue() const; + + virtual VarStatus _getColStatus(int i) const; + virtual VarStatus _getRowStatus(int i) const; + + virtual Value _getPrimalRay(int i) const; + virtual Value _getDualRay(int i) const; ///\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; + virtual ProblemType _getPrimalType() const; + virtual ProblemType _getDualType() 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); + ///Solve with primal simplex + SolveExitStatus solvePrimal(); + + ///Solve with dual simplex + SolveExitStatus solveDual(); + ///Turns on or off the presolver ///Turns on (\c b is \c true) or off (\c b is \c false) the presolver @@ -124,15 +173,86 @@ ///The presolver is off by default. void presolver(bool b); - ///Pointer to the underlying GLPK data structure. - LPX *lpx() {return lp;} + ///Enum for \c messageLevel() parameter + enum MessageLevel { + /// no output (default value) + MESSAGE_NO_OUTPUT = 0, + /// error messages only + MESSAGE_ERROR_MESSAGE = 1, + /// normal output + MESSAGE_NORMAL_OUTPUT = 2, + /// full output (includes informational messages) + MESSAGE_FULL_OUTPUT = 3 + }; - ///Returns the constraint identifier understood by GLPK. - int lpxRow(Row r) { return _lpId(r); } + private: - ///Returns the variable identifier understood by GLPK. - int lpxCol(Col c) { return _lpId(c); } + MessageLevel _message_level; + + 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. + void messageLevel(MessageLevel m); }; + + /// \brief Interface for the GLPK MIP solver + /// + /// This class implements an interface for the GLPK MIP solver. + ///\ingroup lp_group + class MipGlpk : public GlpkBase, public MipSolver { + public: + + ///\e + MipGlpk(); + ///\e + MipGlpk(const MipGlpk&); + + protected: + + virtual MipGlpk* _cloneSolver() const; + virtual MipGlpk* _newSolver() const; + + virtual const char* _solverName() const; + + virtual ColTypes _getColType(int col) const; + virtual void _setColType(int col, ColTypes col_type); + + virtual SolveExitStatus _solve(); + virtual ProblemType _getType() const; + virtual Value _getSol(int i) const; + virtual Value _getSolValue() const; + + ///Enum for \c messageLevel() parameter + enum MessageLevel { + /// no output (default value) + MESSAGE_NO_OUTPUT = 0, + /// error messages only + MESSAGE_ERROR_MESSAGE = 1, + /// normal output + MESSAGE_NORMAL_OUTPUT = 2, + /// full output (includes informational messages) + MESSAGE_FULL_OUTPUT = 3 + }; + + private: + + MessageLevel _message_level; + + 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. + void messageLevel(MessageLevel m); + }; + + } //END OF NAMESPACE LEMON #endif //LEMON_LP_GLPK_H diff --git a/lemon/lp_skeleton.cc b/lemon/lp_skeleton.cc --- a/lemon/lp_skeleton.cc +++ b/lemon/lp_skeleton.cc @@ -22,166 +22,113 @@ ///\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() + int SkeletonSolverBase::_addCol() { return ++col_num; } - int LpSkeleton::_addRow() + int SkeletonSolverBase::_addRow() { return ++row_num; } - void LpSkeleton::_eraseCol(int ) { + void SkeletonSolverBase::_eraseCol(int) {} + void SkeletonSolverBase::_eraseRow(int) {} + + void SkeletonSolverBase::_getColName(int, std::string &) const {} + void SkeletonSolverBase::_setColName(int, const std::string &) {} + int SkeletonSolverBase::_colByName(const std::string&) const { return -1; } + + void SkeletonSolverBase::_getRowName(int, std::string &) const {} + void SkeletonSolverBase::_setRowName(int, const std::string &) {} + int SkeletonSolverBase::_rowByName(const std::string&) const { return -1; } + + void SkeletonSolverBase::_setRowCoeffs(int, ExprIterator, ExprIterator) {} + void SkeletonSolverBase::_getRowCoeffs(int, InsertIterator) const {} + + void SkeletonSolverBase::_setColCoeffs(int, ExprIterator, ExprIterator) {} + void SkeletonSolverBase::_getColCoeffs(int, InsertIterator) const {} + + void SkeletonSolverBase::_setCoeff(int, int, Value) {} + SkeletonSolverBase::Value SkeletonSolverBase::_getCoeff(int, int) const + { return 0; } + + void SkeletonSolverBase::_setColLowerBound(int, Value) {} + SkeletonSolverBase::Value SkeletonSolverBase::_getColLowerBound(int) const + { return 0; } + + void SkeletonSolverBase::_setColUpperBound(int, Value) {} + SkeletonSolverBase::Value SkeletonSolverBase::_getColUpperBound(int) const + { return 0; } + + void SkeletonSolverBase::_setRowLowerBound(int, Value) {} + SkeletonSolverBase::Value SkeletonSolverBase::_getRowLowerBound(int) const + { return 0; } + + void SkeletonSolverBase::_setRowUpperBound(int, Value) {} + SkeletonSolverBase::Value SkeletonSolverBase::_getRowUpperBound(int) const + { return 0; } + + void SkeletonSolverBase::_setObjCoeffs(ExprIterator, ExprIterator) {} + void SkeletonSolverBase::_getObjCoeffs(InsertIterator) const {}; + + void SkeletonSolverBase::_setObjCoeff(int, Value) {} + SkeletonSolverBase::Value SkeletonSolverBase::_getObjCoeff(int) const + { return 0; } + + void SkeletonSolverBase::_setSense(Sense) {} + SkeletonSolverBase::Sense SkeletonSolverBase::_getSense() const + { return MIN; } + + void SkeletonSolverBase::_clear() { + row_num = col_num = 0; } - void LpSkeleton::_eraseRow(int) { - } + LpSkeleton::SolveExitStatus LpSkeleton::_solve() { return SOLVED; } - void LpSkeleton::_getColName(int, std::string &) const { - } + LpSkeleton::Value LpSkeleton::_getPrimal(int) const { return 0; } + LpSkeleton::Value LpSkeleton::_getDual(int) const { return 0; } + LpSkeleton::Value LpSkeleton::_getPrimalValue() const { return 0; } + LpSkeleton::Value LpSkeleton::_getPrimalRay(int) const { return 0; } + LpSkeleton::Value LpSkeleton::_getDualRay(int) const { return 0; } - void LpSkeleton::_setColName(int, const std::string &) { - } + LpSkeleton::ProblemType LpSkeleton::_getPrimalType() const + { return UNDEFINED; } - int LpSkeleton::_colByName(const std::string&) const { return -1; } + LpSkeleton::ProblemType LpSkeleton::_getDualType() const + { return UNDEFINED; } + LpSkeleton::VarStatus LpSkeleton::_getColStatus(int) const + { return BASIC; } - void LpSkeleton::_setRowCoeffs(int, ConstRowIterator, ConstRowIterator) { - } + LpSkeleton::VarStatus LpSkeleton::_getRowStatus(int) const + { return BASIC; } - void LpSkeleton::_getRowCoeffs(int, RowIterator) const { - } + LpSkeleton* LpSkeleton::_newSolver() const + { return static_cast(0); } - void LpSkeleton::_setColCoeffs(int, ConstColIterator, ConstColIterator) { - } + LpSkeleton* LpSkeleton::_cloneSolver() const + { return static_cast(0); } - void LpSkeleton::_getColCoeffs(int, ColIterator) const { - } + const char* LpSkeleton::_solverName() const { return "LpSkeleton"; } - void LpSkeleton::_setCoeff(int, int, Value ) - { - } + MipSkeleton::SolveExitStatus MipSkeleton::_solve() + { return SOLVED; } - LpSkeleton::Value LpSkeleton::_getCoeff(int, int) const - { - return 0; - } + MipSkeleton::Value MipSkeleton::_getSol(int) const { return 0; } + MipSkeleton::Value MipSkeleton::_getSolValue() const { return 0; } + MipSkeleton::ProblemType MipSkeleton::_getType() const + { return UNDEFINED; } - void LpSkeleton::_setColLowerBound(int, Value) - { - } + MipSkeleton* MipSkeleton::_newSolver() const + { return static_cast(0); } - LpSkeleton::Value LpSkeleton::_getColLowerBound(int) const - { - return 0; - } + MipSkeleton* MipSkeleton::_cloneSolver() const + { return static_cast(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; - } + const char* MipSkeleton::_solverName() const { return "MipSkeleton"; } } //namespace lemon diff --git a/lemon/lp_skeleton.h b/lemon/lp_skeleton.h --- a/lemon/lp_skeleton.h +++ b/lemon/lp_skeleton.h @@ -26,15 +26,14 @@ namespace lemon { ///A skeleton class to implement LP solver interfaces - class LpSkeleton :public LpSolverBase { + class SkeletonSolverBase : public virtual LpBase { int col_num,row_num; protected: - ///\e - virtual LpSolverBase* _newLp(); - ///\e - virtual LpSolverBase* _copyLp(); + SkeletonSolverBase() + : col_num(-1), row_num(-1) {} + /// \e virtual int _addCol(); /// \e @@ -43,21 +42,29 @@ virtual void _eraseCol(int i); /// \e virtual void _eraseRow(int i); + /// \e - virtual void _getColName(int col, std::string & name) const; + virtual void _getColName(int col, std::string& name) const; /// \e - virtual void _setColName(int col, const std::string & name); + 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); + virtual void _getRowName(int row, std::string& name) const; /// \e - virtual void _getRowCoeffs(int i, RowIterator b) const; + virtual void _setRowName(int row, const std::string& name); /// \e - virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e); + virtual int _rowByName(const std::string& name) const; + /// \e - virtual void _getColCoeffs(int i, ColIterator b) const; + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e); + /// \e + virtual void _getRowCoeffs(int i, InsertIterator b) const; + /// \e + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e); + /// \e + virtual void _getColCoeffs(int i, InsertIterator b) const; /// Set one element of the coefficient matrix virtual void _setCoeff(int row, int col, Value value); @@ -87,42 +94,103 @@ /// 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 + /// The lower bound of a constraint (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); + /// Value or -\ref INF. + virtual void _setRowLowerBound(int i, Value value); /// \e + /// The lower bound of a constraint (row) is an + /// extended number of type Value, i.e. a finite number of type + /// Value or -\ref INF. + virtual Value _getRowLowerBound(int i) const; - /// 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; + /// 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 or \ref INF. + virtual void _setRowUpperBound(int i, Value value); /// \e + /// The upper bound of a constraint (row) is an + /// extended number of type Value, i.e. a finite number of type + /// Value or \ref INF. + virtual Value _getRowUpperBound(int i) const; /// \e - virtual void _clearObj(); + virtual void _setObjCoeffs(ExprIterator b, ExprIterator e); + /// \e + virtual void _getObjCoeffs(InsertIterator b) const; + /// \e virtual void _setObjCoeff(int i, Value obj_coef); - /// \e virtual Value _getObjCoeff(int i) const; ///\e + virtual void _setSense(Sense); + ///\e + virtual Sense _getSense() const; + + ///\e + virtual void _clear(); + + }; + + /// \brief Interface for a skeleton LP solver + /// + /// This class implements an interface for a skeleton LP solver. + ///\ingroup lp_group + class LpSkeleton : public SkeletonSolverBase, public LpSolver { + public: + LpSkeleton() : SkeletonSolverBase(), LpSolver() {} + + protected: + + ///\e + virtual SolveExitStatus _solve(); + + ///\e + virtual Value _getPrimal(int i) const; + ///\e + virtual Value _getDual(int i) const; + + ///\e + virtual Value _getPrimalValue() const; + + ///\e + virtual Value _getPrimalRay(int i) const; + ///\e + virtual Value _getDualRay(int i) const; + + ///\e + virtual ProblemType _getPrimalType() const; + ///\e + virtual ProblemType _getDualType() const; + + ///\e + virtual VarStatus _getColStatus(int i) const; + ///\e + virtual VarStatus _getRowStatus(int i) const; + + ///\e + virtual LpSkeleton* _newSolver() const; + ///\e + virtual LpSkeleton* _cloneSolver() const; + ///\e + virtual const char* _solverName() const; + + }; + + /// \brief Interface for a skeleton MIP solver + /// + /// This class implements an interface for a skeleton MIP solver. + ///\ingroup lp_group + class MipSkeleton : public SkeletonSolverBase, public MipSolver { + public: + MipSkeleton() : SkeletonSolverBase(), MipSolver() {} + + protected: + ///\e ///\bug Wrong interface /// @@ -132,50 +200,28 @@ ///\bug Wrong interface /// - virtual Value _getPrimal(int i) const; + virtual Value _getSol(int i) const; ///\e ///\bug Wrong interface /// - virtual Value _getDual(int i) const; + virtual Value _getSolValue() const; ///\e ///\bug Wrong interface /// - virtual Value _getPrimalValue() const; + virtual ProblemType _getType() const; ///\e - - ///\bug Wrong interface - /// - virtual SolutionStatus _getPrimalStatus() const; - - ////e - virtual SolutionStatus _getDualStatus() const; - + virtual MipSkeleton* _newSolver() const; ///\e - virtual ProblemTypes _getProblemType() const; + virtual MipSkeleton* _cloneSolver() const; + ///\e + virtual const char* _solverName() 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 diff --git a/lemon/lp_soplex.cc b/lemon/lp_soplex.cc --- a/lemon/lp_soplex.cc +++ b/lemon/lp_soplex.cc @@ -16,8 +16,8 @@ * */ -#include -#include +#include +#include #include @@ -26,54 +26,53 @@ ///\brief Implementation of the LEMON-SOPLEX lp solver interface. namespace lemon { - LpSoplex::LpSoplex() : LpSolverBase() { - rows.setIdHandler(relocateIdHandler); - cols.setIdHandler(relocateIdHandler); + LpSoplex::LpSoplex() { soplex = new soplex::SoPlex; - solved = false; } LpSoplex::~LpSoplex() { delete soplex; } - LpSoplex::LpSoplex(const LpSoplex& lp) : LpSolverBase() { + LpSoplex::LpSoplex(const LpSoplex& lp) { 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; + _col_names = lp._col_names; + _col_names_ref = lp._col_names_ref; - primal_value = lp.primal_value; - dual_value = lp.dual_value; + _row_names = lp._row_names; + _row_names_ref = lp._row_names_ref; } - LpSolverBase* LpSoplex::_newLp() { + void LpSoplex::_clear_temporals() { + _primal_values.clear(); + _dual_values.clear(); + } + + LpSoplex* LpSoplex::_newSolver() const { LpSoplex* newlp = new LpSoplex(); return newlp; } - LpSolverBase* LpSoplex::_copyLp() { + LpSoplex* LpSoplex::_cloneSolver() const { LpSoplex* newlp = new LpSoplex(*this); return newlp; } + const char* LpSoplex::_solverName() const { return "LpSoplex"; } + 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; + _col_names.push_back(std::string()); return soplex->nCols() - 1; } @@ -84,8 +83,7 @@ r.setRhs(soplex::infinity); soplex->addRow(r); - dual_value.push_back(0.0); - solved = false; + _row_names.push_back(std::string()); return soplex->nRows() - 1; } @@ -93,56 +91,84 @@ 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; + _col_names_ref.erase(_col_names[i]); + _col_names[i] = _col_names.back(); + _col_names_ref[_col_names.back()] = i; + _col_names.pop_back(); } void LpSoplex::_eraseRow(int i) { soplex->removeRow(i); - dual_value[i] = dual_value.back(); - dual_value.pop_back(); - solved = false; + _row_names_ref.erase(_row_names[i]); + _row_names[i] = _row_names.back(); + _row_names_ref[_row_names.back()] = i; + _row_names.pop_back(); + } + + void LpSoplex::_eraseColId(int i) { + cols.eraseIndex(i); + cols.relocateIndex(i, cols.maxIndex()); + } + void LpSoplex::_eraseRowId(int i) { + rows.eraseIndex(i); + rows.relocateIndex(i, rows.maxIndex()); } void LpSoplex::_getColName(int c, std::string &name) const { - name = colNames[c]; + name = _col_names[c]; } void LpSoplex::_setColName(int c, const std::string &name) { - invColNames.erase(colNames[c]); - colNames[c] = name; + _col_names_ref.erase(_col_names[c]); + _col_names[c] = name; if (!name.empty()) { - invColNames.insert(std::make_pair(name, c)); + _col_names_ref.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()) { + _col_names_ref.find(name); + if (it != _col_names_ref.end()) { return it->second; } else { return -1; } } + void LpSoplex::_getRowName(int r, std::string &name) const { + name = _row_names[r]; + } - void LpSoplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e) { + void LpSoplex::_setRowName(int r, const std::string &name) { + _row_names_ref.erase(_row_names[r]); + _row_names[r] = name; + if (!name.empty()) { + _row_names_ref.insert(std::make_pair(name, r)); + } + } + + int LpSoplex::_rowByName(const std::string& name) const { + std::map::const_iterator it = + _row_names_ref.find(name); + if (it != _row_names_ref.end()) { + return it->second; + } else { + return -1; + } + } + + + void LpSoplex::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) { for (int j = 0; j < soplex->nCols(); ++j) { soplex->changeElement(i, j, 0.0); } - for(ConstRowIterator it = b; it != e; ++it) { + for(ExprIterator it = b; it != e; ++it) { soplex->changeElement(i, it->first, it->second); } - solved = false; } - void LpSoplex::_getRowCoeffs(int i, RowIterator b) const { + void LpSoplex::_getRowCoeffs(int i, InsertIterator 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)); @@ -150,17 +176,16 @@ } } - void LpSoplex::_setColCoeffs(int j, ConstColIterator b, ConstColIterator e) { + void LpSoplex::_setColCoeffs(int j, ExprIterator b, ExprIterator e) { for (int i = 0; i < soplex->nRows(); ++i) { soplex->changeElement(i, j, 0.0); } - for(ConstColIterator it = b; it != e; ++it) { + for(ExprIterator it = b; it != e; ++it) { soplex->changeElement(it->first, j, it->second); } - solved = false; } - void LpSoplex::_getColCoeffs(int i, ColIterator b) const { + void LpSoplex::_getColCoeffs(int i, InsertIterator 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)); @@ -170,7 +195,6 @@ 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 { @@ -178,8 +202,8 @@ } void LpSoplex::_setColLowerBound(int i, Value value) { + LEMON_ASSERT(value != INF, "Invalid bound"); soplex->changeLower(i, value != -INF ? value : -soplex::infinity); - solved = false; } LpSoplex::Value LpSoplex::_getColLowerBound(int i) const { @@ -188,8 +212,8 @@ } void LpSoplex::_setColUpperBound(int i, Value value) { + LEMON_ASSERT(value != -INF, "Invalid bound"); soplex->changeUpper(i, value != INF ? value : soplex::infinity); - solved = false; } LpSoplex::Value LpSoplex::_getColUpperBound(int i) const { @@ -197,48 +221,63 @@ 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::_setRowLowerBound(int i, Value lb) { + LEMON_ASSERT(lb != INF, "Invalid bound"); + soplex->changeRange(i, lb != -INF ? lb : -soplex::infinity, soplex->rhs(i)); } - 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; + + LpSoplex::Value LpSoplex::_getRowLowerBound(int i) const { + double res = soplex->lhs(i); + return res == -soplex::infinity ? -INF : res; + } + + void LpSoplex::_setRowUpperBound(int i, Value ub) { + LEMON_ASSERT(ub != -INF, "Invalid bound"); + soplex->changeRange(i, soplex->lhs(i), ub != INF ? ub : soplex::infinity); + } + + LpSoplex::Value LpSoplex::_getRowUpperBound(int i) const { + double res = soplex->rhs(i); + return res == soplex::infinity ? INF : res; + } + + void LpSoplex::_setObjCoeffs(ExprIterator b, ExprIterator e) { + for (int j = 0; j < soplex->nCols(); ++j) { + soplex->changeObj(j, 0.0); + } + for (ExprIterator it = b; it != e; ++it) { + soplex->changeObj(it->first, it->second); + } + } + + void LpSoplex::_getObjCoeffs(InsertIterator b) const { + for (int j = 0; j < soplex->nCols(); ++j) { + Value coef = soplex->obj(j); + if (coef != 0.0) { + *b = std::make_pair(j, coef); + ++b; + } + } } 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() { - LpSoplex::SolveExitStatus LpSoplex::_solve() { + _clear_temporals(); + 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; @@ -246,28 +285,87 @@ } LpSoplex::Value LpSoplex::_getPrimal(int i) const { - return primal_value[i]; + if (_primal_values.empty()) { + _primal_values.resize(soplex->nCols()); + soplex::Vector pv(_primal_values.size(), &_primal_values.front()); + soplex->getPrimal(pv); + } + return _primal_values[i]; } LpSoplex::Value LpSoplex::_getDual(int i) const { - return dual_value[i]; + if (_dual_values.empty()) { + _dual_values.resize(soplex->nRows()); + soplex::Vector dv(_dual_values.size(), &_dual_values.front()); + soplex->getDual(dv); + } + return _dual_values[i]; } LpSoplex::Value LpSoplex::_getPrimalValue() const { return soplex->objValue(); } - bool LpSoplex::_isBasicCol(int i) const { - return soplex->getBasisColStatus(i) == soplex::SPxSolver::BASIC; + LpSoplex::VarStatus LpSoplex::_getColStatus(int i) const { + switch (soplex->getBasisColStatus(i)) { + case soplex::SPxSolver::BASIC: + return BASIC; + case soplex::SPxSolver::ON_UPPER: + return UPPER; + case soplex::SPxSolver::ON_LOWER: + return LOWER; + case soplex::SPxSolver::FIXED: + return FIXED; + case soplex::SPxSolver::ZERO: + return FREE; + default: + LEMON_ASSERT(false, "Wrong column status"); + return VarStatus(); + } } - LpSoplex::SolutionStatus LpSoplex::_getPrimalStatus() const { - if (!solved) return UNDEFINED; + LpSoplex::VarStatus LpSoplex::_getRowStatus(int i) const { + switch (soplex->getBasisRowStatus(i)) { + case soplex::SPxSolver::BASIC: + return BASIC; + case soplex::SPxSolver::ON_UPPER: + return UPPER; + case soplex::SPxSolver::ON_LOWER: + return LOWER; + case soplex::SPxSolver::FIXED: + return FIXED; + case soplex::SPxSolver::ZERO: + return FREE; + default: + LEMON_ASSERT(false, "Wrong row status"); + return VarStatus(); + } + } + + LpSoplex::Value LpSoplex::_getPrimalRay(int i) const { + if (_primal_ray.empty()) { + _primal_ray.resize(soplex->nCols()); + soplex::Vector pv(_primal_ray.size(), &_primal_ray.front()); + soplex->getDualfarkas(pv); + } + return _primal_ray[i]; + } + + LpSoplex::Value LpSoplex::_getDualRay(int i) const { + if (_dual_ray.empty()) { + _dual_ray.resize(soplex->nRows()); + soplex::Vector dv(_dual_ray.size(), &_dual_ray.front()); + soplex->getDualfarkas(dv); + } + return _dual_ray[i]; + } + + LpSoplex::ProblemType LpSoplex::_getPrimalType() const { switch (soplex->status()) { case soplex::SPxSolver::OPTIMAL: return OPTIMAL; case soplex::SPxSolver::UNBOUNDED: - return INFINITE; + return UNBOUNDED; case soplex::SPxSolver::INFEASIBLE: return INFEASIBLE; default: @@ -275,42 +373,51 @@ } } - LpSoplex::SolutionStatus LpSoplex::_getDualStatus() const { - if (!solved) return UNDEFINED; + LpSoplex::ProblemType LpSoplex::_getDualType() const { switch (soplex->status()) { case soplex::SPxSolver::OPTIMAL: return OPTIMAL; case soplex::SPxSolver::UNBOUNDED: + return UNBOUNDED; + case soplex::SPxSolver::INFEASIBLE: 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::_setSense(Sense sense) { + switch (sense) { + case MIN: + soplex->changeSense(soplex::SPxSolver::MINIMIZE); + break; + case MAX: + soplex->changeSense(soplex::SPxSolver::MAXIMIZE); } } - 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; + LpSoplex::Sense LpSoplex::_getSense() const { + switch (soplex->spxSense()) { + case soplex::SPxSolver::MAXIMIZE: + return MAX; + case soplex::SPxSolver::MINIMIZE: + return MIN; + default: + LEMON_ASSERT(false, "Wrong sense."); + return LpSoplex::Sense(); + } } + void LpSoplex::_clear() { + soplex->clear(); + _col_names.clear(); + _col_names_ref.clear(); + _row_names.clear(); + _row_names_ref.clear(); + cols.clear(); + rows.clear(); + _clear_temporals(); + } } //namespace lemon diff --git a/lemon/lp_soplex.h b/lemon/lp_soplex.h --- a/lemon/lp_soplex.h +++ b/lemon/lp_soplex.h @@ -43,26 +43,30 @@ /// 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; + class LpSoplex : public LpSolver { + private: soplex::SoPlex* soplex; - bool solved; - std::vector colNames; - std::map invColNames; + std::vector _col_names; + std::map _col_names_ref; - std::vector primal_value; - std::vector dual_value; + std::vector _row_names; + std::map _row_names_ref; + private: + + // these values cannot be retrieved element by element + mutable std::vector _primal_values; + mutable std::vector _dual_values; + + mutable std::vector _primal_ray; + mutable std::vector _dual_ray; + + void _clear_temporals(); public: - typedef LpSolverBase Parent; - - /// \e LpSoplex(); /// \e @@ -72,48 +76,75 @@ protected: - virtual LpSolverBase* _newLp(); - virtual LpSolverBase* _copyLp(); + virtual LpSoplex* _newSolver() const; + virtual LpSoplex* _cloneSolver() const; + + virtual const char* _solverName() const; 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 void _eraseColId(int i); + virtual void _eraseRowId(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 _getRowName(int row, std::string& name) const; + virtual void _setRowName(int row, const std::string& name); + virtual int _rowByName(const std::string& name) const; + + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e); + virtual void _getRowCoeffs(int i, InsertIterator b) const; + + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e); + virtual void _getColCoeffs(int i, InsertIterator 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 _setRowLowerBound(int i, Value value); + virtual Value _getRowLowerBound(int i) const; + virtual void _setRowUpperBound(int i, Value value); + virtual Value _getRowUpperBound(int i) const; + + virtual void _setObjCoeffs(ExprIterator b, ExprIterator e); + virtual void _getObjCoeffs(InsertIterator b) const; + virtual void _setObjCoeff(int i, Value obj_coef); virtual Value _getObjCoeff(int i) const; - virtual void _clearObj(); + + virtual void _setSense(Sense sense); + virtual Sense _getSense() const; 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 Value _getPrimalRay(int i) const; + virtual Value _getDualRay(int i) const; + virtual VarStatus _getColStatus(int i) const; + virtual VarStatus _getRowStatus(int i) const; - virtual void _setMax(); - virtual void _setMin(); - virtual bool _isMax() const; + virtual ProblemType _getPrimalType() const; + virtual ProblemType _getDualType() const; + + virtual void _clear(); }; + } //END OF NAMESPACE LEMON #endif //LEMON_LP_SOPLEX_H diff --git a/lemon/mip_cplex.cc b/lemon/mip_cplex.cc deleted file mode 100644 --- a/lemon/mip_cplex.cc +++ /dev/null @@ -1,141 +0,0 @@ -/* -*- 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 deleted file mode 100644 --- a/lemon/mip_cplex.h +++ /dev/null @@ -1,61 +0,0 @@ -/* -*- 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 deleted file mode 100644 --- a/lemon/mip_glpk.cc +++ /dev/null @@ -1,154 +0,0 @@ -/* -*- 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 deleted file mode 100644 --- a/lemon/mip_glpk.h +++ /dev/null @@ -1,59 +0,0 @@ -/* -*- 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/m4/lx_check_clp.m4 b/m4/lx_check_clp.m4 new file mode 100644 --- /dev/null +++ b/m4/lx_check_clp.m4 @@ -0,0 +1,73 @@ +AC_DEFUN([LX_CHECK_CLP], +[ + AC_ARG_WITH([clp], +AS_HELP_STRING([--with-clp@<:@=PREFIX@:>@], [search for CLP under PREFIX or under the default search paths if PREFIX is not given @<:@default@:>@]) +AS_HELP_STRING([--without-clp], [disable checking for CLP]), + [], [with_clp=yes]) + + AC_ARG_WITH([clp-includedir], +AS_HELP_STRING([--with-clp-includedir=DIR], [search for CLP headers in DIR]), + [], [with_clp_includedir=no]) + + AC_ARG_WITH([clp-libdir], +AS_HELP_STRING([--with-clp-libdir=DIR], [search for CLP libraries in DIR]), + [], [with_clp_libdir=no]) + + lx_clp_found=no + if test x"$with_clp" != x"no"; then + AC_MSG_CHECKING([for CLP]) + + if test x"$with_clp_includedir" != x"no"; then + CLP_CXXFLAGS="-I$with_clp_includedir" + elif test x"$with_clp" != x"yes"; then + CLP_CXXFLAGS="-I$with_clp/include" + fi + + if test x"$with_clp_libdir" != x"no"; then + CLP_LDFLAGS="-L$with_clp_libdir" + elif test x"$with_clp" != x"yes"; then + CLP_LDFLAGS="-L$with_clp/lib" + fi + CLP_LIBS="-lClp -lCoinUtils -lm" + + lx_save_cxxflags="$CXXFLAGS" + lx_save_ldflags="$LDFLAGS" + lx_save_libs="$LIBS" + CXXFLAGS="$CLP_CXXFLAGS" + LDFLAGS="$CLP_LDFLAGS" + LIBS="$CLP_LIBS" + + lx_clp_test_prog=' + #include + + int main(int argc, char** argv) + { + ClpModel clp; + return 0; + }' + + AC_LANG_PUSH(C++) + AC_LINK_IFELSE([$lx_clp_test_prog], [lx_clp_found=yes], [lx_clp_found=no]) + AC_LANG_POP(C++) + + CXXFLAGS="$lx_save_cxxflags" + LDFLAGS="$lx_save_ldflags" + LIBS="$lx_save_libs" + + if test x"$lx_clp_found" = x"yes"; then + AC_DEFINE([HAVE_CLP], [1], [Define to 1 if you have CLP.]) + lx_lp_found=yes + AC_DEFINE([HAVE_LP], [1], [Define to 1 if you have any LP solver.]) + AC_MSG_RESULT([yes]) + else + CLP_CXXFLAGS="" + CLP_LDFLAGS="" + CLP_LIBS="" + AC_MSG_RESULT([no]) + fi + fi + CLP_LIBS="$CLP_LDFLAGS $CLP_LIBS" + AC_SUBST(CLP_CXXFLAGS) + AC_SUBST(CLP_LIBS) + AM_CONDITIONAL([HAVE_CLP], [test x"$lx_clp_found" = x"yes"]) +]) diff --git a/m4/lx_check_glpk.m4 b/m4/lx_check_glpk.m4 --- a/m4/lx_check_glpk.m4 +++ b/m4/lx_check_glpk.m4 @@ -42,6 +42,11 @@ #include } + #if (GLP_MAJOR_VERSION < 4) \ + || (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION < 33) + #error Supported GLPK versions: 4.33 or above + #endif + int main(int argc, char** argv) { LPX *lp; diff --git a/test/lp_test.cc b/test/lp_test.cc --- a/test/lp_test.cc +++ b/test/lp_test.cc @@ -37,14 +37,16 @@ #include #endif +#ifdef HAVE_CLP +#include +#endif + using namespace lemon; -void lpTest(LpSolverBase & lp) +void lpTest(LpSolver& lp) { - - - typedef LpSolverBase LP; + typedef LpSolver LP; std::vector x(10); // for(int i=0;i<10;i++) x.push_back(lp.addCol()); @@ -52,7 +54,6 @@ lp.colLowerBound(x,1); lp.colUpperBound(x,1); lp.colBounds(x,1,2); -#ifndef GYORSITAS std::vector y(10); lp.addColSet(y); @@ -86,11 +87,11 @@ p5=lp.addCol(); e[p1]=2; - e.constComp()=12; + *e=12; e[p1]+=2; - e.constComp()+=12; + *e+=12; e[p1]-=2; - e.constComp()-=12; + *e-=12; e=2; e=2.2; @@ -170,11 +171,11 @@ e[x[3]]=2; e[x[3]]=4; e[x[3]]=1; - e.constComp()=12; + *e=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(-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); @@ -183,21 +184,6 @@ 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<(e)[p2]==0, buf.str()); } @@ -247,36 +233,33 @@ ); } -#endif } -void solveAndCheck(LpSolverBase& lp, LpSolverBase::SolutionStatus stat, +void solveAndCheck(LpSolver& lp, LpSolver::ProblemType stat, double exp_opt) { using std::string; lp.solve(); - //int decimal,sign; + std::ostringstream buf; - buf << "Primalstatus should be: " << int(stat); + buf << "PrimalType should be: " << int(stat) << int(lp.primalType()); - // itoa(stat,buf1, 10); - check(lp.primalStatus()==stat, buf.str()); + check(lp.primalType()==stat, buf.str()); - if (stat == LpSolverBase::OPTIMAL) { + if (stat == LpSolver::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) + check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str()); } } -void aTest(LpSolverBase & lp) +void aTest(LpSolver & lp) { - typedef LpSolverBase LP; + typedef LpSolver LP; //The following example is very simple - typedef LpSolverBase::Row Row; - typedef LpSolverBase::Col Col; + typedef LpSolver::Row Row; + typedef LpSolver::Col Col; Col x1 = lp.addCol(); @@ -284,7 +267,7 @@ //Constraints - Row upright=lp.addRow(x1+x2 <=1); + Row upright=lp.addRow(x1+2*x2 <=1); lp.addRow(x1+x2 >=-1); lp.addRow(x1-x2 <=1); lp.addRow(x1-x2 >=-1); @@ -294,129 +277,126 @@ //Objective function lp.obj(x1+x2); - lp.max(); + lp.sense(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.sense() == lp.MAX,"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->sense() == clp->MAX,"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."); + check(clp->colLowerBound(x1)==0, + "The lower bound for variable x1 should be 0."); + check(clp->colUpperBound(x1)==LpSolver::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."); + check(lp.rowLowerBound(upright)==-LpSolver::INF, + "The lower bound for the first row should be -infty."); + check(lp.rowUpperBound(upright)==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."); + 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."); + 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); + solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt); //Minimization - lp.min(); + lp.sense(lp.MIN); expected_opt=0; - solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt); + solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt); //Vertex (-1,0) instead of (0,0) - lp.colLowerBound(x1, -LpSolverBase::INF); + lp.colLowerBound(x1, -LpSolver::INF); expected_opt=-1; - solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt); + solveAndCheck(lp, LpSolver::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); + lp.erase(upright); + lp.sense(lp.MAX); + expected_opt=LpSolver::INF; + solveAndCheck(lp, LpSolver::UNBOUNDED, 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 = "< +#include #endif #ifdef HAVE_GLPK -#include +#include #endif using namespace lemon; -void solveAndCheck(MipSolverBase& lp, MipSolverBase::SolutionStatus stat, +void solveAndCheck(MipSolver& mip, MipSolver::ProblemType stat, double exp_opt) { using std::string; - lp.solve(); + mip.solve(); //int decimal,sign; std::ostringstream buf; - buf << "Primalstatus should be: " << int(stat) - <<" and it is "<