1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/athos/lp/expression.h Tue Mar 22 11:45:47 2005 +0000
1.3 @@ -0,0 +1,197 @@
1.4 +// -*- c++ -*-
1.5 +#ifndef LEMON_EXPRESSION_H
1.6 +#define LEMON_EXPRESSION_H
1.7 +
1.8 +#include <iostream>
1.9 +#include <map>
1.10 +#include <limits>
1.11 +
1.12 +namespace lemon {
1.13 +
1.14 + /*! \brief Linear expression
1.15 +
1.16 + \c Expr<_Col,_Value> implements a class of linear expressions with the
1.17 + operations of addition and multiplication with scalar.
1.18 +
1.19 + \author Marton Makai
1.20 + */
1.21 + template <typename _Col, typename _Value>
1.22 + class Expr;
1.23 +
1.24 + template <typename _Col, typename _Value>
1.25 + class Expr {
1.26 +// protected:
1.27 + public:
1.28 + typedef
1.29 + typename std::map<_Col, _Value> Data;
1.30 + Data data;
1.31 + public:
1.32 + void simplify() {
1.33 + for (typename Data::iterator i=data.begin();
1.34 + i!=data.end(); ++i) {
1.35 + if ((*i).second==0) data.erase(i);
1.36 + }
1.37 + }
1.38 + Expr() { }
1.39 + Expr(_Col _col) {
1.40 + data.insert(std::make_pair(_col, 1));
1.41 + }
1.42 + Expr& operator*=(_Value _value) {
1.43 + for (typename Data::iterator i=data.begin();
1.44 + i!=data.end(); ++i) {
1.45 + (*i).second *= _value;
1.46 + }
1.47 + simplify();
1.48 + return *this;
1.49 + }
1.50 + Expr& operator+=(const Expr<_Col, _Value>& expr) {
1.51 + for (typename Data::const_iterator j=expr.data.begin();
1.52 + j!=expr.data.end(); ++j) {
1.53 + typename Data::iterator i=data.find((*j).first);
1.54 + if (i==data.end()) {
1.55 + data.insert(std::make_pair((*j).first, (*j).second));
1.56 + } else {
1.57 + (*i).second+=(*j).second;
1.58 + }
1.59 + }
1.60 + simplify();
1.61 + return *this;
1.62 + }
1.63 + Expr& operator-=(const Expr<_Col, _Value>& expr) {
1.64 + for (typename Data::const_iterator j=expr.data.begin();
1.65 + j!=expr.data.end(); ++j) {
1.66 + typename Data::iterator i=data.find((*j).first);
1.67 + if (i==data.end()) {
1.68 + data.insert(std::make_pair((*j).first, -(*j).second));
1.69 + } else {
1.70 + (*i).second+=-(*j).second;
1.71 + }
1.72 + }
1.73 + simplify();
1.74 + return *this;
1.75 + }
1.76 + template <typename _C, typename _V>
1.77 + friend std::ostream& operator<<(std::ostream& os,
1.78 + const Expr<_C, _V>& expr);
1.79 + };
1.80 +
1.81 + template <typename _Col, typename _Value>
1.82 + Expr<_Col, _Value> operator*(_Value _value, _Col _col) {
1.83 + Expr<_Col, _Value> tmp(_col);
1.84 + tmp*=_value;
1.85 + tmp.simplify();
1.86 + return tmp;
1.87 + }
1.88 +
1.89 + template <typename _Col, typename _Value>
1.90 + Expr<_Col, _Value> operator*(_Value _value,
1.91 + const Expr<_Col, _Value>& expr) {
1.92 + Expr<_Col, _Value> tmp(expr);
1.93 + tmp*=_value;
1.94 + tmp.simplify();
1.95 + return tmp;
1.96 + }
1.97 +
1.98 + template <typename _Col, typename _Value>
1.99 + Expr<_Col, _Value> operator+(const Expr<_Col, _Value>& expr1,
1.100 + const Expr<_Col, _Value>& expr2) {
1.101 + Expr<_Col, _Value> tmp(expr1);
1.102 + tmp+=expr2;
1.103 + tmp.simplify();
1.104 + return tmp;
1.105 + }
1.106 +
1.107 + template <typename _Col, typename _Value>
1.108 + Expr<_Col, _Value> operator-(const Expr<_Col, _Value>& expr1,
1.109 + const Expr<_Col, _Value>& expr2) {
1.110 + Expr<_Col, _Value> tmp(expr1);
1.111 + tmp-=expr2;
1.112 + tmp.simplify();
1.113 + return tmp;
1.114 + }
1.115 +
1.116 + template <typename _Col, typename _Value>
1.117 + std::ostream& operator<<(std::ostream& os,
1.118 + const Expr<_Col, _Value>& expr) {
1.119 + for (typename Expr<_Col, _Value>::Data::const_iterator i=
1.120 + expr.data.begin();
1.121 + i!=expr.data.end(); ++i) {
1.122 + os << (*i).second << "*" << (*i).first << " ";
1.123 + }
1.124 + return os;
1.125 + }
1.126 +
1.127 + template <typename _Col, typename _Value>
1.128 + class LConstr {
1.129 + // protected:
1.130 + public:
1.131 + Expr<_Col, _Value> expr;
1.132 + _Value lo;
1.133 + public:
1.134 + LConstr(const Expr<_Col, _Value>& _expr, _Value _lo) :
1.135 + expr(_expr), lo(_lo) { }
1.136 + };
1.137 +
1.138 + template <typename _Col, typename _Value>
1.139 + LConstr<_Col, _Value>
1.140 + operator<=(_Value lo, const Expr<_Col, _Value>& expr) {
1.141 + return LConstr<_Col, _Value>(expr, lo);
1.142 + }
1.143 +
1.144 + template <typename _Col, typename _Value>
1.145 + class UConstr {
1.146 + // protected:
1.147 + public:
1.148 + Expr<_Col, _Value> expr;
1.149 + _Value up;
1.150 + public:
1.151 + UConstr(const Expr<_Col, _Value>& _expr, _Value _up) :
1.152 + expr(_expr), up(_up) { }
1.153 + };
1.154 +
1.155 + template <typename _Col, typename _Value>
1.156 + UConstr<_Col, _Value>
1.157 + operator<=(const Expr<_Col, _Value>& expr, _Value up) {
1.158 + return UConstr<_Col, _Value>(expr, up);
1.159 + }
1.160 +
1.161 + template <typename _Col, typename _Value>
1.162 + class Constr {
1.163 + // protected:
1.164 + public:
1.165 + Expr<_Col, _Value> expr;
1.166 + _Value lo, up;
1.167 + public:
1.168 + Constr(const Expr<_Col, _Value>& _expr, _Value _lo, _Value _up) :
1.169 + expr(_expr), lo(_lo), up(_up) { }
1.170 + Constr(const LConstr<_Col, _Value>& _lconstr) :
1.171 + expr(_lconstr.expr),
1.172 + lo(_lconstr.lo),
1.173 + up(std::numeric_limits<_Value>::infinity()) { }
1.174 + Constr(const UConstr<_Col, _Value>& _uconstr) :
1.175 + expr(_uconstr.expr),
1.176 + lo(-std::numeric_limits<_Value>::infinity()),
1.177 + up(_uconstr.up) { }
1.178 + };
1.179 +
1.180 + template <typename _Col, typename _Value>
1.181 + Constr<_Col, _Value>
1.182 + operator<=(const LConstr<_Col, _Value>& lconstr, _Value up) {
1.183 + return Constr<_Col, _Value>(lconstr.expr, lconstr.lo, up);
1.184 + }
1.185 +
1.186 + template <typename _Col, typename _Value>
1.187 + Constr<_Col, _Value>
1.188 + operator<=(_Value lo, const UConstr<_Col, _Value>& uconstr) {
1.189 + return Constr<_Col, _Value>(uconstr.expr, lo, uconstr.up);
1.190 + }
1.191 +
1.192 + template <typename _Col, typename _Value>
1.193 + Constr<_Col, _Value>
1.194 + operator==(const Expr<_Col, _Value>& expr, _Value value) {
1.195 + return Constr<_Col, _Value>(expr, value, value);
1.196 + }
1.197 +
1.198 +} //namespace lemon
1.199 +
1.200 +#endif //LEMON_EXPRESSION_H
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/work/athos/lp/expression_test.cc Tue Mar 22 11:45:47 2005 +0000
2.3 @@ -0,0 +1,23 @@
2.4 +#include <expression.h>
2.5 +#include <iostream>
2.6 +#include <string>
2.7 +
2.8 +using std::cout;
2.9 +using std::endl;
2.10 +using std::string;
2.11 +using namespace lemon;
2.12 +
2.13 +int main() {
2.14 + Expr<string, double> b;
2.15 + cout << b << endl;
2.16 + Expr<string, double> c("f");
2.17 + cout << c << endl;
2.18 + Expr<string, double> d=8.0*string("g");
2.19 + cout << d << endl;
2.20 + c*=5;
2.21 + cout << c << endl;
2.22 + Expr<string, double> e=c;
2.23 + e+=8.9*9.0*string("l");
2.24 + cout << e << endl;
2.25 + cout << c+d << endl;
2.26 +}
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/src/work/athos/lp/lp_solver_base.h Tue Mar 22 11:45:47 2005 +0000
3.3 @@ -0,0 +1,1116 @@
3.4 +// -*- c++ -*-
3.5 +#ifndef LEMON_LP_SOLVER_BASE_H
3.6 +#define LEMON_LP_SOLVER_BASE_H
3.7 +
3.8 +///\ingroup misc
3.9 +///\file
3.10 +
3.11 +// #include <stdio.h>
3.12 +#include <stdlib.h>
3.13 +#include <iostream>
3.14 +#include <map>
3.15 +#include <limits>
3.16 +// #include <stdio>
3.17 +//#include <stdlib>
3.18 +extern "C" {
3.19 +#include "glpk.h"
3.20 +}
3.21 +
3.22 +#include <iostream>
3.23 +#include <vector>
3.24 +#include <string>
3.25 +#include <list>
3.26 +#include <memory>
3.27 +#include <utility>
3.28 +
3.29 +#include <lemon/invalid.h>
3.30 +#include <expression.h>
3.31 +//#include <stp.h>
3.32 +//#include <lemon/max_flow.h>
3.33 +//#include <augmenting_flow.h>
3.34 +//#include <iter_map.h>
3.35 +
3.36 +using std::cout;
3.37 +using std::cin;
3.38 +using std::endl;
3.39 +
3.40 +namespace lemon {
3.41 +
3.42 + /// \addtogroup misc
3.43 + /// @{
3.44 +
3.45 + /// \brief A partitioned vector with iterable classes.
3.46 + ///
3.47 + /// This class implements a container in which the data is stored in an
3.48 + /// stl vector, the range is partitioned into sets and each set is
3.49 + /// doubly linked in a list.
3.50 + /// That is, each class is iterable by lemon iterators, and any member of
3.51 + /// the vector can bo moved to an other class.
3.52 + template <typename T>
3.53 + class IterablePartition {
3.54 + protected:
3.55 + struct Node {
3.56 + T data;
3.57 + int prev; //invalid az -1
3.58 + int next;
3.59 + };
3.60 + std::vector<Node> nodes;
3.61 + struct Tip {
3.62 + int first;
3.63 + int last;
3.64 + };
3.65 + std::vector<Tip> tips;
3.66 + public:
3.67 + /// The classes are indexed by integers from \c 0 to \c classNum()-1.
3.68 + int classNum() const { return tips.size(); }
3.69 + /// This lemon style iterator iterates through a class.
3.70 + class Class;
3.71 + /// Constructor. The number of classes is to be given which is fixed
3.72 + /// over the life of the container.
3.73 + /// The partition classes are indexed from 0 to class_num-1.
3.74 + IterablePartition(int class_num) {
3.75 + for (int i=0; i<class_num; ++i) {
3.76 + Tip t;
3.77 + t.first=t.last=-1;
3.78 + tips.push_back(t);
3.79 + }
3.80 + }
3.81 + protected:
3.82 + void befuz(Class it, int class_id) {
3.83 + if (tips[class_id].first==-1) {
3.84 + if (tips[class_id].last==-1) {
3.85 + nodes[it.i].prev=nodes[it.i].next=-1;
3.86 + tips[class_id].first=tips[class_id].last=it.i;
3.87 + }
3.88 + } else {
3.89 + nodes[it.i].prev=tips[class_id].last;
3.90 + nodes[it.i].next=-1;
3.91 + nodes[tips[class_id].last].next=it.i;
3.92 + tips[class_id].last=it.i;
3.93 + }
3.94 + }
3.95 + void kifuz(Class it, int class_id) {
3.96 + if (tips[class_id].first==it.i) {
3.97 + if (tips[class_id].last==it.i) {
3.98 + tips[class_id].first=tips[class_id].last=-1;
3.99 + } else {
3.100 + tips[class_id].first=nodes[it.i].next;
3.101 + nodes[nodes[it.i].next].prev=-1;
3.102 + }
3.103 + } else {
3.104 + if (tips[class_id].last==it.i) {
3.105 + tips[class_id].last=nodes[it.i].prev;
3.106 + nodes[nodes[it.i].prev].next=-1;
3.107 + } else {
3.108 + nodes[nodes[it.i].next].prev=nodes[it.i].prev;
3.109 + nodes[nodes[it.i].prev].next=nodes[it.i].next;
3.110 + }
3.111 + }
3.112 + }
3.113 + public:
3.114 + /// A new element with data \c t is pushed into the vector and into class
3.115 + /// \c class_id.
3.116 + Class push_back(const T& t, int class_id) {
3.117 + Node n;
3.118 + n.data=t;
3.119 + nodes.push_back(n);
3.120 + int i=nodes.size()-1;
3.121 + befuz(i, class_id);
3.122 + return i;
3.123 + }
3.124 + /// A member is moved to an other class.
3.125 + void set(Class it, int old_class_id, int new_class_id) {
3.126 + kifuz(it.i, old_class_id);
3.127 + befuz(it.i, new_class_id);
3.128 + }
3.129 + /// Returns the data pointed by \c it.
3.130 + T& operator[](Class it) { return nodes[it.i].data; }
3.131 + /// Returns the data pointed by \c it.
3.132 + const T& operator[](Class it) const { return nodes[it.i].data; }
3.133 + ///.
3.134 + class Class {
3.135 + friend class IterablePartition;
3.136 + protected:
3.137 + int i;
3.138 + public:
3.139 + /// Default constructor.
3.140 + Class() { }
3.141 + /// This constructor constructs an iterator which points
3.142 + /// to the member of th container indexed by the integer _i.
3.143 + Class(const int& _i) : i(_i) { }
3.144 + /// Invalid constructor.
3.145 + Class(const Invalid&) : i(-1) { }
3.146 + friend bool operator<(const Class& x, const Class& y);
3.147 + friend std::ostream& operator<<(std::ostream& os,
3.148 + const Class& it);
3.149 + bool operator==(const Class& node) const {return i == node.i;}
3.150 + bool operator!=(const Class& node) const {return i != node.i;}
3.151 + };
3.152 + friend bool operator<(const Class& x, const Class& y) {
3.153 + return (x.i < y.i);
3.154 + }
3.155 + friend std::ostream& operator<<(std::ostream& os,
3.156 + const Class& it) {
3.157 + os << it.i;
3.158 + return os;
3.159 + }
3.160 + /// First member of class \c class_id.
3.161 + Class& first(Class& it, int class_id) const {
3.162 + it.i=tips[class_id].first;
3.163 + return it;
3.164 + }
3.165 + /// Next member.
3.166 + Class& next(Class& it) const {
3.167 + it.i=nodes[it.i].next;
3.168 + return it;
3.169 + }
3.170 + /// True iff the iterator is valid.
3.171 + bool valid(const Class& it) const { return it.i!=-1; }
3.172 +
3.173 + class ClassIt : public Class {
3.174 + const IterablePartition* iterable_partition;
3.175 + public:
3.176 + ClassIt() { }
3.177 + ClassIt(Invalid i) : Class(i) { }
3.178 + ClassIt(const IterablePartition& _iterable_partition,
3.179 + const int& i) : iterable_partition(&_iterable_partition) {
3.180 + _iterable_partition.first(*this, i);
3.181 + }
3.182 + ClassIt(const IterablePartition& _iterable_partition,
3.183 + const Class& _class) :
3.184 + Class(_class), iterable_partition(&_iterable_partition) { }
3.185 + ClassIt& operator++() {
3.186 + iterable_partition->next(*this);
3.187 + return *this;
3.188 + }
3.189 + };
3.190 +
3.191 + };
3.192 +
3.193 +
3.194 + /*! \e
3.195 + \todo kellenene uj iterable structure bele, mert ez nem az igazi
3.196 + \todo A[x,y]-t cserel. Jobboldal, baloldal csere.
3.197 + \todo LEKERDEZESEK!!!
3.198 + \todo DOKSI!!!! Doxygen group!!!
3.199 + The aim of this class is to give a general surface to different
3.200 + solvers, i.e. it makes possible to write algorithms using LP's,
3.201 + in which the solver can be changed to an other one easily.
3.202 + \nosubgrouping
3.203 + */
3.204 + template <typename _Value>
3.205 + class LPSolverBase {
3.206 +
3.207 + /*! @name Uncategorized functions and types (public members)
3.208 + */
3.209 + //@{
3.210 + public:
3.211 +
3.212 + //UNCATEGORIZED
3.213 +
3.214 + /// \e
3.215 + typedef IterablePartition<int> Rows;
3.216 + /// \e
3.217 + typedef IterablePartition<int> Cols;
3.218 + /// \e
3.219 + typedef _Value Value;
3.220 + /// \e
3.221 + typedef Rows::Class Row;
3.222 + /// \e
3.223 + typedef Cols::Class Col;
3.224 + public:
3.225 + /// \e
3.226 + IterablePartition<int> row_iter_map;
3.227 + /// \e
3.228 + IterablePartition<int> col_iter_map;
3.229 + /// \e
3.230 + std::vector<Row> int_row_map;
3.231 + /// \e
3.232 + std::vector<Col> int_col_map;
3.233 + /// \e
3.234 + const int VALID_CLASS;
3.235 + /// \e
3.236 + const int INVALID_CLASS;
3.237 + /// \e
3.238 + static const _Value INF;
3.239 + public:
3.240 + /// \e
3.241 + LPSolverBase() : row_iter_map(2),
3.242 + col_iter_map(2),
3.243 + VALID_CLASS(0), INVALID_CLASS(1) { }
3.244 + /// \e
3.245 + virtual ~LPSolverBase() { }
3.246 + //@}
3.247 +
3.248 + /*! @name Medium level interface (public members)
3.249 + These functions appear in the low level and also in the high level
3.250 + interfaces thus these each of these functions have to be implemented
3.251 + only once in the different interfaces.
3.252 + This means that these functions have to be reimplemented for all of the
3.253 + different lp solvers. These are basic functions, and have the same
3.254 + parameter lists in the low and high level interfaces.
3.255 + */
3.256 + //@{
3.257 + public:
3.258 +
3.259 + //UNCATEGORIZED FUNCTIONS
3.260 +
3.261 + /// \e
3.262 + virtual void setMinimize() = 0;
3.263 + /// \e
3.264 + virtual void setMaximize() = 0;
3.265 +
3.266 + //SOLVER FUNCTIONS
3.267 +
3.268 + /// \e
3.269 + virtual void solveSimplex() = 0;
3.270 + /// \e
3.271 + virtual void solvePrimalSimplex() = 0;
3.272 + /// \e
3.273 + virtual void solveDualSimplex() = 0;
3.274 +
3.275 + //SOLUTION RETRIEVING
3.276 +
3.277 + /// \e
3.278 + virtual _Value getObjVal() = 0;
3.279 +
3.280 + //OTHER FUNCTIONS
3.281 +
3.282 + /// \e
3.283 + virtual int rowNum() const = 0;
3.284 + /// \e
3.285 + virtual int colNum() const = 0;
3.286 + /// \e
3.287 + virtual int warmUp() = 0;
3.288 + /// \e
3.289 + virtual void printWarmUpStatus(int i) = 0;
3.290 + /// \e
3.291 + virtual int getPrimalStatus() = 0;
3.292 + /// \e
3.293 + virtual void printPrimalStatus(int i) = 0;
3.294 + /// \e
3.295 + virtual int getDualStatus() = 0;
3.296 + /// \e
3.297 + virtual void printDualStatus(int i) = 0;
3.298 + /// Returns the status of the slack variable assigned to row \c row.
3.299 + virtual int getRowStat(const Row& row) = 0;
3.300 + /// \e
3.301 + virtual void printRowStatus(int i) = 0;
3.302 + /// Returns the status of the variable assigned to column \c col.
3.303 + virtual int getColStat(const Col& col) = 0;
3.304 + /// \e
3.305 + virtual void printColStatus(int i) = 0;
3.306 +
3.307 + //@}
3.308 +
3.309 + /*! @name Low level interface (protected members)
3.310 + Problem manipulating functions in the low level interface
3.311 + */
3.312 + //@{
3.313 + protected:
3.314 +
3.315 + //MATRIX MANIPULATING FUNCTIONS
3.316 +
3.317 + /// \e
3.318 + virtual int _addCol() = 0;
3.319 + /// \e
3.320 + virtual int _addRow() = 0;
3.321 + /// \e
3.322 + virtual void _eraseCol(int i) = 0;
3.323 + /// \e
3.324 + virtual void _eraseRow(int i) = 0;
3.325 + /// \e
3.326 + virtual void _setRowCoeffs(int i,
3.327 + const std::vector<std::pair<int, _Value> >& coeffs) = 0;
3.328 + /// \e
3.329 + /// This routine modifies \c coeffs only by the \c push_back method.
3.330 + virtual void _getRowCoeffs(int i,
3.331 + std::vector<std::pair<int, _Value> >& coeffs) = 0;
3.332 + /// \e
3.333 + virtual void _setColCoeffs(int i,
3.334 + const std::vector<std::pair<int, _Value> >& coeffs) = 0;
3.335 + /// \e
3.336 + /// This routine modifies \c coeffs only by the \c push_back method.
3.337 + virtual void _getColCoeffs(int i,
3.338 + std::vector<std::pair<int, _Value> >& coeffs) = 0;
3.339 + /// \e
3.340 + virtual void _setCoeff(int col, int row, _Value value) = 0;
3.341 + /// \e
3.342 + virtual _Value _getCoeff(int col, int row) = 0;
3.343 + // public:
3.344 + // /// \e
3.345 + // enum Bound { FREE, LOWER, UPPER, DOUBLE, FIXED };
3.346 + protected:
3.347 + /// \e
3.348 + /// The lower bound of a variable (column) have to be given by an
3.349 + /// extended number of type _Value, i.e. a finite number of type
3.350 + /// _Value or -INF.
3.351 + virtual void _setColLowerBound(int i, _Value value) = 0;
3.352 + /// \e
3.353 + /// The lower bound of a variable (column) is an
3.354 + /// extended number of type _Value, i.e. a finite number of type
3.355 + /// _Value or -INF.
3.356 + virtual _Value _getColLowerBound(int i) = 0;
3.357 + /// \e
3.358 + /// The upper bound of a variable (column) have to be given by an
3.359 + /// extended number of type _Value, i.e. a finite number of type
3.360 + /// _Value or INF.
3.361 + virtual void _setColUpperBound(int i, _Value value) = 0;
3.362 + /// \e
3.363 + /// The upper bound of a variable (column) is an
3.364 + /// extended number of type _Value, i.e. a finite number of type
3.365 + /// _Value or INF.
3.366 + virtual _Value _getColUpperBound(int i) = 0;
3.367 + /// \e
3.368 + /// The lower bound of a linear expression (row) have to be given by an
3.369 + /// extended number of type _Value, i.e. a finite number of type
3.370 + /// _Value or -INF.
3.371 + virtual void _setRowLowerBound(int i, _Value value) = 0;
3.372 + /// \e
3.373 + /// The lower bound of a linear expression (row) is an
3.374 + /// extended number of type _Value, i.e. a finite number of type
3.375 + /// _Value or -INF.
3.376 + virtual _Value _getRowLowerBound(int i) = 0;
3.377 + /// \e
3.378 + /// The upper bound of a linear expression (row) have to be given by an
3.379 + /// extended number of type _Value, i.e. a finite number of type
3.380 + /// _Value or INF.
3.381 + virtual void _setRowUpperBound(int i, _Value value) = 0;
3.382 + /// \e
3.383 + /// The upper bound of a linear expression (row) is an
3.384 + /// extended number of type _Value, i.e. a finite number of type
3.385 + /// _Value or INF.
3.386 + virtual _Value _getRowUpperBound(int i) = 0;
3.387 + /// \e
3.388 + virtual void _setObjCoeff(int i, _Value obj_coef) = 0;
3.389 + /// \e
3.390 + virtual _Value _getObjCoeff(int i) = 0;
3.391 +
3.392 + //SOLUTION RETRIEVING
3.393 +
3.394 + /// \e
3.395 + virtual _Value _getPrimal(int i) = 0;
3.396 + //@}
3.397 +
3.398 + /*! @name High level interface (public members)
3.399 + Problem manipulating functions in the high level interface
3.400 + */
3.401 + //@{
3.402 + public:
3.403 +
3.404 + //MATRIX MANIPULATING FUNCTIONS
3.405 +
3.406 + /// \e
3.407 + Col addCol() {
3.408 + int i=_addCol();
3.409 + Col col;
3.410 + col_iter_map.first(col, INVALID_CLASS);
3.411 + if (col_iter_map.valid(col)) { //van hasznalhato hely
3.412 + col_iter_map.set(col, INVALID_CLASS, VALID_CLASS);
3.413 + col_iter_map[col]=i;
3.414 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
3.415 + col=col_iter_map.push_back(i, VALID_CLASS);
3.416 + }
3.417 + int_col_map.push_back(col);
3.418 + return col;
3.419 + }
3.420 + /// \e
3.421 + Row addRow() {
3.422 + int i=_addRow();
3.423 + Row row;
3.424 + row_iter_map.first(row, INVALID_CLASS);
3.425 + if (row_iter_map.valid(row)) { //van hasznalhato hely
3.426 + row_iter_map.set(row, INVALID_CLASS, VALID_CLASS);
3.427 + row_iter_map[row]=i;
3.428 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
3.429 + row=row_iter_map.push_back(i, VALID_CLASS);
3.430 + }
3.431 + int_row_map.push_back(row);
3.432 + return row;
3.433 + }
3.434 + /// \e
3.435 + void eraseCol(const Col& col) {
3.436 + col_iter_map.set(col, VALID_CLASS, INVALID_CLASS);
3.437 + int cols[2];
3.438 + cols[1]=col_iter_map[col];
3.439 + _eraseCol(cols[1]);
3.440 + col_iter_map[col]=0; //glpk specifikus, de kell ez??
3.441 + Col it;
3.442 + for (col_iter_map.first(it, VALID_CLASS);
3.443 + col_iter_map.valid(it); col_iter_map.next(it)) {
3.444 + if (col_iter_map[it]>cols[1]) --col_iter_map[it];
3.445 + }
3.446 + int_col_map.erase(int_col_map.begin()+cols[1]);
3.447 + }
3.448 + /// \e
3.449 + void eraseRow(const Row& row) {
3.450 + row_iter_map.set(row, VALID_CLASS, INVALID_CLASS);
3.451 + int rows[2];
3.452 + rows[1]=row_iter_map[row];
3.453 + _eraseRow(rows[1]);
3.454 + row_iter_map[row]=0; //glpk specifikus, de kell ez??
3.455 + Row it;
3.456 + for (row_iter_map.first(it, VALID_CLASS);
3.457 + row_iter_map.valid(it); row_iter_map.next(it)) {
3.458 + if (row_iter_map[it]>rows[1]) --row_iter_map[it];
3.459 + }
3.460 + int_row_map.erase(int_row_map.begin()+rows[1]);
3.461 + }
3.462 + /// \e
3.463 + void setCoeff(Col col, Row row, _Value value) {
3.464 + _setCoeff(col_iter_map[col], row_iter_map[row], value);
3.465 + }
3.466 + /// \e
3.467 + _Value getCoeff(Col col, Row row) {
3.468 + return _getCoeff(col_iter_map[col], row_iter_map[row], value);
3.469 + }
3.470 + /// \e
3.471 + void setColLowerBound(Col col, _Value lo) {
3.472 + _setColLowerBound(col_iter_map[col], lo);
3.473 + }
3.474 + /// \e
3.475 + _Value getColLowerBound(Col col) {
3.476 + return _getColLowerBound(col_iter_map[col]);
3.477 + }
3.478 + /// \e
3.479 + void setColUpperBound(Col col, _Value up) {
3.480 + _setColUpperBound(col_iter_map[col], up);
3.481 + }
3.482 + /// \e
3.483 + _Value getColUpperBound(Col col) {
3.484 + return _getColUpperBound(col_iter_map[col]);
3.485 + }
3.486 + /// \e
3.487 + void setRowLowerBound(Row row, _Value lo) {
3.488 + _setRowLowerBound(row_iter_map[row], lo);
3.489 + }
3.490 + /// \e
3.491 + _Value getRowLowerBound(Row row) {
3.492 + return _getRowLowerBound(row_iter_map[row]);
3.493 + }
3.494 + /// \e
3.495 + void setRowUpperBound(Row row, _Value up) {
3.496 + _setRowUpperBound(row_iter_map[row], up);
3.497 + }
3.498 + /// \e
3.499 + _Value getRowUpperBound(Row row) {
3.500 + return _getRowUpperBound(row_iter_map[row]);
3.501 + }
3.502 + /// \e
3.503 + void setObjCoeff(const Col& col, _Value obj_coef) {
3.504 + _setObjCoeff(col_iter_map[col], obj_coef);
3.505 + }
3.506 + /// \e
3.507 + _Value getObjCoeff(const Col& col) {
3.508 + return _getObjCoeff(col_iter_map[col]);
3.509 + }
3.510 +
3.511 + //SOLUTION RETRIEVING FUNCTIONS
3.512 +
3.513 + /// \e
3.514 + _Value getPrimal(const Col& col) {
3.515 + return _getPrimal(col_iter_map[col]);
3.516 + }
3.517 +
3.518 + //@}
3.519 +
3.520 + /*! @name User friend interface
3.521 + Problem manipulating functions in the user friend interface
3.522 + */
3.523 + //@{
3.524 +
3.525 + //EXPRESSION TYPES
3.526 +
3.527 + /// \e
3.528 + typedef Expr<Col, _Value> Expression;
3.529 + /// \e
3.530 + typedef Expr<Row, _Value> DualExpression;
3.531 + /// \e
3.532 + typedef Constr<Col, _Value> Constraint;
3.533 +
3.534 + //MATRIX MANIPULATING FUNCTIONS
3.535 +
3.536 + /// \e
3.537 + void setRowCoeffs(Row row, const Expression& expr) {
3.538 + std::vector<std::pair<int, _Value> > row_coeffs;
3.539 + for(typename Expression::Data::const_iterator i=expr.data.begin();
3.540 + i!=expr.data.end(); ++i) {
3.541 + row_coeffs.push_back(std::make_pair
3.542 + (col_iter_map[(*i).first], (*i).second));
3.543 + }
3.544 + _setRowCoeffs(row_iter_map[row], row_coeffs);
3.545 + }
3.546 + /// \e
3.547 + void setRow(Row row, const Constraint& constr) {
3.548 + setRowCoeffs(row, constr.expr);
3.549 + setRowLowerBound(row, constr.lo);
3.550 + setRowUpperBound(row, constr.up);
3.551 + }
3.552 + /// \e
3.553 + Row addRow(const Constraint& constr) {
3.554 + Row row=addRow();
3.555 + setRowCoeffs(row, constr.expr);
3.556 + setRowLowerBound(row, constr.lo);
3.557 + setRowUpperBound(row, constr.up);
3.558 + return row;
3.559 + }
3.560 + /// \e
3.561 + /// This routine modifies \c expr by only adding to it.
3.562 + void getRowCoeffs(Row row, Expression& expr) {
3.563 + std::vector<std::pair<int, _Value> > row_coeffs;
3.564 + _getRowCoeffs(row_iter_map[row], row_coeffs);
3.565 + for(typename std::vector<std::pair<int, _Value> >::const_iterator
3.566 + i=row_coeffs.begin(); i!=row_coeffs.end(); ++i) {
3.567 + expr+= (*i).second*int_col_map[(*i).first];
3.568 + }
3.569 + }
3.570 + /// \e
3.571 + void setColCoeffs(Col col, const DualExpression& expr) {
3.572 + std::vector<std::pair<int, _Value> > col_coeffs;
3.573 + for(typename DualExpression::Data::const_iterator i=expr.data.begin();
3.574 + i!=expr.data.end(); ++i) {
3.575 + col_coeffs.push_back(std::make_pair
3.576 + (row_iter_map[(*i).first], (*i).second));
3.577 + }
3.578 + _setColCoeffs(col_iter_map[col], col_coeffs);
3.579 + }
3.580 + /// \e
3.581 + /// This routine modifies \c expr by only adding to it.
3.582 + void getColCoeffs(Col col, DualExpression& expr) {
3.583 + std::vector<std::pair<int, _Value> > col_coeffs;
3.584 + _getColCoeffs(col_iter_map[col], col_coeffs);
3.585 + for(typename std::vector<std::pair<int, _Value> >::const_iterator
3.586 + i=col_coeffs.begin(); i!=col_coeffs.end(); ++i) {
3.587 + expr+= (*i).second*int_row_map[(*i).first];
3.588 + }
3.589 + }
3.590 + /// \e
3.591 + void setObjCoeffs(const Expression& expr) {
3.592 + // writing zero everywhere
3.593 + for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
3.594 + setObjCoeff(it, 0.0);
3.595 + // writing the data needed
3.596 + for(typename Expression::Data::const_iterator i=expr.data.begin();
3.597 + i!=expr.data.end(); ++i) {
3.598 + setObjCoeff((*i).first, (*i).second);
3.599 + }
3.600 + }
3.601 + /// \e
3.602 + /// This routine modifies \c expr by only adding to it.
3.603 + void getObjCoeffs(Expression& expr) {
3.604 + for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
3.605 + expr+=getObjCoeff(it)*it;
3.606 + }
3.607 + //@}
3.608 +
3.609 +
3.610 + /*! @name MIP functions and types (public members)
3.611 + */
3.612 + //@{
3.613 + public:
3.614 + /// \e
3.615 + virtual void solveBandB() = 0;
3.616 + /// \e
3.617 + virtual void setLP() = 0;
3.618 + /// \e
3.619 + virtual void setMIP() = 0;
3.620 + protected:
3.621 + /// \e
3.622 + virtual void _setColCont(int i) = 0;
3.623 + /// \e
3.624 + virtual void _setColInt(int i) = 0;
3.625 + /// \e
3.626 + virtual _Value _getMIPPrimal(int i) = 0;
3.627 + public:
3.628 + /// \e
3.629 + void setColCont(Col col) {
3.630 + _setColCont(col_iter_map[col]);
3.631 + }
3.632 + /// \e
3.633 + void setColInt(Col col) {
3.634 + _setColInt(col_iter_map[col]);
3.635 + }
3.636 + /// \e
3.637 + _Value getMIPPrimal(Col col) {
3.638 + return _getMIPPrimal(col_iter_map[col]);
3.639 + }
3.640 + //@}
3.641 + };
3.642 +
3.643 + template <typename _Value>
3.644 + const _Value LPSolverBase<_Value>::INF=std::numeric_limits<_Value>::infinity();
3.645 +
3.646 +
3.647 + /// \brief Wrapper for GLPK solver
3.648 + ///
3.649 + /// This class implements a lemon wrapper for GLPK.
3.650 + class LPGLPK : public LPSolverBase<double> {
3.651 + public:
3.652 + typedef LPSolverBase<double> Parent;
3.653 +
3.654 + public:
3.655 + /// \e
3.656 + LPX* lp;
3.657 +
3.658 + public:
3.659 + /// \e
3.660 + LPGLPK() : Parent(),
3.661 + lp(lpx_create_prob()) {
3.662 + int_row_map.push_back(Row());
3.663 + int_col_map.push_back(Col());
3.664 + lpx_set_int_parm(lp, LPX_K_DUAL, 1);
3.665 + }
3.666 + /// \e
3.667 + ~LPGLPK() {
3.668 + lpx_delete_prob(lp);
3.669 + }
3.670 +
3.671 + //MATRIX INDEPEDENT MANIPULATING FUNCTIONS
3.672 +
3.673 + /// \e
3.674 + void setMinimize() {
3.675 + lpx_set_obj_dir(lp, LPX_MIN);
3.676 + }
3.677 + /// \e
3.678 + void setMaximize() {
3.679 + lpx_set_obj_dir(lp, LPX_MAX);
3.680 + }
3.681 +
3.682 + //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS
3.683 +
3.684 + protected:
3.685 + /// \e
3.686 + int _addCol() {
3.687 + int i=lpx_add_cols(lp, 1);
3.688 + _setColLowerBound(i, -INF);
3.689 + _setColUpperBound(i, INF);
3.690 + return i;
3.691 + }
3.692 + /// \e
3.693 + int _addRow() {
3.694 + int i=lpx_add_rows(lp, 1);
3.695 + return i;
3.696 + }
3.697 + /// \e
3.698 + virtual void _setRowCoeffs(int i,
3.699 + const std::vector<std::pair<int, double> >& coeffs) {
3.700 + int mem_length=1+colNum();
3.701 + int* indices = new int[mem_length];
3.702 + double* doubles = new double[mem_length];
3.703 + int length=0;
3.704 + for (std::vector<std::pair<int, double> >::
3.705 + const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
3.706 + ++length;
3.707 + indices[length]=it->first;
3.708 + doubles[length]=it->second;
3.709 + }
3.710 + lpx_set_mat_row(lp, i, length, indices, doubles);
3.711 + delete [] indices;
3.712 + delete [] doubles;
3.713 + }
3.714 + /// \e
3.715 + virtual void _getRowCoeffs(int i,
3.716 + std::vector<std::pair<int, double> >& coeffs) {
3.717 + int mem_length=1+colNum();
3.718 + int* indices = new int[mem_length];
3.719 + double* doubles = new double[mem_length];
3.720 + int length=lpx_get_mat_row(lp, i, indices, doubles);
3.721 + for (int i=1; i<=length; ++i) {
3.722 + coeffs.push_back(std::make_pair(indices[i], doubles[i]));
3.723 + }
3.724 + delete [] indices;
3.725 + delete [] doubles;
3.726 + }
3.727 + /// \e
3.728 + virtual void _setColCoeffs(int i,
3.729 + const std::vector<std::pair<int, double> >& coeffs) {
3.730 + int mem_length=1+rowNum();
3.731 + int* indices = new int[mem_length];
3.732 + double* doubles = new double[mem_length];
3.733 + int length=0;
3.734 + for (std::vector<std::pair<int, double> >::
3.735 + const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
3.736 + ++length;
3.737 + indices[length]=it->first;
3.738 + doubles[length]=it->second;
3.739 + }
3.740 + lpx_set_mat_col(lp, i, length, indices, doubles);
3.741 + delete [] indices;
3.742 + delete [] doubles;
3.743 + }
3.744 + /// \e
3.745 + virtual void _getColCoeffs(int i,
3.746 + std::vector<std::pair<int, double> >& coeffs) {
3.747 + int mem_length=1+rowNum();
3.748 + int* indices = new int[mem_length];
3.749 + double* doubles = new double[mem_length];
3.750 + int length=lpx_get_mat_col(lp, i, indices, doubles);
3.751 + for (int i=1; i<=length; ++i) {
3.752 + coeffs.push_back(std::make_pair(indices[i], doubles[i]));
3.753 + }
3.754 + delete [] indices;
3.755 + delete [] doubles;
3.756 + }
3.757 + /// \e
3.758 + virtual void _eraseCol(int i) {
3.759 + int cols[2];
3.760 + cols[1]=i;
3.761 + lpx_del_cols(lp, 1, cols);
3.762 + }
3.763 + virtual void _eraseRow(int i) {
3.764 + int rows[2];
3.765 + rows[1]=i;
3.766 + lpx_del_rows(lp, 1, rows);
3.767 + }
3.768 + void _setCoeff(int col, int row, double value) {
3.769 + /// FIXME not yet implemented
3.770 + }
3.771 + double _getCoeff(int col, int row) {
3.772 + /// FIXME not yet implemented
3.773 + return 0.0;
3.774 + }
3.775 + virtual void _setColLowerBound(int i, double lo) {
3.776 + if (lo==INF) {
3.777 + //FIXME error
3.778 + }
3.779 + int b=lpx_get_col_type(lp, i);
3.780 + double up=lpx_get_col_ub(lp, i);
3.781 + if (lo==-INF) {
3.782 + switch (b) {
3.783 + case LPX_FR:
3.784 + case LPX_LO:
3.785 + lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
3.786 + break;
3.787 + case LPX_UP:
3.788 + break;
3.789 + case LPX_DB:
3.790 + case LPX_FX:
3.791 + lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
3.792 + break;
3.793 + default: ;
3.794 + //FIXME error
3.795 + }
3.796 + } else {
3.797 + switch (b) {
3.798 + case LPX_FR:
3.799 + case LPX_LO:
3.800 + lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
3.801 + break;
3.802 + case LPX_UP:
3.803 + case LPX_DB:
3.804 + case LPX_FX:
3.805 + if (lo==up)
3.806 + lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
3.807 + else
3.808 + lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
3.809 + break;
3.810 + default: ;
3.811 + //FIXME error
3.812 + }
3.813 + }
3.814 + }
3.815 + virtual double _getColLowerBound(int i) {
3.816 + int b=lpx_get_col_type(lp, i);
3.817 + switch (b) {
3.818 + case LPX_FR:
3.819 + return -INF;
3.820 + case LPX_LO:
3.821 + return lpx_get_col_lb(lp, i);
3.822 + case LPX_UP:
3.823 + return -INF;
3.824 + case LPX_DB:
3.825 + case LPX_FX:
3.826 + return lpx_get_col_lb(lp, i);
3.827 + default: ;
3.828 + //FIXME error
3.829 + return 0.0;
3.830 + }
3.831 + }
3.832 + virtual void _setColUpperBound(int i, double up) {
3.833 + if (up==-INF) {
3.834 + //FIXME error
3.835 + }
3.836 + int b=lpx_get_col_type(lp, i);
3.837 + double lo=lpx_get_col_lb(lp, i);
3.838 + if (up==INF) {
3.839 + switch (b) {
3.840 + case LPX_FR:
3.841 + case LPX_LO:
3.842 + break;
3.843 + case LPX_UP:
3.844 + lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
3.845 + break;
3.846 + case LPX_DB:
3.847 + case LPX_FX:
3.848 + lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
3.849 + break;
3.850 + default: ;
3.851 + //FIXME error
3.852 + }
3.853 + } else {
3.854 + switch (b) {
3.855 + case LPX_FR:
3.856 + lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
3.857 + case LPX_LO:
3.858 + if (lo==up)
3.859 + lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
3.860 + else
3.861 + lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
3.862 + break;
3.863 + case LPX_UP:
3.864 + lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
3.865 + break;
3.866 + case LPX_DB:
3.867 + case LPX_FX:
3.868 + if (lo==up)
3.869 + lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
3.870 + else
3.871 + lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
3.872 + break;
3.873 + default: ;
3.874 + //FIXME error
3.875 + }
3.876 + }
3.877 + }
3.878 + virtual double _getColUpperBound(int i) {
3.879 + int b=lpx_get_col_type(lp, i);
3.880 + switch (b) {
3.881 + case LPX_FR:
3.882 + case LPX_LO:
3.883 + return INF;
3.884 + case LPX_UP:
3.885 + case LPX_DB:
3.886 + case LPX_FX:
3.887 + return lpx_get_col_ub(lp, i);
3.888 + default: ;
3.889 + //FIXME error
3.890 + return 0.0;
3.891 + }
3.892 + }
3.893 + virtual void _setRowLowerBound(int i, double lo) {
3.894 + if (lo==INF) {
3.895 + //FIXME error
3.896 + }
3.897 + int b=lpx_get_row_type(lp, i);
3.898 + double up=lpx_get_row_ub(lp, i);
3.899 + if (lo==-INF) {
3.900 + switch (b) {
3.901 + case LPX_FR:
3.902 + case LPX_LO:
3.903 + lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
3.904 + break;
3.905 + case LPX_UP:
3.906 + break;
3.907 + case LPX_DB:
3.908 + case LPX_FX:
3.909 + lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
3.910 + break;
3.911 + default: ;
3.912 + //FIXME error
3.913 + }
3.914 + } else {
3.915 + switch (b) {
3.916 + case LPX_FR:
3.917 + case LPX_LO:
3.918 + lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
3.919 + break;
3.920 + case LPX_UP:
3.921 + case LPX_DB:
3.922 + case LPX_FX:
3.923 + if (lo==up)
3.924 + lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
3.925 + else
3.926 + lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
3.927 + break;
3.928 + default: ;
3.929 + //FIXME error
3.930 + }
3.931 + }
3.932 + }
3.933 + virtual double _getRowLowerBound(int i) {
3.934 + int b=lpx_get_row_type(lp, i);
3.935 + switch (b) {
3.936 + case LPX_FR:
3.937 + return -INF;
3.938 + case LPX_LO:
3.939 + return lpx_get_row_lb(lp, i);
3.940 + case LPX_UP:
3.941 + return -INF;
3.942 + case LPX_DB:
3.943 + case LPX_FX:
3.944 + return lpx_get_row_lb(lp, i);
3.945 + default: ;
3.946 + //FIXME error
3.947 + return 0.0;
3.948 + }
3.949 + }
3.950 + virtual void _setRowUpperBound(int i, double up) {
3.951 + if (up==-INF) {
3.952 + //FIXME error
3.953 + }
3.954 + int b=lpx_get_row_type(lp, i);
3.955 + double lo=lpx_get_row_lb(lp, i);
3.956 + if (up==INF) {
3.957 + switch (b) {
3.958 + case LPX_FR:
3.959 + case LPX_LO:
3.960 + break;
3.961 + case LPX_UP:
3.962 + lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
3.963 + break;
3.964 + case LPX_DB:
3.965 + case LPX_FX:
3.966 + lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
3.967 + break;
3.968 + default: ;
3.969 + //FIXME error
3.970 + }
3.971 + } else {
3.972 + switch (b) {
3.973 + case LPX_FR:
3.974 + lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
3.975 + case LPX_LO:
3.976 + if (lo==up)
3.977 + lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
3.978 + else
3.979 + lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
3.980 + break;
3.981 + case LPX_UP:
3.982 + lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
3.983 + break;
3.984 + case LPX_DB:
3.985 + case LPX_FX:
3.986 + if (lo==up)
3.987 + lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
3.988 + else
3.989 + lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
3.990 + break;
3.991 + default: ;
3.992 + //FIXME error
3.993 + }
3.994 + }
3.995 + }
3.996 + virtual double _getRowUpperBound(int i) {
3.997 + int b=lpx_get_row_type(lp, i);
3.998 + switch (b) {
3.999 + case LPX_FR:
3.1000 + case LPX_LO:
3.1001 + return INF;
3.1002 + case LPX_UP:
3.1003 + case LPX_DB:
3.1004 + case LPX_FX:
3.1005 + return lpx_get_row_ub(lp, i);
3.1006 + default: ;
3.1007 + //FIXME error
3.1008 + return 0.0;
3.1009 + }
3.1010 + }
3.1011 + /// \e
3.1012 + virtual double _getObjCoeff(int i) {
3.1013 + return lpx_get_obj_coef(lp, i);
3.1014 + }
3.1015 + /// \e
3.1016 + virtual void _setObjCoeff(int i, double obj_coef) {
3.1017 + lpx_set_obj_coef(lp, i, obj_coef);
3.1018 + }
3.1019 + public:
3.1020 + /// \e
3.1021 + void solveSimplex() { lpx_simplex(lp); }
3.1022 + /// \e
3.1023 + void solvePrimalSimplex() { lpx_simplex(lp); }
3.1024 + /// \e
3.1025 + void solveDualSimplex() { lpx_simplex(lp); }
3.1026 + protected:
3.1027 + virtual double _getPrimal(int i) {
3.1028 + return lpx_get_col_prim(lp, i);
3.1029 + }
3.1030 + public:
3.1031 + /// \e
3.1032 + double getObjVal() { return lpx_get_obj_val(lp); }
3.1033 + /// \e
3.1034 + int rowNum() const { return lpx_get_num_rows(lp); }
3.1035 + /// \e
3.1036 + int colNum() const { return lpx_get_num_cols(lp); }
3.1037 + /// \e
3.1038 + int warmUp() { return lpx_warm_up(lp); }
3.1039 + /// \e
3.1040 + void printWarmUpStatus(int i) {
3.1041 + switch (i) {
3.1042 + case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
3.1043 + case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;
3.1044 + case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
3.1045 + case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
3.1046 + }
3.1047 + }
3.1048 + /// \e
3.1049 + int getPrimalStatus() { return lpx_get_prim_stat(lp); }
3.1050 + /// \e
3.1051 + void printPrimalStatus(int i) {
3.1052 + switch (i) {
3.1053 + case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
3.1054 + case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;
3.1055 + case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
3.1056 + case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
3.1057 + }
3.1058 + }
3.1059 + /// \e
3.1060 + int getDualStatus() { return lpx_get_dual_stat(lp); }
3.1061 + /// \e
3.1062 + void printDualStatus(int i) {
3.1063 + switch (i) {
3.1064 + case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
3.1065 + case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;
3.1066 + case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
3.1067 + case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
3.1068 + }
3.1069 + }
3.1070 + /// Returns the status of the slack variable assigned to row \c row.
3.1071 + int getRowStat(const Row& row) {
3.1072 + return lpx_get_row_stat(lp, row_iter_map[row]);
3.1073 + }
3.1074 + /// \e
3.1075 + void printRowStatus(int i) {
3.1076 + switch (i) {
3.1077 + case LPX_BS: cout << "LPX_BS" << endl; break;
3.1078 + case LPX_NL: cout << "LPX_NL" << endl; break;
3.1079 + case LPX_NU: cout << "LPX_NU" << endl; break;
3.1080 + case LPX_NF: cout << "LPX_NF" << endl; break;
3.1081 + case LPX_NS: cout << "LPX_NS" << endl; break;
3.1082 + }
3.1083 + }
3.1084 + /// Returns the status of the variable assigned to column \c col.
3.1085 + int getColStat(const Col& col) {
3.1086 + return lpx_get_col_stat(lp, col_iter_map[col]);
3.1087 + }
3.1088 + /// \e
3.1089 + void printColStatus(int i) {
3.1090 + switch (i) {
3.1091 + case LPX_BS: cout << "LPX_BS" << endl; break;
3.1092 + case LPX_NL: cout << "LPX_NL" << endl; break;
3.1093 + case LPX_NU: cout << "LPX_NU" << endl; break;
3.1094 + case LPX_NF: cout << "LPX_NF" << endl; break;
3.1095 + case LPX_NS: cout << "LPX_NS" << endl; break;
3.1096 + }
3.1097 + }
3.1098 +
3.1099 + // MIP
3.1100 + /// \e
3.1101 + void solveBandB() { lpx_integer(lp); }
3.1102 + /// \e
3.1103 + void setLP() { lpx_set_class(lp, LPX_LP); }
3.1104 + /// \e
3.1105 + void setMIP() { lpx_set_class(lp, LPX_MIP); }
3.1106 + protected:
3.1107 + /// \e
3.1108 + void _setColCont(int i) { lpx_set_col_kind(lp, i, LPX_CV); }
3.1109 + /// \e
3.1110 + void _setColInt(int i) { lpx_set_col_kind(lp, i, LPX_IV); }
3.1111 + /// \e
3.1112 + double _getMIPPrimal(int i) { return lpx_mip_col_val(lp, i); }
3.1113 + };
3.1114 +
3.1115 + /// @}
3.1116 +
3.1117 +} //namespace lemon
3.1118 +
3.1119 +#endif //LEMON_LP_SOLVER_BASE_H
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/src/work/athos/lp/lp_solver_wrapper.h Tue Mar 22 11:45:47 2005 +0000
4.3 @@ -0,0 +1,431 @@
4.4 +// -*- c++ -*-
4.5 +#ifndef LEMON_LP_SOLVER_WRAPPER_H
4.6 +#define LEMON_LP_SOLVER_WRAPPER_H
4.7 +
4.8 +///\ingroup misc
4.9 +///\file
4.10 +///\brief Dijkstra algorithm.
4.11 +
4.12 +// #include <stdio.h>
4.13 +#include <stdlib.h>
4.14 +// #include <stdio>
4.15 +//#include <stdlib>
4.16 +extern "C" {
4.17 +#include "glpk.h"
4.18 +}
4.19 +
4.20 +#include <iostream>
4.21 +#include <vector>
4.22 +#include <string>
4.23 +#include <list>
4.24 +#include <memory>
4.25 +#include <utility>
4.26 +
4.27 +//#include <sage_graph.h>
4.28 +//#include <lemon/list_graph.h>
4.29 +//#include <lemon/graph_wrapper.h>
4.30 +#include <lemon/invalid.h>
4.31 +//#include <bfs_dfs.h>
4.32 +//#include <stp.h>
4.33 +//#include <lemon/max_flow.h>
4.34 +//#include <augmenting_flow.h>
4.35 +//#include <iter_map.h>
4.36 +
4.37 +using std::cout;
4.38 +using std::cin;
4.39 +using std::endl;
4.40 +
4.41 +namespace lemon {
4.42 +
4.43 +
4.44 + /// \addtogroup misc
4.45 + /// @{
4.46 +
4.47 + /// \brief A partitioned vector with iterable classes.
4.48 + ///
4.49 + /// This class implements a container in which the data is stored in an
4.50 + /// stl vector, the range is partitioned into sets and each set is
4.51 + /// doubly linked in a list.
4.52 + /// That is, each class is iterable by lemon iterators, and any member of
4.53 + /// the vector can bo moved to an other class.
4.54 + template <typename T>
4.55 + class IterablePartition {
4.56 + protected:
4.57 + struct Node {
4.58 + T data;
4.59 + int prev; //invalid az -1
4.60 + int next;
4.61 + };
4.62 + std::vector<Node> nodes;
4.63 + struct Tip {
4.64 + int first;
4.65 + int last;
4.66 + };
4.67 + std::vector<Tip> tips;
4.68 + public:
4.69 + /// The classes are indexed by integers from \c 0 to \c classNum()-1.
4.70 + int classNum() const { return tips.size(); }
4.71 + /// This lemon style iterator iterates through a class.
4.72 + class ClassIt;
4.73 + /// Constructor. The number of classes is to be given which is fixed
4.74 + /// over the life of the container.
4.75 + /// The partition classes are indexed from 0 to class_num-1.
4.76 + IterablePartition(int class_num) {
4.77 + for (int i=0; i<class_num; ++i) {
4.78 + Tip t;
4.79 + t.first=t.last=-1;
4.80 + tips.push_back(t);
4.81 + }
4.82 + }
4.83 + protected:
4.84 + void befuz(ClassIt it, int class_id) {
4.85 + if (tips[class_id].first==-1) {
4.86 + if (tips[class_id].last==-1) {
4.87 + nodes[it.i].prev=nodes[it.i].next=-1;
4.88 + tips[class_id].first=tips[class_id].last=it.i;
4.89 + }
4.90 + } else {
4.91 + nodes[it.i].prev=tips[class_id].last;
4.92 + nodes[it.i].next=-1;
4.93 + nodes[tips[class_id].last].next=it.i;
4.94 + tips[class_id].last=it.i;
4.95 + }
4.96 + }
4.97 + void kifuz(ClassIt it, int class_id) {
4.98 + if (tips[class_id].first==it.i) {
4.99 + if (tips[class_id].last==it.i) {
4.100 + tips[class_id].first=tips[class_id].last=-1;
4.101 + } else {
4.102 + tips[class_id].first=nodes[it.i].next;
4.103 + nodes[nodes[it.i].next].prev=-1;
4.104 + }
4.105 + } else {
4.106 + if (tips[class_id].last==it.i) {
4.107 + tips[class_id].last=nodes[it.i].prev;
4.108 + nodes[nodes[it.i].prev].next=-1;
4.109 + } else {
4.110 + nodes[nodes[it.i].next].prev=nodes[it.i].prev;
4.111 + nodes[nodes[it.i].prev].next=nodes[it.i].next;
4.112 + }
4.113 + }
4.114 + }
4.115 + public:
4.116 + /// A new element with data \c t is pushed into the vector and into class
4.117 + /// \c class_id.
4.118 + ClassIt push_back(const T& t, int class_id) {
4.119 + Node n;
4.120 + n.data=t;
4.121 + nodes.push_back(n);
4.122 + int i=nodes.size()-1;
4.123 + befuz(i, class_id);
4.124 + return i;
4.125 + }
4.126 + /// A member is moved to an other class.
4.127 + void set(ClassIt it, int old_class_id, int new_class_id) {
4.128 + kifuz(it.i, old_class_id);
4.129 + befuz(it.i, new_class_id);
4.130 + }
4.131 + /// Returns the data pointed by \c it.
4.132 + T& operator[](ClassIt it) { return nodes[it.i].data; }
4.133 + /// Returns the data pointed by \c it.
4.134 + const T& operator[](ClassIt it) const { return nodes[it.i].data; }
4.135 + ///.
4.136 + class ClassIt {
4.137 + friend class IterablePartition;
4.138 + protected:
4.139 + int i;
4.140 + public:
4.141 + /// Default constructor.
4.142 + ClassIt() { }
4.143 + /// This constructor constructs an iterator which points
4.144 + /// to the member of th container indexed by the integer _i.
4.145 + ClassIt(const int& _i) : i(_i) { }
4.146 + /// Invalid constructor.
4.147 + ClassIt(const Invalid&) : i(-1) { }
4.148 + };
4.149 + /// First member of class \c class_id.
4.150 + ClassIt& first(ClassIt& it, int class_id) const {
4.151 + it.i=tips[class_id].first;
4.152 + return it;
4.153 + }
4.154 + /// Next member.
4.155 + ClassIt& next(ClassIt& it) const {
4.156 + it.i=nodes[it.i].next;
4.157 + return it;
4.158 + }
4.159 + /// True iff the iterator is valid.
4.160 + bool valid(const ClassIt& it) const { return it.i!=-1; }
4.161 + };
4.162 +
4.163 + /// \brief Wrappers for LP solvers
4.164 + ///
4.165 + /// This class implements a lemon wrapper for glpk.
4.166 + /// Later other LP-solvers will be wrapped into lemon.
4.167 + /// The aim of this class is to give a general surface to different
4.168 + /// solvers, i.e. it makes possible to write algorithms using LP's,
4.169 + /// in which the solver can be changed to an other one easily.
4.170 + class LPSolverWrapper {
4.171 + public:
4.172 +
4.173 +// class Row {
4.174 +// protected:
4.175 +// int i;
4.176 +// public:
4.177 +// Row() { }
4.178 +// Row(const Invalid&) : i(0) { }
4.179 +// Row(const int& _i) : i(_i) { }
4.180 +// operator int() const { return i; }
4.181 +// };
4.182 +// class RowIt : public Row {
4.183 +// public:
4.184 +// RowIt(const Row& row) : Row(row) { }
4.185 +// };
4.186 +
4.187 +// class Col {
4.188 +// protected:
4.189 +// int i;
4.190 +// public:
4.191 +// Col() { }
4.192 +// Col(const Invalid&) : i(0) { }
4.193 +// Col(const int& _i) : i(_i) { }
4.194 +// operator int() const { return i; }
4.195 +// };
4.196 +// class ColIt : public Col {
4.197 +// ColIt(const Col& col) : Col(col) { }
4.198 +// };
4.199 +
4.200 + public:
4.201 + ///.
4.202 + LPX* lp;
4.203 + ///.
4.204 + typedef IterablePartition<int>::ClassIt RowIt;
4.205 + ///.
4.206 + IterablePartition<int> row_iter_map;
4.207 + ///.
4.208 + typedef IterablePartition<int>::ClassIt ColIt;
4.209 + ///.
4.210 + IterablePartition<int> col_iter_map;
4.211 + //std::vector<int> row_id_to_lp_row_id;
4.212 + //std::vector<int> col_id_to_lp_col_id;
4.213 + ///.
4.214 + const int VALID_ID;
4.215 + ///.
4.216 + const int INVALID_ID;
4.217 +
4.218 + public:
4.219 + ///.
4.220 + LPSolverWrapper() : lp(lpx_create_prob()),
4.221 + row_iter_map(2),
4.222 + col_iter_map(2),
4.223 + //row_id_to_lp_row_id(), col_id_to_lp_col_id(),
4.224 + VALID_ID(0), INVALID_ID(1) {
4.225 + lpx_set_int_parm(lp, LPX_K_DUAL, 1);
4.226 + }
4.227 + ///.
4.228 + ~LPSolverWrapper() {
4.229 + lpx_delete_prob(lp);
4.230 + }
4.231 + ///.
4.232 + void setMinimize() {
4.233 + lpx_set_obj_dir(lp, LPX_MIN);
4.234 + }
4.235 + ///.
4.236 + void setMaximize() {
4.237 + lpx_set_obj_dir(lp, LPX_MAX);
4.238 + }
4.239 + ///.
4.240 + ColIt addCol() {
4.241 + int i=lpx_add_cols(lp, 1);
4.242 + ColIt col_it;
4.243 + col_iter_map.first(col_it, INVALID_ID);
4.244 + if (col_iter_map.valid(col_it)) { //van hasznalhato hely
4.245 + col_iter_map.set(col_it, INVALID_ID, VALID_ID);
4.246 + col_iter_map[col_it]=i;
4.247 + //col_id_to_lp_col_id[col_iter_map[col_it]]=i;
4.248 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
4.249 + //col_id_to_lp_col_id.push_back(i);
4.250 + //int j=col_id_to_lp_col_id.size()-1;
4.251 + col_it=col_iter_map.push_back(i, VALID_ID);
4.252 + }
4.253 +// edge_index_map.set(e, i);
4.254 +// lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
4.255 +// lpx_set_obj_coef(lp, i, cost[e]);
4.256 + return col_it;
4.257 + }
4.258 + ///.
4.259 + RowIt addRow() {
4.260 + int i=lpx_add_rows(lp, 1);
4.261 + RowIt row_it;
4.262 + row_iter_map.first(row_it, INVALID_ID);
4.263 + if (row_iter_map.valid(row_it)) { //van hasznalhato hely
4.264 + row_iter_map.set(row_it, INVALID_ID, VALID_ID);
4.265 + row_iter_map[row_it]=i;
4.266 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
4.267 + row_it=row_iter_map.push_back(i, VALID_ID);
4.268 + }
4.269 + return row_it;
4.270 + }
4.271 + //pair<RowIt, double>-bol kell megadni egy std range-et
4.272 + ///.
4.273 + template <typename Begin, typename End>
4.274 + void setColCoeffs(const ColIt& col_it,
4.275 + Begin begin, End end) {
4.276 + int mem_length=1+lpx_get_num_rows(lp);
4.277 + int* indices = new int[mem_length];
4.278 + double* doubles = new double[mem_length];
4.279 + int length=0;
4.280 + for ( ; begin!=end; ++begin) {
4.281 + ++length;
4.282 + indices[length]=row_iter_map[begin->first];
4.283 + doubles[length]=begin->second;
4.284 + }
4.285 + lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
4.286 + delete [] indices;
4.287 + delete [] doubles;
4.288 + }
4.289 + //pair<ColIt, double>-bol kell megadni egy std range-et
4.290 + ///.
4.291 + template <typename Begin, typename End>
4.292 + void setRowCoeffs(const RowIt& row_it,
4.293 + Begin begin, End end) {
4.294 + int mem_length=1+lpx_get_num_cols(lp);
4.295 + int* indices = new int[mem_length];
4.296 + double* doubles = new double[mem_length];
4.297 + int length=0;
4.298 + for ( ; begin!=end; ++begin) {
4.299 + ++length;
4.300 + indices[length]=col_iter_map[begin->first];
4.301 + doubles[length]=begin->second;
4.302 + }
4.303 + lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
4.304 + delete [] indices;
4.305 + delete [] doubles;
4.306 + }
4.307 + ///.
4.308 + void eraseCol(const ColIt& col_it) {
4.309 + col_iter_map.set(col_it, VALID_ID, INVALID_ID);
4.310 + int cols[2];
4.311 + cols[1]=col_iter_map[col_it];
4.312 + lpx_del_cols(lp, 1, cols);
4.313 + col_iter_map[col_it]=0; //glpk specifikus
4.314 + ColIt it;
4.315 + for (col_iter_map.first(it, VALID_ID);
4.316 + col_iter_map.valid(it); col_iter_map.next(it)) {
4.317 + if (col_iter_map[it]>cols[1]) --col_iter_map[it];
4.318 + }
4.319 + }
4.320 + ///.
4.321 + void eraseRow(const RowIt& row_it) {
4.322 + row_iter_map.set(row_it, VALID_ID, INVALID_ID);
4.323 + int rows[2];
4.324 + rows[1]=row_iter_map[row_it];
4.325 + lpx_del_rows(lp, 1, rows);
4.326 + row_iter_map[row_it]=0; //glpk specifikus
4.327 + RowIt it;
4.328 + for (row_iter_map.first(it, VALID_ID);
4.329 + row_iter_map.valid(it); row_iter_map.next(it)) {
4.330 + if (row_iter_map[it]>rows[1]) --row_iter_map[it];
4.331 + }
4.332 + }
4.333 + ///.
4.334 + void setColBounds(const ColIt& col_it, int bound_type,
4.335 + double lo, double up) {
4.336 + lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
4.337 + }
4.338 + ///.
4.339 + double getObjCoef(const ColIt& col_it) {
4.340 + return lpx_get_obj_coef(lp, col_iter_map[col_it]);
4.341 + }
4.342 + ///.
4.343 + void setRowBounds(const RowIt& row_it, int bound_type,
4.344 + double lo, double up) {
4.345 + lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
4.346 + }
4.347 + ///.
4.348 + void setObjCoef(const ColIt& col_it, double obj_coef) {
4.349 + lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
4.350 + }
4.351 + ///.
4.352 + void solveSimplex() { lpx_simplex(lp); }
4.353 + ///.
4.354 + void solvePrimalSimplex() { lpx_simplex(lp); }
4.355 + ///.
4.356 + void solveDualSimplex() { lpx_simplex(lp); }
4.357 + ///.
4.358 + double getPrimal(const ColIt& col_it) {
4.359 + return lpx_get_col_prim(lp, col_iter_map[col_it]);
4.360 + }
4.361 + ///.
4.362 + double getObjVal() { return lpx_get_obj_val(lp); }
4.363 + ///.
4.364 + int rowNum() const { return lpx_get_num_rows(lp); }
4.365 + ///.
4.366 + int colNum() const { return lpx_get_num_cols(lp); }
4.367 + ///.
4.368 + int warmUp() { return lpx_warm_up(lp); }
4.369 + ///.
4.370 + void printWarmUpStatus(int i) {
4.371 + switch (i) {
4.372 + case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
4.373 + case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;
4.374 + case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
4.375 + case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
4.376 + }
4.377 + }
4.378 + ///.
4.379 + int getPrimalStatus() { return lpx_get_prim_stat(lp); }
4.380 + ///.
4.381 + void printPrimalStatus(int i) {
4.382 + switch (i) {
4.383 + case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
4.384 + case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;
4.385 + case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
4.386 + case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
4.387 + }
4.388 + }
4.389 + ///.
4.390 + int getDualStatus() { return lpx_get_dual_stat(lp); }
4.391 + ///.
4.392 + void printDualStatus(int i) {
4.393 + switch (i) {
4.394 + case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
4.395 + case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;
4.396 + case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
4.397 + case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
4.398 + }
4.399 + }
4.400 + /// Returns the status of the slack variable assigned to row \c row_it.
4.401 + int getRowStat(const RowIt& row_it) {
4.402 + return lpx_get_row_stat(lp, row_iter_map[row_it]);
4.403 + }
4.404 + ///.
4.405 + void printRowStatus(int i) {
4.406 + switch (i) {
4.407 + case LPX_BS: cout << "LPX_BS" << endl; break;
4.408 + case LPX_NL: cout << "LPX_NL" << endl; break;
4.409 + case LPX_NU: cout << "LPX_NU" << endl; break;
4.410 + case LPX_NF: cout << "LPX_NF" << endl; break;
4.411 + case LPX_NS: cout << "LPX_NS" << endl; break;
4.412 + }
4.413 + }
4.414 + /// Returns the status of the variable assigned to column \c col_it.
4.415 + int getColStat(const ColIt& col_it) {
4.416 + return lpx_get_col_stat(lp, col_iter_map[col_it]);
4.417 + }
4.418 + ///.
4.419 + void printColStatus(int i) {
4.420 + switch (i) {
4.421 + case LPX_BS: cout << "LPX_BS" << endl; break;
4.422 + case LPX_NL: cout << "LPX_NL" << endl; break;
4.423 + case LPX_NU: cout << "LPX_NU" << endl; break;
4.424 + case LPX_NF: cout << "LPX_NF" << endl; break;
4.425 + case LPX_NS: cout << "LPX_NS" << endl; break;
4.426 + }
4.427 + }
4.428 + };
4.429 +
4.430 + /// @}
4.431 +
4.432 +} //namespace lemon
4.433 +
4.434 +#endif //LEMON_LP_SOLVER_WRAPPER_H
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/src/work/athos/lp/magic_square.cc Tue Mar 22 11:45:47 2005 +0000
5.3 @@ -0,0 +1,90 @@
5.4 +// -*- c++ -*-
5.5 +#include <iostream>
5.6 +#include <fstream>
5.7 +
5.8 +#include <lemon/time_measure.h>
5.9 +#include <lp_solver_base.h>
5.10 +
5.11 +using std::cout;
5.12 +using std::endl;
5.13 +using namespace lemon;
5.14 +
5.15 +/*
5.16 + On an 1537Mhz PC, the run times with
5.17 + glpk are the following.
5.18 + for n=3,4, some secondes
5.19 + for n=5, 25 hours
5.20 + */
5.21 +
5.22 +int main(int, char **) {
5.23 + const int n=4;
5.24 + const double row_sum=(1.0+n*n)*n/2;
5.25 + Timer ts;
5.26 + ts.reset();
5.27 + typedef LPGLPK LPSolver;
5.28 + typedef LPSolver::Col Col;
5.29 + LPSolver lp;
5.30 + typedef std::map<std::pair<int, int>, Col> Coords;
5.31 + Coords x;
5.32 + // we create a new variable for each entry
5.33 + // of the magic square
5.34 + for (int i=1; i<=n; ++i) {
5.35 + for (int j=1; j<=n; ++j) {
5.36 + Col col=lp.addCol();
5.37 + x[std::make_pair(i,j)]=col;
5.38 + lp.setColLowerBound(col, 1.0);
5.39 + lp.setColUpperBound(col, double(n*n));
5.40 + }
5.41 + }
5.42 + LPSolver::Expression expr3, expr4;
5.43 + for (int i=1; i<=n; ++i) {
5.44 + LPSolver::Expression expr1, expr2;
5.45 + for (int j=1; j<=n; ++j) {
5.46 + expr1+=x[std::make_pair(i, j)];
5.47 + expr2+=x[std::make_pair(j, i)];
5.48 + }
5.49 + // sum of rows and columns
5.50 + lp.addRow(expr1==row_sum);
5.51 + lp.addRow(expr2==row_sum);
5.52 + expr3+=x[std::make_pair(i, i)];
5.53 + expr4+=x[std::make_pair(i, (n+1)-i)];
5.54 + }
5.55 + // sum of the diagonal entries
5.56 + lp.addRow(expr3==row_sum);
5.57 + lp.addRow(expr4==row_sum);
5.58 + lp.solveSimplex();
5.59 + cout << "elapsed time: " << ts << endl;
5.60 + for (int i=1; i<=n; ++i) {
5.61 + for (int j=1; j<=n; ++j) {
5.62 + cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)])
5.63 + << endl;
5.64 + }
5.65 + }
5.66 + // we make new binary variables for each pair of
5.67 + // entries of the square to achieve that
5.68 + // the values of different entries are different
5.69 + lp.setMIP();
5.70 + for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) {
5.71 + Coords::const_iterator jt=it; ++jt;
5.72 + for(; jt!=x.end(); ++jt) {
5.73 + Col col1=(*it).second;
5.74 + Col col2=(*jt).second;
5.75 + Col col=lp.addCol();
5.76 + lp.setColLowerBound(col, 0.0);
5.77 + lp.setColUpperBound(col, 1.0);
5.78 + lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0);
5.79 + lp.setColInt(col);
5.80 + }
5.81 + }
5.82 + cout << "elapsed time: " << ts << endl;
5.83 + lp.solveSimplex();
5.84 + // let's solve the integer problem
5.85 + lp.solveBandB();
5.86 + cout << "elapsed time: " << ts << endl;
5.87 + for (int i=1; i<=n; ++i) {
5.88 + for (int j=1; j<=n; ++j) {
5.89 + cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)])
5.90 + << endl;
5.91 + }
5.92 + }
5.93 +}
6.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
6.2 +++ b/src/work/athos/lp/makefile Tue Mar 22 11:45:47 2005 +0000
6.3 @@ -0,0 +1,73 @@
6.4 +#INCLUDEDIRS ?= -I.. -I. -I./{marci,jacint,alpar,klao,akos}
6.5 +#GLPKROOT = /usr/local/glpk-4.4
6.6 +INCLUDEDIRS ?= -I/home/marci/boost -I/usr/local/cplex/cplex75/include -I../../{marci,alpar,klao,akos,athos} -I. -I../../.. -I../.. -I..# -I$(GLPKROOT)/include
6.7 +#INCLUDEDIRS ?= -I../.. -I../.. -I../../{marci,jacint,alpar,klao,akos} -I/usr/local/glpk-4.4/include
6.8 +CXXFLAGS = -g -O2 -W -Wall $(INCLUDEDIRS) -ansi -pedantic
6.9 +LDFLAGS = -lglpk#-lcplex -lm -lpthread -lilocplex -L/usr/local/cplex/cplex75/lib/i86_linux2_glibc2.2_gcc3.0/static_mt# -L$(GLPKROOT)/lib
6.10 +
6.11 +BINARIES = magic_square max_flow_expression expression_test max_flow_by_lp# sample sample2 sample11 sample15
6.12 +
6.13 +#include ../makefile
6.14 +
6.15 +# Hat, ez elismerem, hogy nagyon ronda, de mukodik minden altalam
6.16 +# ismert rendszeren :-) (Misi)
6.17 +ifdef GCCVER
6.18 +CXX := g++-$(GCCVER)
6.19 +else
6.20 +CXX := $(shell type -p g++-3.3 || type -p g++-3.2 || type -p g++-3.0 || type -p g++-3 || echo g++)
6.21 +endif
6.22 +
6.23 +ifdef DEBUG
6.24 +CXXFLAGS += -DDEBUG
6.25 +endif
6.26 +
6.27 +CC := $(CXX)
6.28 +
6.29 +
6.30 +all: $(BINARIES)
6.31 +
6.32 +################
6.33 +# Minden binarishoz egy sor, felsorolva, hogy mely object file-okbol
6.34 +# all elo.
6.35 +# Kiveve ha siman file.cc -> file esetrol van szo, amikor is nem kell
6.36 +# irni semmit.
6.37 +
6.38 +#proba: proba.o seged.o
6.39 +
6.40 +################
6.41 +
6.42 +
6.43 +# .depend dep depend:
6.44 +# -$(CXX) $(CXXFLAGS) -M $(BINARIES:=.cc) > .depend
6.45 +
6.46 +#makefile: .depend
6.47 +#sinclude .depend
6.48 +#moert nem megy az eredeti /usr/bin/ld-vel?
6.49 +
6.50 +# %: %.o
6.51 +# $(CXX) -o $@ $< $(LDFLAGS)
6.52 +
6.53 +# %.o: %.cc
6.54 +# $(CXX) $(CXXFLAGS) -c $<
6.55 +
6.56 +%: %.cc
6.57 + $(CXX) $(CXXFLAGS) -o $@ $< $(LDFLAGS)
6.58 +
6.59 +sample11prof: sample11prof.o
6.60 + $(CXX) -pg -o sample11prof sample11prof.o -L$(GLPKROOT)/lib -lglpk
6.61 +sample11prof.o: sample11.cc
6.62 + $(CXX) -pg $(CXXFLAGS) -c -o sample11prof.o sample11.cc
6.63 +
6.64 +# sample.o: sample.cc
6.65 +# $(CXX) $(CXXFLAGS) -c -o sample.o sample.cc
6.66 +
6.67 +# sample2: sample2.o
6.68 +# $(CXX) -o sample2 sample2.o -L/usr/local/glpk-4.4/lib -lglpk
6.69 +# sample2.o: sample2.cc
6.70 +# $(CXX) $(CXXFLAGS) -c -o sample2.o sample2.cc
6.71 +
6.72 +
6.73 +clean:
6.74 + $(RM) *.o $(BINARIES) .depend
6.75 +
6.76 +.PHONY: all clean dep depend
7.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
7.2 +++ b/src/work/athos/lp/max_flow_by_lp.cc Tue Mar 22 11:45:47 2005 +0000
7.3 @@ -0,0 +1,186 @@
7.4 +// -*- c++ -*-
7.5 +#include <iostream>
7.6 +#include <fstream>
7.7 +
7.8 +#include <lemon/smart_graph.h>
7.9 +#include <lemon/list_graph.h>
7.10 +#include <lemon/dimacs.h>
7.11 +#include <lemon/time_measure.h>
7.12 +//#include <graph_wrapper.h>
7.13 +#include <lemon/preflow.h>
7.14 +#include <augmenting_flow.h>
7.15 +//#include <preflow_res.h>
7.16 +//#include <lp_solver_wrapper_2.h>
7.17 +#include <min_cost_gen_flow.h>
7.18 +
7.19 +// Use a DIMACS max flow file as stdin.
7.20 +// max_flow_demo < dimacs_max_flow_file
7.21 +
7.22 +using namespace lemon;
7.23 +
7.24 +int main(int, char **) {
7.25 +
7.26 + typedef ListGraph MutableGraph;
7.27 + typedef ListGraph Graph;
7.28 + typedef Graph::Node Node;
7.29 + typedef Graph::Edge Edge;
7.30 + typedef Graph::EdgeIt EdgeIt;
7.31 +
7.32 + Graph g;
7.33 +
7.34 + Node s, t;
7.35 + Graph::EdgeMap<int> cap(g);
7.36 + //readDimacsMaxFlow(std::cin, g, s, t, cap);
7.37 + readDimacs(std::cin, g, cap, s, t);
7.38 + Timer ts;
7.39 + Graph::EdgeMap<int> flow(g); //0 flow
7.40 + Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
7.41 + max_flow_test(g, s, t, cap, flow);
7.42 + AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
7.43 + augmenting_flow_test(g, s, t, cap, flow);
7.44 +
7.45 + Graph::NodeMap<bool> cut(g);
7.46 +
7.47 + {
7.48 + std::cout << "preflow ..." << std::endl;
7.49 + ts.reset();
7.50 + max_flow_test.run();
7.51 + std::cout << "elapsed time: " << ts << std::endl;
7.52 + std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
7.53 + max_flow_test.minCut(cut);
7.54 +
7.55 + for (EdgeIt e(g); e!=INVALID; ++e) {
7.56 + if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
7.57 + std::cout << "Slackness does not hold!" << std::endl;
7.58 + if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
7.59 + std::cout << "Slackness does not hold!" << std::endl;
7.60 + }
7.61 + }
7.62 +
7.63 +// {
7.64 +// std::cout << "preflow ..." << std::endl;
7.65 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
7.66 +// ts.reset();
7.67 +// max_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
7.68 +// std::cout << "elapsed time: " << ts << std::endl;
7.69 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
7.70 +
7.71 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
7.72 +// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
7.73 +// std::cout << "Slackness does not hold!" << std::endl;
7.74 +// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
7.75 +// std::cout << "Slackness does not hold!" << std::endl;
7.76 +// }
7.77 +// }
7.78 +
7.79 +// {
7.80 +// std::cout << "wrapped preflow ..." << std::endl;
7.81 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
7.82 +// ts.reset();
7.83 +// pre_flow_res.run();
7.84 +// std::cout << "elapsed time: " << ts << std::endl;
7.85 +// std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
7.86 +// }
7.87 +
7.88 + {
7.89 + std::cout << "physical blocking flow augmentation ..." << std::endl;
7.90 + for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
7.91 + ts.reset();
7.92 + int i=0;
7.93 + while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
7.94 + std::cout << "elapsed time: " << ts << std::endl;
7.95 + std::cout << "number of augmentation phases: " << i << std::endl;
7.96 + std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
7.97 +
7.98 + for (EdgeIt e(g); e!=INVALID; ++e) {
7.99 + if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
7.100 + std::cout << "Slackness does not hold!" << std::endl;
7.101 + if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
7.102 + std::cout << "Slackness does not hold!" << std::endl;
7.103 + }
7.104 + }
7.105 +
7.106 +// {
7.107 +// std::cout << "faster physical blocking flow augmentation ..." << std::endl;
7.108 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
7.109 +// ts.reset();
7.110 +// int i=0;
7.111 +// while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
7.112 +// std::cout << "elapsed time: " << ts << std::endl;
7.113 +// std::cout << "number of augmentation phases: " << i << std::endl;
7.114 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
7.115 +// }
7.116 +
7.117 + {
7.118 + std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
7.119 + for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
7.120 + ts.reset();
7.121 + int i=0;
7.122 + while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
7.123 + std::cout << "elapsed time: " << ts << std::endl;
7.124 + std::cout << "number of augmentation phases: " << i << std::endl;
7.125 + std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
7.126 +
7.127 + for (EdgeIt e(g); e!=INVALID; ++e) {
7.128 + if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
7.129 + std::cout << "Slackness does not hold!" << std::endl;
7.130 + if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
7.131 + std::cout << "Slackness does not hold!" << std::endl;
7.132 + }
7.133 + }
7.134 +
7.135 +// {
7.136 +// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
7.137 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
7.138 +// ts.reset();
7.139 +// int i=0;
7.140 +// while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
7.141 +// std::cout << "elapsed time: " << ts << std::endl;
7.142 +// std::cout << "number of augmentation phases: " << i << std::endl;
7.143 +// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
7.144 +
7.145 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
7.146 +// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
7.147 +// std::cout << "Slackness does not hold!" << std::endl;
7.148 +// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
7.149 +// std::cout << "Slackness does not hold!" << std::endl;
7.150 +// }
7.151 +// }
7.152 +
7.153 +// {
7.154 +// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
7.155 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
7.156 +// ts.reset();
7.157 +// int i=0;
7.158 +// while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
7.159 +// std::cout << "elapsed time: " << ts << std::endl;
7.160 +// std::cout << "number of augmentation phases: " << i << std::endl;
7.161 +// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
7.162 +
7.163 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
7.164 +// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
7.165 +// std::cout << "Slackness does not hold!" << std::endl;
7.166 +// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
7.167 +// std::cout << "Slackness does not hold!" << std::endl;
7.168 +// }
7.169 +// }
7.170 +
7.171 + ts.reset();
7.172 +
7.173 + Edge e=g.addEdge(t, s);
7.174 + Graph::EdgeMap<int> cost(g, 0);
7.175 + cost.set(e, -1);
7.176 + cap.set(e, 10000);
7.177 + typedef ConstMap<Node, int> Excess;
7.178 + Excess excess(0);
7.179 + typedef ConstMap<Edge, int> LCap;
7.180 + LCap lcap(0);
7.181 +
7.182 + MinCostGenFlow<Graph, int, Excess, LCap>
7.183 + min_cost(g, excess, lcap, cap, flow, cost);
7.184 + min_cost.feasible();
7.185 + min_cost.runByLP();
7.186 +
7.187 + std::cout << "elapsed time: " << ts << std::endl;
7.188 + std::cout << "flow value: "<< flow[e] << std::endl;
7.189 +}
8.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
8.2 +++ b/src/work/athos/lp/max_flow_expression.cc Tue Mar 22 11:45:47 2005 +0000
8.3 @@ -0,0 +1,109 @@
8.4 +// -*- c++ -*-
8.5 +#include <iostream>
8.6 +#include <fstream>
8.7 +
8.8 +#include <lemon/graph_utils.h>
8.9 +#include <lemon/smart_graph.h>
8.10 +#include <lemon/list_graph.h>
8.11 +#include <lemon/dimacs.h>
8.12 +#include <lemon/time_measure.h>
8.13 +#include <lp_solver_base.h>
8.14 +
8.15 +using std::cout;
8.16 +using std::endl;
8.17 +using namespace lemon;
8.18 +
8.19 +template<typename Edge, typename EdgeIndexMap>
8.20 +class PrimalMap {
8.21 +protected:
8.22 + LPGLPK* lp;
8.23 + EdgeIndexMap* edge_index_map;
8.24 +public:
8.25 + PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) :
8.26 + lp(&_lp), edge_index_map(&_edge_index_map) { }
8.27 + double operator[](Edge e) const {
8.28 + return lp->getPrimal((*edge_index_map)[e]);
8.29 + }
8.30 +};
8.31 +
8.32 +// Use a DIMACS max flow file as stdin.
8.33 +// max_flow_expresion < dimacs_max_flow_file
8.34 +
8.35 +int main(int, char **) {
8.36 +
8.37 + typedef ListGraph Graph;
8.38 + typedef Graph::Node Node;
8.39 + typedef Graph::Edge Edge;
8.40 + typedef Graph::EdgeIt EdgeIt;
8.41 +
8.42 + Graph g;
8.43 +
8.44 + Node s, t;
8.45 + Graph::EdgeMap<int> cap(g);
8.46 + readDimacs(std::cin, g, cap, s, t);
8.47 + Timer ts;
8.48 +
8.49 + typedef LPGLPK LPSolver;
8.50 + LPSolver lp;
8.51 + lp.setMaximize();
8.52 + typedef LPSolver::Col Col;
8.53 + typedef LPSolver::Row Row;
8.54 + typedef Graph::EdgeMap<Col> EdgeIndexMap;
8.55 + typedef Graph::NodeMap<Row> NodeIndexMap;
8.56 + EdgeIndexMap edge_index_map(g);
8.57 + NodeIndexMap node_index_map(g);
8.58 + PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
8.59 +
8.60 + // nonnegativity of flow and capacity function
8.61 + for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
8.62 + Col col=lp.addCol();
8.63 + edge_index_map.set(e, col);
8.64 + // interesting property in GLPK:
8.65 + // if you change the order of the following two lines, the
8.66 + // two runs of GLPK are extremely different
8.67 + lp.setColLowerBound(col, 0);
8.68 + lp.setColUpperBound(col, cap[e]);
8.69 + }
8.70 +
8.71 + for (Graph::NodeIt n(g); n!=INVALID; ++n) {
8.72 + LPSolver::Expression expr;
8.73 + for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
8.74 + expr+=edge_index_map[e];
8.75 + for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
8.76 + expr-=edge_index_map[e];
8.77 + // cost function
8.78 + if (n==s) {
8.79 + lp.setObjCoeffs(expr);
8.80 + }
8.81 + // flow conservation constraints
8.82 + if ((n!=s) && (n!=t)) {
8.83 + node_index_map.set(n, lp.addRow(expr == 0.0));
8.84 + }
8.85 + }
8.86 + lp.solveSimplex();
8.87 + cout << "elapsed time: " << ts << endl;
8.88 +// cout << "rows:" << endl;
8.89 +// for (
8.90 +// LPSolver::Rows::ClassIt i(lp.row_iter_map, 0);
8.91 +// i!=INVALID;
8.92 +// ++i) {
8.93 +// cout << i << " ";
8.94 +// }
8.95 +// cout << endl;
8.96 +// cout << "cols:" << endl;
8.97 +// for (
8.98 +// LPSolver::Cols::ClassIt i(lp.col_iter_map, 0);
8.99 +// i!=INVALID;
8.100 +// ++i) {
8.101 +// cout << i << " ";
8.102 +// }
8.103 +// cout << endl;
8.104 + lp.setMIP();
8.105 + cout << "elapsed time: " << ts << endl;
8.106 + for (LPSolver::Cols::ClassIt it(lp.col_iter_map ,1); it!=INVALID; ++it) {
8.107 + lp.setColInt(it);
8.108 + }
8.109 + cout << "elapsed time: " << ts << endl;
8.110 + lp.solveBandB();
8.111 + cout << "elapsed time: " << ts << endl;
8.112 +}
9.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
9.2 +++ b/src/work/athos/lp/min_cost_gen_flow.h Tue Mar 22 11:45:47 2005 +0000
9.3 @@ -0,0 +1,268 @@
9.4 +// -*- c++ -*-
9.5 +#ifndef LEMON_MIN_COST_GEN_FLOW_H
9.6 +#define LEMON_MIN_COST_GEN_FLOW_H
9.7 +#include <iostream>
9.8 +//#include <fstream>
9.9 +
9.10 +#include <lemon/smart_graph.h>
9.11 +#include <lemon/list_graph.h>
9.12 +//#include <lemon/dimacs.h>
9.13 +//#include <lemon/time_measure.h>
9.14 +//#include <graph_wrapper.h>
9.15 +#include <lemon/preflow.h>
9.16 +#include <lemon/min_cost_flow.h>
9.17 +//#include <augmenting_flow.h>
9.18 +//#include <preflow_res.h>
9.19 +#include <work/marci/merge_node_graph_wrapper.h>
9.20 +#include <work/marci/lp/lp_solver_wrapper_3.h>
9.21 +
9.22 +namespace lemon {
9.23 +
9.24 + template<typename Edge, typename EdgeIndexMap>
9.25 + class PrimalMap {
9.26 + protected:
9.27 + LPGLPK* lp;
9.28 + EdgeIndexMap* edge_index_map;
9.29 + public:
9.30 + PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) :
9.31 + lp(&_lp), edge_index_map(&_edge_index_map) { }
9.32 + double operator[](Edge e) const {
9.33 + return lp->getPrimal((*edge_index_map)[e]);
9.34 + }
9.35 + };
9.36 +
9.37 + // excess: rho-delta egyelore csak =0-ra.
9.38 + template <typename Graph, typename Num,
9.39 + typename Excess=typename Graph::template NodeMap<Num>,
9.40 + typename LCapMap=typename Graph::template EdgeMap<Num>,
9.41 + typename CapMap=typename Graph::template EdgeMap<Num>,
9.42 + typename FlowMap=typename Graph::template EdgeMap<Num>,
9.43 + typename CostMap=typename Graph::template EdgeMap<Num> >
9.44 + class MinCostGenFlow {
9.45 + protected:
9.46 + const Graph& g;
9.47 + const Excess& excess;
9.48 + const LCapMap& lcapacity;
9.49 + const CapMap& capacity;
9.50 + FlowMap& flow;
9.51 + const CostMap& cost;
9.52 + public:
9.53 + MinCostGenFlow(const Graph& _g, const Excess& _excess,
9.54 + const LCapMap& _lcapacity, const CapMap& _capacity,
9.55 + FlowMap& _flow,
9.56 + const CostMap& _cost) :
9.57 + g(_g), excess(_excess), lcapacity(_lcapacity),
9.58 + capacity(_capacity), flow(_flow), cost(_cost) { }
9.59 + bool feasible() {
9.60 + // std::cout << "making new vertices..." << std::endl;
9.61 + typedef ListGraph Graph2;
9.62 + Graph2 g2;
9.63 + typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
9.64 + // std::cout << "merging..." << std::endl;
9.65 + GW gw(g, g2);
9.66 + typename GW::Node s(INVALID, g2.addNode(), true);
9.67 + typename GW::Node t(INVALID, g2.addNode(), true);
9.68 + typedef SmartGraph Graph3;
9.69 + // std::cout << "making extender graph..." << std::endl;
9.70 + typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
9.71 +// {
9.72 +// checkConcept<StaticGraph, GWW>();
9.73 +// }
9.74 + GWW gww(gw);
9.75 + typedef AugmentingGraphWrapper<GW, GWW> GWWW;
9.76 + GWWW gwww(gw, gww);
9.77 +
9.78 + // std::cout << "making new edges..." << std::endl;
9.79 + typename GWWW::template EdgeMap<Num> translated_cap(gwww);
9.80 +
9.81 + for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
9.82 + translated_cap.set(typename GWWW::Edge(e,INVALID,false),
9.83 + capacity[e]-lcapacity[e]);
9.84 + // cout << "t_cap " << gw.id(e) << " "
9.85 + // << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
9.86 + }
9.87 +
9.88 + Num expected=0;
9.89 +
9.90 + // std::cout << "making new edges 2..." << std::endl;
9.91 + for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
9.92 + Num a=0;
9.93 + for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
9.94 + a+=lcapacity[e];
9.95 + for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
9.96 + a-=lcapacity[e];
9.97 + if (excess[n]>a) {
9.98 + typename GWW::Edge e=
9.99 + gww.addEdge(typename GW::Node(n,INVALID,false), t);
9.100 + translated_cap.set(typename GWWW::Edge(INVALID, e, true),
9.101 + excess[n]-a);
9.102 + // std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
9.103 + }
9.104 + if (excess[n]<a) {
9.105 + typename GWW::Edge e=
9.106 + gww.addEdge(s, typename GW::Node(n,INVALID,false));
9.107 + translated_cap.set(typename GWWW::Edge(INVALID, e, true),
9.108 + a-excess[n]);
9.109 + expected+=a-excess[n];
9.110 + // std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
9.111 + }
9.112 + }
9.113 +
9.114 + // std::cout << "preflow..." << std::endl;
9.115 + typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
9.116 + Preflow<GWWW, Num> preflow(gwww, s, t,
9.117 + translated_cap, translated_flow);
9.118 + preflow.run();
9.119 + // std::cout << "fv: " << preflow.flowValue() << std::endl;
9.120 + // std::cout << "expected: " << expected << std::endl;
9.121 +
9.122 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
9.123 + typename GW::Edge ew(e, INVALID, false);
9.124 + typename GWWW::Edge ewww(ew, INVALID, false);
9.125 + flow.set(e, translated_flow[ewww]+lcapacity[e]);
9.126 + }
9.127 + return (preflow.flowValue()>=expected);
9.128 + }
9.129 + // for nonnegative costs
9.130 + bool run() {
9.131 + // std::cout << "making new vertices..." << std::endl;
9.132 + typedef ListGraph Graph2;
9.133 + Graph2 g2;
9.134 + typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
9.135 + // std::cout << "merging..." << std::endl;
9.136 + GW gw(g, g2);
9.137 + typename GW::Node s(INVALID, g2.addNode(), true);
9.138 + typename GW::Node t(INVALID, g2.addNode(), true);
9.139 + typedef SmartGraph Graph3;
9.140 + // std::cout << "making extender graph..." << std::endl;
9.141 + typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
9.142 +// {
9.143 +// checkConcept<StaticGraph, GWW>();
9.144 +// }
9.145 + GWW gww(gw);
9.146 + typedef AugmentingGraphWrapper<GW, GWW> GWWW;
9.147 + GWWW gwww(gw, gww);
9.148 +
9.149 + // std::cout << "making new edges..." << std::endl;
9.150 + typename GWWW::template EdgeMap<Num> translated_cap(gwww);
9.151 +
9.152 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
9.153 + typename GW::Edge ew(e, INVALID, false);
9.154 + typename GWWW::Edge ewww(ew, INVALID, false);
9.155 + translated_cap.set(ewww, capacity[e]-lcapacity[e]);
9.156 + // cout << "t_cap " << g.id(e) << " "
9.157 + // << translated_cap[ewww] << endl;
9.158 + }
9.159 +
9.160 + Num expected=0;
9.161 +
9.162 + // std::cout << "making new edges 2..." << std::endl;
9.163 + for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
9.164 + // std::cout << "node: " << g.id(n) << std::endl;
9.165 + Num a=0;
9.166 + for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) {
9.167 + a+=lcapacity[e];
9.168 + // std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl;
9.169 + }
9.170 + for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) {
9.171 + a-=lcapacity[e];
9.172 + // std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl;
9.173 + }
9.174 + // std::cout << "excess " << g.id(n) << ": " << a << std::endl;
9.175 + if (0>a) {
9.176 + typename GWW::Edge e=
9.177 + gww.addEdge(typename GW::Node(n,INVALID,false), t);
9.178 + translated_cap.set(typename GWWW::Edge(INVALID, e, true),
9.179 + -a);
9.180 + // std::cout << g.id(n) << "->t " << -a << std::endl;
9.181 + }
9.182 + if (0<a) {
9.183 + typename GWW::Edge e=
9.184 + gww.addEdge(s, typename GW::Node(n,INVALID,false));
9.185 + translated_cap.set(typename GWWW::Edge(INVALID, e, true),
9.186 + a);
9.187 + expected+=a;
9.188 + // std::cout << "s->" << g.id(n) << " "<< a <<std:: endl;
9.189 + }
9.190 + }
9.191 +
9.192 + // std::cout << "preflow..." << std::endl;
9.193 + typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0);
9.194 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
9.195 + translated_cost.set(typename GWWW::Edge(
9.196 + typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]);
9.197 + }
9.198 + // typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
9.199 + MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>,
9.200 + typename GWWW::template EdgeMap<Num> >
9.201 + min_cost_flow(gwww, translated_cost, translated_cap,
9.202 + s, t);
9.203 + while (min_cost_flow.augment()) { }
9.204 + std::cout << "fv: " << min_cost_flow.flowValue() << std::endl;
9.205 + std::cout << "expected: " << expected << std::endl;
9.206 +
9.207 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
9.208 + typename GW::Edge ew(e, INVALID, false);
9.209 + typename GWWW::Edge ewww(ew, INVALID, false);
9.210 + // std::cout << g.id(e) << " " << flow[e] << std::endl;
9.211 + flow.set(e, lcapacity[e]+
9.212 + min_cost_flow.getFlow()[ewww]);
9.213 + }
9.214 + return (min_cost_flow.flowValue()>=expected);
9.215 + }
9.216 + void runByLP() {
9.217 + typedef LPGLPK LPSolver;
9.218 + LPSolver lp;
9.219 + lp.setMinimize();
9.220 + typedef LPSolver::ColIt ColIt;
9.221 + typedef LPSolver::RowIt RowIt;
9.222 + typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
9.223 + EdgeIndexMap edge_index_map(g);
9.224 + PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
9.225 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
9.226 + ColIt col_it=lp.addCol();
9.227 + edge_index_map.set(e, col_it);
9.228 + if (lcapacity[e]==capacity[e])
9.229 + lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]);
9.230 + else
9.231 + lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]);
9.232 + lp.setObjCoef(col_it, cost[e]);
9.233 + }
9.234 + LPSolver::ColIt col_it;
9.235 + for (lp.col_iter_map.first(col_it, lp.VALID_CLASS);
9.236 + lp.col_iter_map.valid(col_it);
9.237 + lp.col_iter_map.next(col_it)) {
9.238 +// std::cout << "ize " << lp.col_iter_map[col_it] << std::endl;
9.239 + }
9.240 + for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
9.241 + typename Graph::template EdgeMap<Num> coeffs(g, 0);
9.242 + for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
9.243 + coeffs.set(e, coeffs[e]+1);
9.244 + for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
9.245 + coeffs.set(e, coeffs[e]-1);
9.246 + RowIt row_it=lp.addRow();
9.247 + typename std::vector< std::pair<ColIt, double> > row;
9.248 + //std::cout << "node:" <<g.id(n)<<std::endl;
9.249 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
9.250 + if (coeffs[e]!=0) {
9.251 + //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
9.252 + row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
9.253 + }
9.254 + }
9.255 + //std::cout << std::endl;
9.256 + //std::cout << " " << g.id(n) << " " << row.size() << std::endl;
9.257 + lp.setRowCoeffs(row_it, row.begin(), row.end());
9.258 + lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
9.259 + }
9.260 + lp.solveSimplex();
9.261 + //std::cout << lp.colNum() << std::endl;
9.262 + //std::cout << lp.rowNum() << std::endl;
9.263 + //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
9.264 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e)
9.265 + flow.set(e, lp_flow[e]);
9.266 + }
9.267 + };
9.268 +
9.269 +} // namespace lemon
9.270 +
9.271 +#endif //LEMON_MIN_COST_GEN_FLOW_H