1.1 --- a/src/work/athos/lp/expression.h Tue Mar 22 16:49:30 2005 +0000
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,197 +0,0 @@
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 --- a/src/work/athos/lp/expression_test.cc Tue Mar 22 16:49:30 2005 +0000
2.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
2.3 @@ -1,23 +0,0 @@
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 --- a/src/work/athos/lp/lp_solver_base.h Tue Mar 22 16:49:30 2005 +0000
3.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
3.3 @@ -1,639 +0,0 @@
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 -
3.19 -#include <iostream>
3.20 -#include <vector>
3.21 -#include <string>
3.22 -#include <list>
3.23 -#include <memory>
3.24 -#include <utility>
3.25 -
3.26 -#include <lemon/invalid.h>
3.27 -#include <expression.h>
3.28 -//#include <stp.h>
3.29 -//#include <lemon/max_flow.h>
3.30 -//#include <augmenting_flow.h>
3.31 -//#include <iter_map.h>
3.32 -
3.33 -using std::cout;
3.34 -using std::cin;
3.35 -using std::endl;
3.36 -
3.37 -namespace lemon {
3.38 -
3.39 - /// \addtogroup misc
3.40 - /// @{
3.41 -
3.42 - /// \brief A partitioned vector with iterable classes.
3.43 - ///
3.44 - /// This class implements a container in which the data is stored in an
3.45 - /// stl vector, the range is partitioned into sets and each set is
3.46 - /// doubly linked in a list.
3.47 - /// That is, each class is iterable by lemon iterators, and any member of
3.48 - /// the vector can bo moved to an other class.
3.49 - template <typename T>
3.50 - class IterablePartition {
3.51 - protected:
3.52 - struct Node {
3.53 - T data;
3.54 - int prev; //invalid az -1
3.55 - int next;
3.56 - };
3.57 - std::vector<Node> nodes;
3.58 - struct Tip {
3.59 - int first;
3.60 - int last;
3.61 - };
3.62 - std::vector<Tip> tips;
3.63 - public:
3.64 - /// The classes are indexed by integers from \c 0 to \c classNum()-1.
3.65 - int classNum() const { return tips.size(); }
3.66 - /// This lemon style iterator iterates through a class.
3.67 - class Class;
3.68 - /// Constructor. The number of classes is to be given which is fixed
3.69 - /// over the life of the container.
3.70 - /// The partition classes are indexed from 0 to class_num-1.
3.71 - IterablePartition(int class_num) {
3.72 - for (int i=0; i<class_num; ++i) {
3.73 - Tip t;
3.74 - t.first=t.last=-1;
3.75 - tips.push_back(t);
3.76 - }
3.77 - }
3.78 - protected:
3.79 - void befuz(Class it, int class_id) {
3.80 - if (tips[class_id].first==-1) {
3.81 - if (tips[class_id].last==-1) {
3.82 - nodes[it.i].prev=nodes[it.i].next=-1;
3.83 - tips[class_id].first=tips[class_id].last=it.i;
3.84 - }
3.85 - } else {
3.86 - nodes[it.i].prev=tips[class_id].last;
3.87 - nodes[it.i].next=-1;
3.88 - nodes[tips[class_id].last].next=it.i;
3.89 - tips[class_id].last=it.i;
3.90 - }
3.91 - }
3.92 - void kifuz(Class it, int class_id) {
3.93 - if (tips[class_id].first==it.i) {
3.94 - if (tips[class_id].last==it.i) {
3.95 - tips[class_id].first=tips[class_id].last=-1;
3.96 - } else {
3.97 - tips[class_id].first=nodes[it.i].next;
3.98 - nodes[nodes[it.i].next].prev=-1;
3.99 - }
3.100 - } else {
3.101 - if (tips[class_id].last==it.i) {
3.102 - tips[class_id].last=nodes[it.i].prev;
3.103 - nodes[nodes[it.i].prev].next=-1;
3.104 - } else {
3.105 - nodes[nodes[it.i].next].prev=nodes[it.i].prev;
3.106 - nodes[nodes[it.i].prev].next=nodes[it.i].next;
3.107 - }
3.108 - }
3.109 - }
3.110 - public:
3.111 - /// A new element with data \c t is pushed into the vector and into class
3.112 - /// \c class_id.
3.113 - Class push_back(const T& t, int class_id) {
3.114 - Node n;
3.115 - n.data=t;
3.116 - nodes.push_back(n);
3.117 - int i=nodes.size()-1;
3.118 - befuz(i, class_id);
3.119 - return i;
3.120 - }
3.121 - /// A member is moved to an other class.
3.122 - void set(Class it, int old_class_id, int new_class_id) {
3.123 - kifuz(it.i, old_class_id);
3.124 - befuz(it.i, new_class_id);
3.125 - }
3.126 - /// Returns the data pointed by \c it.
3.127 - T& operator[](Class it) { return nodes[it.i].data; }
3.128 - /// Returns the data pointed by \c it.
3.129 - const T& operator[](Class it) const { return nodes[it.i].data; }
3.130 - ///.
3.131 - class Class {
3.132 - friend class IterablePartition;
3.133 - protected:
3.134 - int i;
3.135 - public:
3.136 - /// Default constructor.
3.137 - Class() { }
3.138 - /// This constructor constructs an iterator which points
3.139 - /// to the member of th container indexed by the integer _i.
3.140 - Class(const int& _i) : i(_i) { }
3.141 - /// Invalid constructor.
3.142 - Class(const Invalid&) : i(-1) { }
3.143 - friend bool operator<(const Class& x, const Class& y);
3.144 - friend std::ostream& operator<<(std::ostream& os,
3.145 - const Class& it);
3.146 - bool operator==(const Class& node) const {return i == node.i;}
3.147 - bool operator!=(const Class& node) const {return i != node.i;}
3.148 - };
3.149 - friend bool operator<(const Class& x, const Class& y) {
3.150 - return (x.i < y.i);
3.151 - }
3.152 - friend std::ostream& operator<<(std::ostream& os,
3.153 - const Class& it) {
3.154 - os << it.i;
3.155 - return os;
3.156 - }
3.157 - /// First member of class \c class_id.
3.158 - Class& first(Class& it, int class_id) const {
3.159 - it.i=tips[class_id].first;
3.160 - return it;
3.161 - }
3.162 - /// Next member.
3.163 - Class& next(Class& it) const {
3.164 - it.i=nodes[it.i].next;
3.165 - return it;
3.166 - }
3.167 - /// True iff the iterator is valid.
3.168 - bool valid(const Class& it) const { return it.i!=-1; }
3.169 -
3.170 - class ClassIt : public Class {
3.171 - const IterablePartition* iterable_partition;
3.172 - public:
3.173 - ClassIt() { }
3.174 - ClassIt(Invalid i) : Class(i) { }
3.175 - ClassIt(const IterablePartition& _iterable_partition,
3.176 - const int& i) : iterable_partition(&_iterable_partition) {
3.177 - _iterable_partition.first(*this, i);
3.178 - }
3.179 - ClassIt(const IterablePartition& _iterable_partition,
3.180 - const Class& _class) :
3.181 - Class(_class), iterable_partition(&_iterable_partition) { }
3.182 - ClassIt& operator++() {
3.183 - iterable_partition->next(*this);
3.184 - return *this;
3.185 - }
3.186 - };
3.187 -
3.188 - };
3.189 -
3.190 -
3.191 - /*! \e
3.192 - \todo kellenene uj iterable structure bele, mert ez nem az igazi
3.193 - \todo A[x,y]-t cserel. Jobboldal, baloldal csere.
3.194 - \todo LEKERDEZESEK!!!
3.195 - \todo DOKSI!!!! Doxygen group!!!
3.196 - The aim of this class is to give a general surface to different
3.197 - solvers, i.e. it makes possible to write algorithms using LP's,
3.198 - in which the solver can be changed to an other one easily.
3.199 - \nosubgrouping
3.200 - */
3.201 - template <typename _Value>
3.202 - class LpSolverBase {
3.203 -
3.204 - /*! @name Uncategorized functions and types (public members)
3.205 - */
3.206 - //@{
3.207 - public:
3.208 -
3.209 - //UNCATEGORIZED
3.210 -
3.211 - /// \e
3.212 - typedef IterablePartition<int> Rows;
3.213 - /// \e
3.214 - typedef IterablePartition<int> Cols;
3.215 - /// \e
3.216 - typedef _Value Value;
3.217 - /// \e
3.218 - typedef Rows::Class Row;
3.219 - /// \e
3.220 - typedef Cols::Class Col;
3.221 - public:
3.222 - /// \e
3.223 - IterablePartition<int> row_iter_map;
3.224 - /// \e
3.225 - IterablePartition<int> col_iter_map;
3.226 - /// \e
3.227 - std::vector<Row> int_row_map;
3.228 - /// \e
3.229 - std::vector<Col> int_col_map;
3.230 - /// \e
3.231 - const int VALID_CLASS;
3.232 - /// \e
3.233 - const int INVALID_CLASS;
3.234 - /// \e
3.235 - static const Value INF;
3.236 - public:
3.237 - /// \e
3.238 - LpSolverBase() : row_iter_map(2),
3.239 - col_iter_map(2),
3.240 - VALID_CLASS(0), INVALID_CLASS(1) { }
3.241 - /// \e
3.242 - virtual ~LpSolverBase() { }
3.243 - //@}
3.244 -
3.245 - /*! @name Medium level interface (public members)
3.246 - These functions appear in the low level and also in the high level
3.247 - interfaces thus these each of these functions have to be implemented
3.248 - only once in the different interfaces.
3.249 - This means that these functions have to be reimplemented for all of the
3.250 - different lp solvers. These are basic functions, and have the same
3.251 - parameter lists in the low and high level interfaces.
3.252 - */
3.253 - //@{
3.254 - public:
3.255 -
3.256 - //UNCATEGORIZED FUNCTIONS
3.257 -
3.258 - /// \e
3.259 - virtual void setMinimize() = 0;
3.260 - /// \e
3.261 - virtual void setMaximize() = 0;
3.262 -
3.263 - //SOLVER FUNCTIONS
3.264 -
3.265 - /// \e
3.266 - virtual void solveSimplex() = 0;
3.267 - /// \e
3.268 - virtual void solvePrimalSimplex() = 0;
3.269 - /// \e
3.270 - virtual void solveDualSimplex() = 0;
3.271 -
3.272 - //SOLUTION RETRIEVING
3.273 -
3.274 - /// \e
3.275 - virtual Value getObjVal() = 0;
3.276 -
3.277 - //OTHER FUNCTIONS
3.278 -
3.279 - /// \e
3.280 - virtual int rowNum() const = 0;
3.281 - /// \e
3.282 - virtual int colNum() const = 0;
3.283 - /// \e
3.284 - virtual int warmUp() = 0;
3.285 - /// \e
3.286 - virtual void printWarmUpStatus(int i) = 0;
3.287 - /// \e
3.288 - virtual int getPrimalStatus() = 0;
3.289 - /// \e
3.290 - virtual void printPrimalStatus(int i) = 0;
3.291 - /// \e
3.292 - virtual int getDualStatus() = 0;
3.293 - /// \e
3.294 - virtual void printDualStatus(int i) = 0;
3.295 - /// Returns the status of the slack variable assigned to row \c row.
3.296 - virtual int getRowStat(const Row& row) = 0;
3.297 - /// \e
3.298 - virtual void printRowStatus(int i) = 0;
3.299 - /// Returns the status of the variable assigned to column \c col.
3.300 - virtual int getColStat(const Col& col) = 0;
3.301 - /// \e
3.302 - virtual void printColStatus(int i) = 0;
3.303 -
3.304 - //@}
3.305 -
3.306 - /*! @name Low level interface (protected members)
3.307 - Problem manipulating functions in the low level interface
3.308 - */
3.309 - //@{
3.310 - protected:
3.311 -
3.312 - //MATRIX MANIPULATING FUNCTIONS
3.313 -
3.314 - /// \e
3.315 - virtual int _addCol() = 0;
3.316 - /// \e
3.317 - virtual int _addRow() = 0;
3.318 - /// \e
3.319 - virtual void _eraseCol(int i) = 0;
3.320 - /// \e
3.321 - virtual void _eraseRow(int i) = 0;
3.322 - /// \e
3.323 - virtual void _setRowCoeffs(int i,
3.324 - const std::vector<std::pair<int, Value> >& coeffs) = 0;
3.325 - /// \e
3.326 - /// This routine modifies \c coeffs only by the \c push_back method.
3.327 - virtual void _getRowCoeffs(int i,
3.328 - std::vector<std::pair<int, Value> >& coeffs) = 0;
3.329 - /// \e
3.330 - virtual void _setColCoeffs(int i,
3.331 - const std::vector<std::pair<int, Value> >& coeffs) = 0;
3.332 - /// \e
3.333 - /// This routine modifies \c coeffs only by the \c push_back method.
3.334 - virtual void _getColCoeffs(int i,
3.335 - std::vector<std::pair<int, Value> >& coeffs) = 0;
3.336 - /// \e
3.337 - virtual void _setCoeff(int col, int row, Value value) = 0;
3.338 - /// \e
3.339 - virtual Value _getCoeff(int col, int row) = 0;
3.340 - // public:
3.341 - // /// \e
3.342 - // enum Bound { FREE, LOWER, UPPER, DOUBLE, FIXED };
3.343 - protected:
3.344 - /// \e
3.345 - /// The lower bound of a variable (column) have to be given by an
3.346 - /// extended number of type Value, i.e. a finite number of type
3.347 - /// Value or -INF.
3.348 - virtual void _setColLowerBound(int i, Value value) = 0;
3.349 - /// \e
3.350 - /// The lower bound of a variable (column) is an
3.351 - /// extended number of type Value, i.e. a finite number of type
3.352 - /// Value or -INF.
3.353 - virtual Value _getColLowerBound(int i) = 0;
3.354 - /// \e
3.355 - /// The upper bound of a variable (column) have to be given by an
3.356 - /// extended number of type Value, i.e. a finite number of type
3.357 - /// Value or INF.
3.358 - virtual void _setColUpperBound(int i, Value value) = 0;
3.359 - /// \e
3.360 - /// The upper bound of a variable (column) is an
3.361 - /// extended number of type Value, i.e. a finite number of type
3.362 - /// Value or INF.
3.363 - virtual Value _getColUpperBound(int i) = 0;
3.364 - /// \e
3.365 - /// The lower bound of a linear expression (row) have to be given by an
3.366 - /// extended number of type Value, i.e. a finite number of type
3.367 - /// Value or -INF.
3.368 - virtual void _setRowLowerBound(int i, Value value) = 0;
3.369 - /// \e
3.370 - /// The lower bound of a linear expression (row) is an
3.371 - /// extended number of type Value, i.e. a finite number of type
3.372 - /// Value or -INF.
3.373 - virtual Value _getRowLowerBound(int i) = 0;
3.374 - /// \e
3.375 - /// The upper bound of a linear expression (row) have to be given by an
3.376 - /// extended number of type Value, i.e. a finite number of type
3.377 - /// Value or INF.
3.378 - virtual void _setRowUpperBound(int i, Value value) = 0;
3.379 - /// \e
3.380 - /// The upper bound of a linear expression (row) is an
3.381 - /// extended number of type Value, i.e. a finite number of type
3.382 - /// Value or INF.
3.383 - virtual Value _getRowUpperBound(int i) = 0;
3.384 - /// \e
3.385 - virtual void _setObjCoeff(int i, Value obj_coef) = 0;
3.386 - /// \e
3.387 - virtual Value _getObjCoeff(int i) = 0;
3.388 -
3.389 - //SOLUTION RETRIEVING
3.390 -
3.391 - /// \e
3.392 - virtual Value _getPrimal(int i) = 0;
3.393 - //@}
3.394 -
3.395 - /*! @name High level interface (public members)
3.396 - Problem manipulating functions in the high level interface
3.397 - */
3.398 - //@{
3.399 - public:
3.400 -
3.401 - //MATRIX MANIPULATING FUNCTIONS
3.402 -
3.403 - /// \e
3.404 - Col addCol() {
3.405 - int i=_addCol();
3.406 - Col col;
3.407 - col_iter_map.first(col, INVALID_CLASS);
3.408 - if (col_iter_map.valid(col)) { //van hasznalhato hely
3.409 - col_iter_map.set(col, INVALID_CLASS, VALID_CLASS);
3.410 - col_iter_map[col]=i;
3.411 - } else { //a cucc vegere kell inzertalni mert nincs szabad hely
3.412 - col=col_iter_map.push_back(i, VALID_CLASS);
3.413 - }
3.414 - int_col_map.push_back(col);
3.415 - return col;
3.416 - }
3.417 - /// \e
3.418 - Row addRow() {
3.419 - int i=_addRow();
3.420 - Row row;
3.421 - row_iter_map.first(row, INVALID_CLASS);
3.422 - if (row_iter_map.valid(row)) { //van hasznalhato hely
3.423 - row_iter_map.set(row, INVALID_CLASS, VALID_CLASS);
3.424 - row_iter_map[row]=i;
3.425 - } else { //a cucc vegere kell inzertalni mert nincs szabad hely
3.426 - row=row_iter_map.push_back(i, VALID_CLASS);
3.427 - }
3.428 - int_row_map.push_back(row);
3.429 - return row;
3.430 - }
3.431 - /// \e
3.432 - void eraseCol(const Col& col) {
3.433 - col_iter_map.set(col, VALID_CLASS, INVALID_CLASS);
3.434 - int cols[2];
3.435 - cols[1]=col_iter_map[col];
3.436 - _eraseCol(cols[1]);
3.437 - col_iter_map[col]=0; //glpk specifikus, de kell ez??
3.438 - Col it;
3.439 - for (col_iter_map.first(it, VALID_CLASS);
3.440 - col_iter_map.valid(it); col_iter_map.next(it)) {
3.441 - if (col_iter_map[it]>cols[1]) --col_iter_map[it];
3.442 - }
3.443 - int_col_map.erase(int_col_map.begin()+cols[1]);
3.444 - }
3.445 - /// \e
3.446 - void eraseRow(const Row& row) {
3.447 - row_iter_map.set(row, VALID_CLASS, INVALID_CLASS);
3.448 - int rows[2];
3.449 - rows[1]=row_iter_map[row];
3.450 - _eraseRow(rows[1]);
3.451 - row_iter_map[row]=0; //glpk specifikus, de kell ez??
3.452 - Row it;
3.453 - for (row_iter_map.first(it, VALID_CLASS);
3.454 - row_iter_map.valid(it); row_iter_map.next(it)) {
3.455 - if (row_iter_map[it]>rows[1]) --row_iter_map[it];
3.456 - }
3.457 - int_row_map.erase(int_row_map.begin()+rows[1]);
3.458 - }
3.459 - /// \e
3.460 - void setCoeff(Col col, Row row, Value value) {
3.461 - _setCoeff(col_iter_map[col], row_iter_map[row], value);
3.462 - }
3.463 - /// \e
3.464 - Value getCoeff(Col col, Row row) {
3.465 - return _getCoeff(col_iter_map[col], row_iter_map[row], value);
3.466 - }
3.467 - /// \e
3.468 - void setColLowerBound(Col col, Value lo) {
3.469 - _setColLowerBound(col_iter_map[col], lo);
3.470 - }
3.471 - /// \e
3.472 - Value getColLowerBound(Col col) {
3.473 - return _getColLowerBound(col_iter_map[col]);
3.474 - }
3.475 - /// \e
3.476 - void setColUpperBound(Col col, Value up) {
3.477 - _setColUpperBound(col_iter_map[col], up);
3.478 - }
3.479 - /// \e
3.480 - Value getColUpperBound(Col col) {
3.481 - return _getColUpperBound(col_iter_map[col]);
3.482 - }
3.483 - /// \e
3.484 - void setRowLowerBound(Row row, Value lo) {
3.485 - _setRowLowerBound(row_iter_map[row], lo);
3.486 - }
3.487 - /// \e
3.488 - Value getRowLowerBound(Row row) {
3.489 - return _getRowLowerBound(row_iter_map[row]);
3.490 - }
3.491 - /// \e
3.492 - void setRowUpperBound(Row row, Value up) {
3.493 - _setRowUpperBound(row_iter_map[row], up);
3.494 - }
3.495 - /// \e
3.496 - Value getRowUpperBound(Row row) {
3.497 - return _getRowUpperBound(row_iter_map[row]);
3.498 - }
3.499 - /// \e
3.500 - void setObjCoeff(const Col& col, Value obj_coef) {
3.501 - _setObjCoeff(col_iter_map[col], obj_coef);
3.502 - }
3.503 - /// \e
3.504 - Value getObjCoeff(const Col& col) {
3.505 - return _getObjCoeff(col_iter_map[col]);
3.506 - }
3.507 -
3.508 - //SOLUTION RETRIEVING FUNCTIONS
3.509 -
3.510 - /// \e
3.511 - Value getPrimal(const Col& col) {
3.512 - return _getPrimal(col_iter_map[col]);
3.513 - }
3.514 -
3.515 - //@}
3.516 -
3.517 - /*! @name User friend interface
3.518 - Problem manipulating functions in the user friend interface
3.519 - */
3.520 - //@{
3.521 -
3.522 - //EXPRESSION TYPES
3.523 -
3.524 - /// \e
3.525 - typedef Expr<Col, Value> Expression;
3.526 - /// \e
3.527 - typedef Expr<Row, Value> DualExpression;
3.528 - /// \e
3.529 - typedef Constr<Col, Value> Constraint;
3.530 -
3.531 - //MATRIX MANIPULATING FUNCTIONS
3.532 -
3.533 - /// \e
3.534 - void setRowCoeffs(Row row, const Expression& expr) {
3.535 - std::vector<std::pair<int, Value> > row_coeffs;
3.536 - for(typename Expression::Data::const_iterator i=expr.data.begin();
3.537 - i!=expr.data.end(); ++i) {
3.538 - row_coeffs.push_back(std::make_pair
3.539 - (col_iter_map[(*i).first], (*i).second));
3.540 - }
3.541 - _setRowCoeffs(row_iter_map[row], row_coeffs);
3.542 - }
3.543 - /// \e
3.544 - void setRow(Row row, const Constraint& constr) {
3.545 - setRowCoeffs(row, constr.expr);
3.546 - setRowLowerBound(row, constr.lo);
3.547 - setRowUpperBound(row, constr.up);
3.548 - }
3.549 - /// \e
3.550 - Row addRow(const Constraint& constr) {
3.551 - Row row=addRow();
3.552 - setRowCoeffs(row, constr.expr);
3.553 - setRowLowerBound(row, constr.lo);
3.554 - setRowUpperBound(row, constr.up);
3.555 - return row;
3.556 - }
3.557 - /// \e
3.558 - /// This routine modifies \c expr by only adding to it.
3.559 - void getRowCoeffs(Row row, Expression& expr) {
3.560 - std::vector<std::pair<int, Value> > row_coeffs;
3.561 - _getRowCoeffs(row_iter_map[row], row_coeffs);
3.562 - for(typename std::vector<std::pair<int, Value> >::const_iterator
3.563 - i=row_coeffs.begin(); i!=row_coeffs.end(); ++i) {
3.564 - expr+= (*i).second*int_col_map[(*i).first];
3.565 - }
3.566 - }
3.567 - /// \e
3.568 - void setColCoeffs(Col col, const DualExpression& expr) {
3.569 - std::vector<std::pair<int, Value> > col_coeffs;
3.570 - for(typename DualExpression::Data::const_iterator i=expr.data.begin();
3.571 - i!=expr.data.end(); ++i) {
3.572 - col_coeffs.push_back(std::make_pair
3.573 - (row_iter_map[(*i).first], (*i).second));
3.574 - }
3.575 - _setColCoeffs(col_iter_map[col], col_coeffs);
3.576 - }
3.577 - /// \e
3.578 - /// This routine modifies \c expr by only adding to it.
3.579 - void getColCoeffs(Col col, DualExpression& expr) {
3.580 - std::vector<std::pair<int, Value> > col_coeffs;
3.581 - _getColCoeffs(col_iter_map[col], col_coeffs);
3.582 - for(typename std::vector<std::pair<int, Value> >::const_iterator
3.583 - i=col_coeffs.begin(); i!=col_coeffs.end(); ++i) {
3.584 - expr+= (*i).second*int_row_map[(*i).first];
3.585 - }
3.586 - }
3.587 - /// \e
3.588 - void setObjCoeffs(const Expression& expr) {
3.589 - // writing zero everywhere
3.590 - for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
3.591 - setObjCoeff(it, 0.0);
3.592 - // writing the data needed
3.593 - for(typename Expression::Data::const_iterator i=expr.data.begin();
3.594 - i!=expr.data.end(); ++i) {
3.595 - setObjCoeff((*i).first, (*i).second);
3.596 - }
3.597 - }
3.598 - /// \e
3.599 - /// This routine modifies \c expr by only adding to it.
3.600 - void getObjCoeffs(Expression& expr) {
3.601 - for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
3.602 - expr+=getObjCoeff(it)*it;
3.603 - }
3.604 - //@}
3.605 -
3.606 -
3.607 - /*! @name MIP functions and types (public members)
3.608 - */
3.609 - //@{
3.610 - public:
3.611 - /// \e
3.612 - virtual void solveBandB() = 0;
3.613 - /// \e
3.614 - virtual void setLP() = 0;
3.615 - /// \e
3.616 - virtual void setMIP() = 0;
3.617 - protected:
3.618 - /// \e
3.619 - virtual void _setColCont(int i) = 0;
3.620 - /// \e
3.621 - virtual void _setColInt(int i) = 0;
3.622 - /// \e
3.623 - virtual Value _getMIPPrimal(int i) = 0;
3.624 - public:
3.625 - /// \e
3.626 - void setColCont(Col col) {
3.627 - _setColCont(col_iter_map[col]);
3.628 - }
3.629 - /// \e
3.630 - void setColInt(Col col) {
3.631 - _setColInt(col_iter_map[col]);
3.632 - }
3.633 - /// \e
3.634 - Value getMIPPrimal(Col col) {
3.635 - return _getMIPPrimal(col_iter_map[col]);
3.636 - }
3.637 - //@}
3.638 - };
3.639 -
3.640 -} //namespace lemon
3.641 -
3.642 -#endif //LEMON_LP_SOLVER_BASE_H
4.1 --- a/src/work/athos/lp/lp_solver_glpk.h Tue Mar 22 16:49:30 2005 +0000
4.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
4.3 @@ -1,545 +0,0 @@
4.4 -// -*- c++ -*-
4.5 -#ifndef LEMON_LP_SOLVER_GLPK_H
4.6 -#define LEMON_LP_SOLVER_GLPK_H
4.7 -
4.8 -///\ingroup misc
4.9 -///\file
4.10 -
4.11 -// #include <stdio.h>
4.12 -/* #include <stdlib.h> */
4.13 -/* #include <iostream> */
4.14 -/* #include <map> */
4.15 -/* #include <limits> */
4.16 -// #include <stdio>
4.17 -//#include <stdlib>
4.18 -extern "C" {
4.19 -#include "glpk.h"
4.20 -}
4.21 -
4.22 -/* #include <iostream> */
4.23 -/* #include <vector> */
4.24 -/* #include <string> */
4.25 -/* #include <list> */
4.26 -/* #include <memory> */
4.27 -/* #include <utility> */
4.28 -
4.29 -//#include <lemon/invalid.h>
4.30 -//#include <expression.h>
4.31 -#include <lp_solver_base.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 - template <typename Value>
4.45 - const Value LpSolverBase<Value>::INF=std::numeric_limits<Value>::infinity();
4.46 -
4.47 -
4.48 - /// \brief Wrapper for GLPK solver
4.49 - ///
4.50 - /// This class implements a lemon wrapper for GLPK.
4.51 - class LpGlpk : public LpSolverBase<double> {
4.52 - public:
4.53 - typedef LpSolverBase<double> Parent;
4.54 -
4.55 - public:
4.56 - /// \e
4.57 - LPX* lp;
4.58 -
4.59 - public:
4.60 - /// \e
4.61 - LpGlpk() : Parent(),
4.62 - lp(lpx_create_prob()) {
4.63 - int_row_map.push_back(Row());
4.64 - int_col_map.push_back(Col());
4.65 - lpx_set_int_parm(lp, LPX_K_DUAL, 1);
4.66 - }
4.67 - /// \e
4.68 - ~LpGlpk() {
4.69 - lpx_delete_prob(lp);
4.70 - }
4.71 -
4.72 - //MATRIX INDEPEDENT MANIPULATING FUNCTIONS
4.73 -
4.74 - /// \e
4.75 - void setMinimize() {
4.76 - lpx_set_obj_dir(lp, LPX_MIN);
4.77 - }
4.78 - /// \e
4.79 - void setMaximize() {
4.80 - lpx_set_obj_dir(lp, LPX_MAX);
4.81 - }
4.82 -
4.83 - //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS
4.84 -
4.85 - protected:
4.86 - /// \e
4.87 - int _addCol() {
4.88 - int i=lpx_add_cols(lp, 1);
4.89 - _setColLowerBound(i, -INF);
4.90 - _setColUpperBound(i, INF);
4.91 - return i;
4.92 - }
4.93 - /// \e
4.94 - int _addRow() {
4.95 - int i=lpx_add_rows(lp, 1);
4.96 - return i;
4.97 - }
4.98 - /// \e
4.99 - virtual void _setRowCoeffs(int i,
4.100 - const std::vector<std::pair<int, double> >& coeffs) {
4.101 - int mem_length=1+colNum();
4.102 - int* indices = new int[mem_length];
4.103 - double* doubles = new double[mem_length];
4.104 - int length=0;
4.105 - for (std::vector<std::pair<int, double> >::
4.106 - const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
4.107 - ++length;
4.108 - indices[length]=it->first;
4.109 - doubles[length]=it->second;
4.110 - }
4.111 - lpx_set_mat_row(lp, i, length, indices, doubles);
4.112 - delete [] indices;
4.113 - delete [] doubles;
4.114 - }
4.115 - /// \e
4.116 - virtual void _getRowCoeffs(int i,
4.117 - std::vector<std::pair<int, double> >& coeffs) {
4.118 - int mem_length=1+colNum();
4.119 - int* indices = new int[mem_length];
4.120 - double* doubles = new double[mem_length];
4.121 - int length=lpx_get_mat_row(lp, i, indices, doubles);
4.122 - for (int i=1; i<=length; ++i) {
4.123 - coeffs.push_back(std::make_pair(indices[i], doubles[i]));
4.124 - }
4.125 - delete [] indices;
4.126 - delete [] doubles;
4.127 - }
4.128 - /// \e
4.129 - virtual void _setColCoeffs(int i,
4.130 - const std::vector<std::pair<int, double> >& coeffs) {
4.131 - int mem_length=1+rowNum();
4.132 - int* indices = new int[mem_length];
4.133 - double* doubles = new double[mem_length];
4.134 - int length=0;
4.135 - for (std::vector<std::pair<int, double> >::
4.136 - const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
4.137 - ++length;
4.138 - indices[length]=it->first;
4.139 - doubles[length]=it->second;
4.140 - }
4.141 - lpx_set_mat_col(lp, i, length, indices, doubles);
4.142 - delete [] indices;
4.143 - delete [] doubles;
4.144 - }
4.145 - /// \e
4.146 - virtual void _getColCoeffs(int i,
4.147 - std::vector<std::pair<int, double> >& coeffs) {
4.148 - int mem_length=1+rowNum();
4.149 - int* indices = new int[mem_length];
4.150 - double* doubles = new double[mem_length];
4.151 - int length=lpx_get_mat_col(lp, i, indices, doubles);
4.152 - for (int i=1; i<=length; ++i) {
4.153 - coeffs.push_back(std::make_pair(indices[i], doubles[i]));
4.154 - }
4.155 - delete [] indices;
4.156 - delete [] doubles;
4.157 - }
4.158 - /// \e
4.159 - virtual void _eraseCol(int i) {
4.160 - int cols[2];
4.161 - cols[1]=i;
4.162 - lpx_del_cols(lp, 1, cols);
4.163 - }
4.164 - virtual void _eraseRow(int i) {
4.165 - int rows[2];
4.166 - rows[1]=i;
4.167 - lpx_del_rows(lp, 1, rows);
4.168 - }
4.169 - void _setCoeff(int col, int row, double value) {
4.170 - ///FIXME Of course this is not efficient at all, but GLPK knows not more.
4.171 - int change_this;
4.172 - bool get_set_row;
4.173 - //The only thing we can do is optimize on whether working with a row
4.174 - //or a coloumn
4.175 - int row_num = rowNum();
4.176 - int col_num = colNum();
4.177 - if (col_num<row_num){
4.178 - //There are more rows than coloumns
4.179 - get_set_row=true;
4.180 - int mem_length=1+row_num;
4.181 - int* indices = new int[mem_length];
4.182 - double* doubles = new double[mem_length];
4.183 - int length=lpx_get_mat_col(lp, i, indices, doubles);
4.184 - }else{
4.185 - get_set_row=false;
4.186 - int mem_length=1+col_num;
4.187 - int* indices = new int[mem_length];
4.188 - double* doubles = new double[mem_length];
4.189 - int length=lpx_get_mat_row(lp, i, indices, doubles);
4.190 - }
4.191 - //Itten
4.192 -int* indices = new int[mem_length];
4.193 - double* doubles = new double[mem_length];
4.194 - int length=lpx_get_mat_col(lp, i, indices, doubles);
4.195 -
4.196 - delete [] indices;
4.197 - delete [] doubles;
4.198 -
4.199 - }
4.200 - double _getCoeff(int col, int row) {
4.201 - /// FIXME not yet implemented
4.202 - return 0.0;
4.203 - }
4.204 - virtual void _setColLowerBound(int i, double lo) {
4.205 - if (lo==INF) {
4.206 - //FIXME error
4.207 - }
4.208 - int b=lpx_get_col_type(lp, i);
4.209 - double up=lpx_get_col_ub(lp, i);
4.210 - if (lo==-INF) {
4.211 - switch (b) {
4.212 - case LPX_FR:
4.213 - case LPX_LO:
4.214 - lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
4.215 - break;
4.216 - case LPX_UP:
4.217 - break;
4.218 - case LPX_DB:
4.219 - case LPX_FX:
4.220 - lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
4.221 - break;
4.222 - default: ;
4.223 - //FIXME error
4.224 - }
4.225 - } else {
4.226 - switch (b) {
4.227 - case LPX_FR:
4.228 - case LPX_LO:
4.229 - lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
4.230 - break;
4.231 - case LPX_UP:
4.232 - case LPX_DB:
4.233 - case LPX_FX:
4.234 - if (lo==up)
4.235 - lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
4.236 - else
4.237 - lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
4.238 - break;
4.239 - default: ;
4.240 - //FIXME error
4.241 - }
4.242 - }
4.243 - }
4.244 - virtual double _getColLowerBound(int i) {
4.245 - int b=lpx_get_col_type(lp, i);
4.246 - switch (b) {
4.247 - case LPX_FR:
4.248 - return -INF;
4.249 - case LPX_LO:
4.250 - return lpx_get_col_lb(lp, i);
4.251 - case LPX_UP:
4.252 - return -INF;
4.253 - case LPX_DB:
4.254 - case LPX_FX:
4.255 - return lpx_get_col_lb(lp, i);
4.256 - default: ;
4.257 - //FIXME error
4.258 - return 0.0;
4.259 - }
4.260 - }
4.261 - virtual void _setColUpperBound(int i, double up) {
4.262 - if (up==-INF) {
4.263 - //FIXME error
4.264 - }
4.265 - int b=lpx_get_col_type(lp, i);
4.266 - double lo=lpx_get_col_lb(lp, i);
4.267 - if (up==INF) {
4.268 - switch (b) {
4.269 - case LPX_FR:
4.270 - case LPX_LO:
4.271 - break;
4.272 - case LPX_UP:
4.273 - lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
4.274 - break;
4.275 - case LPX_DB:
4.276 - case LPX_FX:
4.277 - lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
4.278 - break;
4.279 - default: ;
4.280 - //FIXME error
4.281 - }
4.282 - } else {
4.283 - switch (b) {
4.284 - case LPX_FR:
4.285 - lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
4.286 - case LPX_LO:
4.287 - if (lo==up)
4.288 - lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
4.289 - else
4.290 - lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
4.291 - break;
4.292 - case LPX_UP:
4.293 - lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
4.294 - break;
4.295 - case LPX_DB:
4.296 - case LPX_FX:
4.297 - if (lo==up)
4.298 - lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
4.299 - else
4.300 - lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
4.301 - break;
4.302 - default: ;
4.303 - //FIXME error
4.304 - }
4.305 - }
4.306 - }
4.307 - virtual double _getColUpperBound(int i) {
4.308 - int b=lpx_get_col_type(lp, i);
4.309 - switch (b) {
4.310 - case LPX_FR:
4.311 - case LPX_LO:
4.312 - return INF;
4.313 - case LPX_UP:
4.314 - case LPX_DB:
4.315 - case LPX_FX:
4.316 - return lpx_get_col_ub(lp, i);
4.317 - default: ;
4.318 - //FIXME error
4.319 - return 0.0;
4.320 - }
4.321 - }
4.322 - virtual void _setRowLowerBound(int i, double lo) {
4.323 - if (lo==INF) {
4.324 - //FIXME error
4.325 - }
4.326 - int b=lpx_get_row_type(lp, i);
4.327 - double up=lpx_get_row_ub(lp, i);
4.328 - if (lo==-INF) {
4.329 - switch (b) {
4.330 - case LPX_FR:
4.331 - case LPX_LO:
4.332 - lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
4.333 - break;
4.334 - case LPX_UP:
4.335 - break;
4.336 - case LPX_DB:
4.337 - case LPX_FX:
4.338 - lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
4.339 - break;
4.340 - default: ;
4.341 - //FIXME error
4.342 - }
4.343 - } else {
4.344 - switch (b) {
4.345 - case LPX_FR:
4.346 - case LPX_LO:
4.347 - lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
4.348 - break;
4.349 - case LPX_UP:
4.350 - case LPX_DB:
4.351 - case LPX_FX:
4.352 - if (lo==up)
4.353 - lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
4.354 - else
4.355 - lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
4.356 - break;
4.357 - default: ;
4.358 - //FIXME error
4.359 - }
4.360 - }
4.361 - }
4.362 - virtual double _getRowLowerBound(int i) {
4.363 - int b=lpx_get_row_type(lp, i);
4.364 - switch (b) {
4.365 - case LPX_FR:
4.366 - return -INF;
4.367 - case LPX_LO:
4.368 - return lpx_get_row_lb(lp, i);
4.369 - case LPX_UP:
4.370 - return -INF;
4.371 - case LPX_DB:
4.372 - case LPX_FX:
4.373 - return lpx_get_row_lb(lp, i);
4.374 - default: ;
4.375 - //FIXME error
4.376 - return 0.0;
4.377 - }
4.378 - }
4.379 - virtual void _setRowUpperBound(int i, double up) {
4.380 - if (up==-INF) {
4.381 - //FIXME error
4.382 - }
4.383 - int b=lpx_get_row_type(lp, i);
4.384 - double lo=lpx_get_row_lb(lp, i);
4.385 - if (up==INF) {
4.386 - switch (b) {
4.387 - case LPX_FR:
4.388 - case LPX_LO:
4.389 - break;
4.390 - case LPX_UP:
4.391 - lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
4.392 - break;
4.393 - case LPX_DB:
4.394 - case LPX_FX:
4.395 - lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
4.396 - break;
4.397 - default: ;
4.398 - //FIXME error
4.399 - }
4.400 - } else {
4.401 - switch (b) {
4.402 - case LPX_FR:
4.403 - lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
4.404 - case LPX_LO:
4.405 - if (lo==up)
4.406 - lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
4.407 - else
4.408 - lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
4.409 - break;
4.410 - case LPX_UP:
4.411 - lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
4.412 - break;
4.413 - case LPX_DB:
4.414 - case LPX_FX:
4.415 - if (lo==up)
4.416 - lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
4.417 - else
4.418 - lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
4.419 - break;
4.420 - default: ;
4.421 - //FIXME error
4.422 - }
4.423 - }
4.424 - }
4.425 - virtual double _getRowUpperBound(int i) {
4.426 - int b=lpx_get_row_type(lp, i);
4.427 - switch (b) {
4.428 - case LPX_FR:
4.429 - case LPX_LO:
4.430 - return INF;
4.431 - case LPX_UP:
4.432 - case LPX_DB:
4.433 - case LPX_FX:
4.434 - return lpx_get_row_ub(lp, i);
4.435 - default: ;
4.436 - //FIXME error
4.437 - return 0.0;
4.438 - }
4.439 - }
4.440 - /// \e
4.441 - virtual double _getObjCoeff(int i) {
4.442 - return lpx_get_obj_coef(lp, i);
4.443 - }
4.444 - /// \e
4.445 - virtual void _setObjCoeff(int i, double obj_coef) {
4.446 - lpx_set_obj_coef(lp, i, obj_coef);
4.447 - }
4.448 - public:
4.449 - /// \e
4.450 - void solveSimplex() { lpx_simplex(lp); }
4.451 - /// \e
4.452 - void solvePrimalSimplex() { lpx_simplex(lp); }
4.453 - /// \e
4.454 - void solveDualSimplex() { lpx_simplex(lp); }
4.455 - protected:
4.456 - virtual double _getPrimal(int i) {
4.457 - return lpx_get_col_prim(lp, i);
4.458 - }
4.459 - public:
4.460 - /// \e
4.461 - double getObjVal() { return lpx_get_obj_val(lp); }
4.462 - /// \e
4.463 - int rowNum() const { return lpx_get_num_rows(lp); }
4.464 - /// \e
4.465 - int colNum() const { return lpx_get_num_cols(lp); }
4.466 - /// \e
4.467 - int warmUp() { return lpx_warm_up(lp); }
4.468 - /// \e
4.469 - void printWarmUpStatus(int i) {
4.470 - switch (i) {
4.471 - case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
4.472 - case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;
4.473 - case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
4.474 - case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
4.475 - }
4.476 - }
4.477 - /// \e
4.478 - int getPrimalStatus() { return lpx_get_prim_stat(lp); }
4.479 - /// \e
4.480 - void printPrimalStatus(int i) {
4.481 - switch (i) {
4.482 - case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
4.483 - case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;
4.484 - case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
4.485 - case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
4.486 - }
4.487 - }
4.488 - /// \e
4.489 - int getDualStatus() { return lpx_get_dual_stat(lp); }
4.490 - /// \e
4.491 - void printDualStatus(int i) {
4.492 - switch (i) {
4.493 - case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
4.494 - case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;
4.495 - case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
4.496 - case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
4.497 - }
4.498 - }
4.499 - /// Returns the status of the slack variable assigned to row \c row.
4.500 - int getRowStat(const Row& row) {
4.501 - return lpx_get_row_stat(lp, row_iter_map[row]);
4.502 - }
4.503 - /// \e
4.504 - void printRowStatus(int i) {
4.505 - switch (i) {
4.506 - case LPX_BS: cout << "LPX_BS" << endl; break;
4.507 - case LPX_NL: cout << "LPX_NL" << endl; break;
4.508 - case LPX_NU: cout << "LPX_NU" << endl; break;
4.509 - case LPX_NF: cout << "LPX_NF" << endl; break;
4.510 - case LPX_NS: cout << "LPX_NS" << endl; break;
4.511 - }
4.512 - }
4.513 - /// Returns the status of the variable assigned to column \c col.
4.514 - int getColStat(const Col& col) {
4.515 - return lpx_get_col_stat(lp, col_iter_map[col]);
4.516 - }
4.517 - /// \e
4.518 - void printColStatus(int i) {
4.519 - switch (i) {
4.520 - case LPX_BS: cout << "LPX_BS" << endl; break;
4.521 - case LPX_NL: cout << "LPX_NL" << endl; break;
4.522 - case LPX_NU: cout << "LPX_NU" << endl; break;
4.523 - case LPX_NF: cout << "LPX_NF" << endl; break;
4.524 - case LPX_NS: cout << "LPX_NS" << endl; break;
4.525 - }
4.526 - }
4.527 -
4.528 - // MIP
4.529 - /// \e
4.530 - void solveBandB() { lpx_integer(lp); }
4.531 - /// \e
4.532 - void setLP() { lpx_set_class(lp, LPX_LP); }
4.533 - /// \e
4.534 - void setMIP() { lpx_set_class(lp, LPX_MIP); }
4.535 - protected:
4.536 - /// \e
4.537 - void _setColCont(int i) { lpx_set_col_kind(lp, i, LPX_CV); }
4.538 - /// \e
4.539 - void _setColInt(int i) { lpx_set_col_kind(lp, i, LPX_IV); }
4.540 - /// \e
4.541 - double _getMIPPrimal(int i) { return lpx_mip_col_val(lp, i); }
4.542 - };
4.543 -
4.544 - /// @}
4.545 -
4.546 -} //namespace lemon
4.547 -
4.548 -#endif //LEMON_LP_SOLVER_GLPK_H
5.1 --- a/src/work/athos/lp/lp_solver_wrapper.h Tue Mar 22 16:49:30 2005 +0000
5.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
5.3 @@ -1,431 +0,0 @@
5.4 -// -*- c++ -*-
5.5 -#ifndef LEMON_LP_SOLVER_WRAPPER_H
5.6 -#define LEMON_LP_SOLVER_WRAPPER_H
5.7 -
5.8 -///\ingroup misc
5.9 -///\file
5.10 -///\brief Dijkstra algorithm.
5.11 -
5.12 -// #include <stdio.h>
5.13 -#include <stdlib.h>
5.14 -// #include <stdio>
5.15 -//#include <stdlib>
5.16 -extern "C" {
5.17 -#include "glpk.h"
5.18 -}
5.19 -
5.20 -#include <iostream>
5.21 -#include <vector>
5.22 -#include <string>
5.23 -#include <list>
5.24 -#include <memory>
5.25 -#include <utility>
5.26 -
5.27 -//#include <sage_graph.h>
5.28 -//#include <lemon/list_graph.h>
5.29 -//#include <lemon/graph_wrapper.h>
5.30 -#include <lemon/invalid.h>
5.31 -//#include <bfs_dfs.h>
5.32 -//#include <stp.h>
5.33 -//#include <lemon/max_flow.h>
5.34 -//#include <augmenting_flow.h>
5.35 -//#include <iter_map.h>
5.36 -
5.37 -using std::cout;
5.38 -using std::cin;
5.39 -using std::endl;
5.40 -
5.41 -namespace lemon {
5.42 -
5.43 -
5.44 - /// \addtogroup misc
5.45 - /// @{
5.46 -
5.47 - /// \brief A partitioned vector with iterable classes.
5.48 - ///
5.49 - /// This class implements a container in which the data is stored in an
5.50 - /// stl vector, the range is partitioned into sets and each set is
5.51 - /// doubly linked in a list.
5.52 - /// That is, each class is iterable by lemon iterators, and any member of
5.53 - /// the vector can bo moved to an other class.
5.54 - template <typename T>
5.55 - class IterablePartition {
5.56 - protected:
5.57 - struct Node {
5.58 - T data;
5.59 - int prev; //invalid az -1
5.60 - int next;
5.61 - };
5.62 - std::vector<Node> nodes;
5.63 - struct Tip {
5.64 - int first;
5.65 - int last;
5.66 - };
5.67 - std::vector<Tip> tips;
5.68 - public:
5.69 - /// The classes are indexed by integers from \c 0 to \c classNum()-1.
5.70 - int classNum() const { return tips.size(); }
5.71 - /// This lemon style iterator iterates through a class.
5.72 - class ClassIt;
5.73 - /// Constructor. The number of classes is to be given which is fixed
5.74 - /// over the life of the container.
5.75 - /// The partition classes are indexed from 0 to class_num-1.
5.76 - IterablePartition(int class_num) {
5.77 - for (int i=0; i<class_num; ++i) {
5.78 - Tip t;
5.79 - t.first=t.last=-1;
5.80 - tips.push_back(t);
5.81 - }
5.82 - }
5.83 - protected:
5.84 - void befuz(ClassIt it, int class_id) {
5.85 - if (tips[class_id].first==-1) {
5.86 - if (tips[class_id].last==-1) {
5.87 - nodes[it.i].prev=nodes[it.i].next=-1;
5.88 - tips[class_id].first=tips[class_id].last=it.i;
5.89 - }
5.90 - } else {
5.91 - nodes[it.i].prev=tips[class_id].last;
5.92 - nodes[it.i].next=-1;
5.93 - nodes[tips[class_id].last].next=it.i;
5.94 - tips[class_id].last=it.i;
5.95 - }
5.96 - }
5.97 - void kifuz(ClassIt it, int class_id) {
5.98 - if (tips[class_id].first==it.i) {
5.99 - if (tips[class_id].last==it.i) {
5.100 - tips[class_id].first=tips[class_id].last=-1;
5.101 - } else {
5.102 - tips[class_id].first=nodes[it.i].next;
5.103 - nodes[nodes[it.i].next].prev=-1;
5.104 - }
5.105 - } else {
5.106 - if (tips[class_id].last==it.i) {
5.107 - tips[class_id].last=nodes[it.i].prev;
5.108 - nodes[nodes[it.i].prev].next=-1;
5.109 - } else {
5.110 - nodes[nodes[it.i].next].prev=nodes[it.i].prev;
5.111 - nodes[nodes[it.i].prev].next=nodes[it.i].next;
5.112 - }
5.113 - }
5.114 - }
5.115 - public:
5.116 - /// A new element with data \c t is pushed into the vector and into class
5.117 - /// \c class_id.
5.118 - ClassIt push_back(const T& t, int class_id) {
5.119 - Node n;
5.120 - n.data=t;
5.121 - nodes.push_back(n);
5.122 - int i=nodes.size()-1;
5.123 - befuz(i, class_id);
5.124 - return i;
5.125 - }
5.126 - /// A member is moved to an other class.
5.127 - void set(ClassIt it, int old_class_id, int new_class_id) {
5.128 - kifuz(it.i, old_class_id);
5.129 - befuz(it.i, new_class_id);
5.130 - }
5.131 - /// Returns the data pointed by \c it.
5.132 - T& operator[](ClassIt it) { return nodes[it.i].data; }
5.133 - /// Returns the data pointed by \c it.
5.134 - const T& operator[](ClassIt it) const { return nodes[it.i].data; }
5.135 - ///.
5.136 - class ClassIt {
5.137 - friend class IterablePartition;
5.138 - protected:
5.139 - int i;
5.140 - public:
5.141 - /// Default constructor.
5.142 - ClassIt() { }
5.143 - /// This constructor constructs an iterator which points
5.144 - /// to the member of th container indexed by the integer _i.
5.145 - ClassIt(const int& _i) : i(_i) { }
5.146 - /// Invalid constructor.
5.147 - ClassIt(const Invalid&) : i(-1) { }
5.148 - };
5.149 - /// First member of class \c class_id.
5.150 - ClassIt& first(ClassIt& it, int class_id) const {
5.151 - it.i=tips[class_id].first;
5.152 - return it;
5.153 - }
5.154 - /// Next member.
5.155 - ClassIt& next(ClassIt& it) const {
5.156 - it.i=nodes[it.i].next;
5.157 - return it;
5.158 - }
5.159 - /// True iff the iterator is valid.
5.160 - bool valid(const ClassIt& it) const { return it.i!=-1; }
5.161 - };
5.162 -
5.163 - /// \brief Wrappers for LP solvers
5.164 - ///
5.165 - /// This class implements a lemon wrapper for glpk.
5.166 - /// Later other LP-solvers will be wrapped into lemon.
5.167 - /// The aim of this class is to give a general surface to different
5.168 - /// solvers, i.e. it makes possible to write algorithms using LP's,
5.169 - /// in which the solver can be changed to an other one easily.
5.170 - class LPSolverWrapper {
5.171 - public:
5.172 -
5.173 -// class Row {
5.174 -// protected:
5.175 -// int i;
5.176 -// public:
5.177 -// Row() { }
5.178 -// Row(const Invalid&) : i(0) { }
5.179 -// Row(const int& _i) : i(_i) { }
5.180 -// operator int() const { return i; }
5.181 -// };
5.182 -// class RowIt : public Row {
5.183 -// public:
5.184 -// RowIt(const Row& row) : Row(row) { }
5.185 -// };
5.186 -
5.187 -// class Col {
5.188 -// protected:
5.189 -// int i;
5.190 -// public:
5.191 -// Col() { }
5.192 -// Col(const Invalid&) : i(0) { }
5.193 -// Col(const int& _i) : i(_i) { }
5.194 -// operator int() const { return i; }
5.195 -// };
5.196 -// class ColIt : public Col {
5.197 -// ColIt(const Col& col) : Col(col) { }
5.198 -// };
5.199 -
5.200 - public:
5.201 - ///.
5.202 - LPX* lp;
5.203 - ///.
5.204 - typedef IterablePartition<int>::ClassIt RowIt;
5.205 - ///.
5.206 - IterablePartition<int> row_iter_map;
5.207 - ///.
5.208 - typedef IterablePartition<int>::ClassIt ColIt;
5.209 - ///.
5.210 - IterablePartition<int> col_iter_map;
5.211 - //std::vector<int> row_id_to_lp_row_id;
5.212 - //std::vector<int> col_id_to_lp_col_id;
5.213 - ///.
5.214 - const int VALID_ID;
5.215 - ///.
5.216 - const int INVALID_ID;
5.217 -
5.218 - public:
5.219 - ///.
5.220 - LPSolverWrapper() : lp(lpx_create_prob()),
5.221 - row_iter_map(2),
5.222 - col_iter_map(2),
5.223 - //row_id_to_lp_row_id(), col_id_to_lp_col_id(),
5.224 - VALID_ID(0), INVALID_ID(1) {
5.225 - lpx_set_int_parm(lp, LPX_K_DUAL, 1);
5.226 - }
5.227 - ///.
5.228 - ~LPSolverWrapper() {
5.229 - lpx_delete_prob(lp);
5.230 - }
5.231 - ///.
5.232 - void setMinimize() {
5.233 - lpx_set_obj_dir(lp, LPX_MIN);
5.234 - }
5.235 - ///.
5.236 - void setMaximize() {
5.237 - lpx_set_obj_dir(lp, LPX_MAX);
5.238 - }
5.239 - ///.
5.240 - ColIt addCol() {
5.241 - int i=lpx_add_cols(lp, 1);
5.242 - ColIt col_it;
5.243 - col_iter_map.first(col_it, INVALID_ID);
5.244 - if (col_iter_map.valid(col_it)) { //van hasznalhato hely
5.245 - col_iter_map.set(col_it, INVALID_ID, VALID_ID);
5.246 - col_iter_map[col_it]=i;
5.247 - //col_id_to_lp_col_id[col_iter_map[col_it]]=i;
5.248 - } else { //a cucc vegere kell inzertalni mert nincs szabad hely
5.249 - //col_id_to_lp_col_id.push_back(i);
5.250 - //int j=col_id_to_lp_col_id.size()-1;
5.251 - col_it=col_iter_map.push_back(i, VALID_ID);
5.252 - }
5.253 -// edge_index_map.set(e, i);
5.254 -// lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
5.255 -// lpx_set_obj_coef(lp, i, cost[e]);
5.256 - return col_it;
5.257 - }
5.258 - ///.
5.259 - RowIt addRow() {
5.260 - int i=lpx_add_rows(lp, 1);
5.261 - RowIt row_it;
5.262 - row_iter_map.first(row_it, INVALID_ID);
5.263 - if (row_iter_map.valid(row_it)) { //van hasznalhato hely
5.264 - row_iter_map.set(row_it, INVALID_ID, VALID_ID);
5.265 - row_iter_map[row_it]=i;
5.266 - } else { //a cucc vegere kell inzertalni mert nincs szabad hely
5.267 - row_it=row_iter_map.push_back(i, VALID_ID);
5.268 - }
5.269 - return row_it;
5.270 - }
5.271 - //pair<RowIt, double>-bol kell megadni egy std range-et
5.272 - ///.
5.273 - template <typename Begin, typename End>
5.274 - void setColCoeffs(const ColIt& col_it,
5.275 - Begin begin, End end) {
5.276 - int mem_length=1+lpx_get_num_rows(lp);
5.277 - int* indices = new int[mem_length];
5.278 - double* doubles = new double[mem_length];
5.279 - int length=0;
5.280 - for ( ; begin!=end; ++begin) {
5.281 - ++length;
5.282 - indices[length]=row_iter_map[begin->first];
5.283 - doubles[length]=begin->second;
5.284 - }
5.285 - lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
5.286 - delete [] indices;
5.287 - delete [] doubles;
5.288 - }
5.289 - //pair<ColIt, double>-bol kell megadni egy std range-et
5.290 - ///.
5.291 - template <typename Begin, typename End>
5.292 - void setRowCoeffs(const RowIt& row_it,
5.293 - Begin begin, End end) {
5.294 - int mem_length=1+lpx_get_num_cols(lp);
5.295 - int* indices = new int[mem_length];
5.296 - double* doubles = new double[mem_length];
5.297 - int length=0;
5.298 - for ( ; begin!=end; ++begin) {
5.299 - ++length;
5.300 - indices[length]=col_iter_map[begin->first];
5.301 - doubles[length]=begin->second;
5.302 - }
5.303 - lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
5.304 - delete [] indices;
5.305 - delete [] doubles;
5.306 - }
5.307 - ///.
5.308 - void eraseCol(const ColIt& col_it) {
5.309 - col_iter_map.set(col_it, VALID_ID, INVALID_ID);
5.310 - int cols[2];
5.311 - cols[1]=col_iter_map[col_it];
5.312 - lpx_del_cols(lp, 1, cols);
5.313 - col_iter_map[col_it]=0; //glpk specifikus
5.314 - ColIt it;
5.315 - for (col_iter_map.first(it, VALID_ID);
5.316 - col_iter_map.valid(it); col_iter_map.next(it)) {
5.317 - if (col_iter_map[it]>cols[1]) --col_iter_map[it];
5.318 - }
5.319 - }
5.320 - ///.
5.321 - void eraseRow(const RowIt& row_it) {
5.322 - row_iter_map.set(row_it, VALID_ID, INVALID_ID);
5.323 - int rows[2];
5.324 - rows[1]=row_iter_map[row_it];
5.325 - lpx_del_rows(lp, 1, rows);
5.326 - row_iter_map[row_it]=0; //glpk specifikus
5.327 - RowIt it;
5.328 - for (row_iter_map.first(it, VALID_ID);
5.329 - row_iter_map.valid(it); row_iter_map.next(it)) {
5.330 - if (row_iter_map[it]>rows[1]) --row_iter_map[it];
5.331 - }
5.332 - }
5.333 - ///.
5.334 - void setColBounds(const ColIt& col_it, int bound_type,
5.335 - double lo, double up) {
5.336 - lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
5.337 - }
5.338 - ///.
5.339 - double getObjCoef(const ColIt& col_it) {
5.340 - return lpx_get_obj_coef(lp, col_iter_map[col_it]);
5.341 - }
5.342 - ///.
5.343 - void setRowBounds(const RowIt& row_it, int bound_type,
5.344 - double lo, double up) {
5.345 - lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
5.346 - }
5.347 - ///.
5.348 - void setObjCoef(const ColIt& col_it, double obj_coef) {
5.349 - lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
5.350 - }
5.351 - ///.
5.352 - void solveSimplex() { lpx_simplex(lp); }
5.353 - ///.
5.354 - void solvePrimalSimplex() { lpx_simplex(lp); }
5.355 - ///.
5.356 - void solveDualSimplex() { lpx_simplex(lp); }
5.357 - ///.
5.358 - double getPrimal(const ColIt& col_it) {
5.359 - return lpx_get_col_prim(lp, col_iter_map[col_it]);
5.360 - }
5.361 - ///.
5.362 - double getObjVal() { return lpx_get_obj_val(lp); }
5.363 - ///.
5.364 - int rowNum() const { return lpx_get_num_rows(lp); }
5.365 - ///.
5.366 - int colNum() const { return lpx_get_num_cols(lp); }
5.367 - ///.
5.368 - int warmUp() { return lpx_warm_up(lp); }
5.369 - ///.
5.370 - void printWarmUpStatus(int i) {
5.371 - switch (i) {
5.372 - case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
5.373 - case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;
5.374 - case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
5.375 - case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
5.376 - }
5.377 - }
5.378 - ///.
5.379 - int getPrimalStatus() { return lpx_get_prim_stat(lp); }
5.380 - ///.
5.381 - void printPrimalStatus(int i) {
5.382 - switch (i) {
5.383 - case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
5.384 - case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;
5.385 - case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
5.386 - case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
5.387 - }
5.388 - }
5.389 - ///.
5.390 - int getDualStatus() { return lpx_get_dual_stat(lp); }
5.391 - ///.
5.392 - void printDualStatus(int i) {
5.393 - switch (i) {
5.394 - case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
5.395 - case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;
5.396 - case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
5.397 - case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
5.398 - }
5.399 - }
5.400 - /// Returns the status of the slack variable assigned to row \c row_it.
5.401 - int getRowStat(const RowIt& row_it) {
5.402 - return lpx_get_row_stat(lp, row_iter_map[row_it]);
5.403 - }
5.404 - ///.
5.405 - void printRowStatus(int i) {
5.406 - switch (i) {
5.407 - case LPX_BS: cout << "LPX_BS" << endl; break;
5.408 - case LPX_NL: cout << "LPX_NL" << endl; break;
5.409 - case LPX_NU: cout << "LPX_NU" << endl; break;
5.410 - case LPX_NF: cout << "LPX_NF" << endl; break;
5.411 - case LPX_NS: cout << "LPX_NS" << endl; break;
5.412 - }
5.413 - }
5.414 - /// Returns the status of the variable assigned to column \c col_it.
5.415 - int getColStat(const ColIt& col_it) {
5.416 - return lpx_get_col_stat(lp, col_iter_map[col_it]);
5.417 - }
5.418 - ///.
5.419 - void printColStatus(int i) {
5.420 - switch (i) {
5.421 - case LPX_BS: cout << "LPX_BS" << endl; break;
5.422 - case LPX_NL: cout << "LPX_NL" << endl; break;
5.423 - case LPX_NU: cout << "LPX_NU" << endl; break;
5.424 - case LPX_NF: cout << "LPX_NF" << endl; break;
5.425 - case LPX_NS: cout << "LPX_NS" << endl; break;
5.426 - }
5.427 - }
5.428 - };
5.429 -
5.430 - /// @}
5.431 -
5.432 -} //namespace lemon
5.433 -
5.434 -#endif //LEMON_LP_SOLVER_WRAPPER_H
6.1 --- a/src/work/athos/lp/magic_square.cc Tue Mar 22 16:49:30 2005 +0000
6.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
6.3 @@ -1,99 +0,0 @@
6.4 -// -*- c++ -*-
6.5 -#include <iostream>
6.6 -#include <fstream>
6.7 -
6.8 -#include <lemon/time_measure.h>
6.9 -#include <lp_solver_glpk.h>
6.10 -
6.11 -using std::cout;
6.12 -using std::endl;
6.13 -using namespace lemon;
6.14 -
6.15 -/*
6.16 - On an 1537Mhz PC, the run times with
6.17 - glpk are the following.
6.18 - for n=3,4, some secondes
6.19 - for n=5, 25 hours
6.20 - */
6.21 -
6.22 -int main(int, char **) {
6.23 - const int n=3;
6.24 - const double row_sum=(1.0+n*n)*n/2;
6.25 - Timer ts;
6.26 - ts.reset();
6.27 - typedef LpGlpk LPSolver;
6.28 - typedef LPSolver::Col Col;
6.29 - LPSolver lp;
6.30 - typedef std::map<std::pair<int, int>, Col> Coords;
6.31 - Coords x;
6.32 - // we create a new variable for each entry
6.33 - // of the magic square
6.34 - for (int i=1; i<=n; ++i) {
6.35 - for (int j=1; j<=n; ++j) {
6.36 - Col col=lp.addCol();
6.37 - x[std::make_pair(i,j)]=col;
6.38 - lp.setColLowerBound(col, 1.0);
6.39 - lp.setColUpperBound(col, double(n*n));
6.40 - }
6.41 - }
6.42 - LPSolver::Expression expr3, expr4;
6.43 - for (int i=1; i<=n; ++i) {
6.44 - LPSolver::Expression expr1, expr2;
6.45 - for (int j=1; j<=n; ++j) {
6.46 - expr1+=x[std::make_pair(i, j)];
6.47 - expr2+=x[std::make_pair(j, i)];
6.48 - }
6.49 -
6.50 - // sum of rows and columns
6.51 - lp.addRow(expr1==row_sum);
6.52 - lp.addRow(expr2==row_sum);
6.53 - cout <<"Expr1: "<<expr1<<endl;
6.54 - cout <<"Expr2: "<<expr2<<endl;
6.55 -
6.56 - expr3+=x[std::make_pair(i, i)];
6.57 - expr4+=x[std::make_pair(i, (n+1)-i)];
6.58 - }
6.59 - cout <<"Expr3: "<<expr3<<endl;
6.60 - cout <<"Expr4: "<<expr4<<endl;
6.61 - // sum of the diagonal entries
6.62 - lp.addRow(expr3==row_sum);
6.63 - lp.addRow(expr4==row_sum);
6.64 - lp.solveSimplex();
6.65 - cout << "elapsed time: " << ts << endl;
6.66 - for (int i=1; i<=n; ++i) {
6.67 - for (int j=1; j<=n; ++j) {
6.68 - cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)])
6.69 - << endl;
6.70 - }
6.71 - }
6.72 -
6.73 -
6.74 -
6.75 -// // we make new binary variables for each pair of
6.76 -// // entries of the square to achieve that
6.77 -// // the values of different entries are different
6.78 -// lp.setMIP();
6.79 -// for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) {
6.80 -// Coords::const_iterator jt=it; ++jt;
6.81 -// for(; jt!=x.end(); ++jt) {
6.82 -// Col col1=(*it).second;
6.83 -// Col col2=(*jt).second;
6.84 -// Col col=lp.addCol();
6.85 -// lp.setColLowerBound(col, 0.0);
6.86 -// lp.setColUpperBound(col, 1.0);
6.87 -// lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0);
6.88 -// lp.setColInt(col);
6.89 -// }
6.90 -// }
6.91 -// cout << "elapsed time: " << ts << endl;
6.92 -// lp.solveSimplex();
6.93 -// // let's solve the integer problem
6.94 -// lp.solveBandB();
6.95 -// cout << "elapsed time: " << ts << endl;
6.96 -// for (int i=1; i<=n; ++i) {
6.97 -// for (int j=1; j<=n; ++j) {
6.98 -// cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)])
6.99 -// << endl;
6.100 -// }
6.101 -// }
6.102 -}
7.1 --- a/src/work/athos/lp/makefile Tue Mar 22 16:49:30 2005 +0000
7.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
7.3 @@ -1,73 +0,0 @@
7.4 -#INCLUDEDIRS ?= -I.. -I. -I./{marci,jacint,alpar,klao,akos}
7.5 -#GLPKROOT = /usr/local/glpk-4.4
7.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
7.7 -#INCLUDEDIRS ?= -I../.. -I../.. -I../../{marci,jacint,alpar,klao,akos} -I/usr/local/glpk-4.4/include
7.8 -CXXFLAGS = -g -O2 -W -Wall $(INCLUDEDIRS) -ansi -pedantic
7.9 -LDFLAGS = -lglpk#-lcplex -lm -lpthread -lilocplex -L/usr/local/cplex/cplex75/lib/i86_linux2_glibc2.2_gcc3.0/static_mt# -L$(GLPKROOT)/lib
7.10 -
7.11 -BINARIES = magic_square max_flow_expression #expression_test max_flow_by_lp# sample sample2 sample11 sample15
7.12 -
7.13 -#include ../makefile
7.14 -
7.15 -# Hat, ez elismerem, hogy nagyon ronda, de mukodik minden altalam
7.16 -# ismert rendszeren :-) (Misi)
7.17 -ifdef GCCVER
7.18 -CXX := g++-$(GCCVER)
7.19 -else
7.20 -CXX := $(shell type -p g++-3.3 || type -p g++-3.2 || type -p g++-3.0 || type -p g++-3 || echo g++)
7.21 -endif
7.22 -
7.23 -ifdef DEBUG
7.24 -CXXFLAGS += -DDEBUG
7.25 -endif
7.26 -
7.27 -CC := $(CXX)
7.28 -
7.29 -
7.30 -all: $(BINARIES)
7.31 -
7.32 -################
7.33 -# Minden binarishoz egy sor, felsorolva, hogy mely object file-okbol
7.34 -# all elo.
7.35 -# Kiveve ha siman file.cc -> file esetrol van szo, amikor is nem kell
7.36 -# irni semmit.
7.37 -
7.38 -#proba: proba.o seged.o
7.39 -
7.40 -################
7.41 -
7.42 -
7.43 -# .depend dep depend:
7.44 -# -$(CXX) $(CXXFLAGS) -M $(BINARIES:=.cc) > .depend
7.45 -
7.46 -#makefile: .depend
7.47 -#sinclude .depend
7.48 -#moert nem megy az eredeti /usr/bin/ld-vel?
7.49 -
7.50 -# %: %.o
7.51 -# $(CXX) -o $@ $< $(LDFLAGS)
7.52 -
7.53 -# %.o: %.cc
7.54 -# $(CXX) $(CXXFLAGS) -c $<
7.55 -
7.56 -%: %.cc
7.57 - $(CXX) $(CXXFLAGS) -o $@ $< $(LDFLAGS)
7.58 -
7.59 -sample11prof: sample11prof.o
7.60 - $(CXX) -pg -o sample11prof sample11prof.o -L$(GLPKROOT)/lib -lglpk
7.61 -sample11prof.o: sample11.cc
7.62 - $(CXX) -pg $(CXXFLAGS) -c -o sample11prof.o sample11.cc
7.63 -
7.64 -# sample.o: sample.cc
7.65 -# $(CXX) $(CXXFLAGS) -c -o sample.o sample.cc
7.66 -
7.67 -# sample2: sample2.o
7.68 -# $(CXX) -o sample2 sample2.o -L/usr/local/glpk-4.4/lib -lglpk
7.69 -# sample2.o: sample2.cc
7.70 -# $(CXX) $(CXXFLAGS) -c -o sample2.o sample2.cc
7.71 -
7.72 -
7.73 -clean:
7.74 - $(RM) *.o $(BINARIES) .depend
7.75 -
7.76 -.PHONY: all clean dep depend
8.1 --- a/src/work/athos/lp/max_flow_by_lp.cc Tue Mar 22 16:49:30 2005 +0000
8.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
8.3 @@ -1,186 +0,0 @@
8.4 -// -*- c++ -*-
8.5 -#include <iostream>
8.6 -#include <fstream>
8.7 -
8.8 -#include <lemon/smart_graph.h>
8.9 -#include <lemon/list_graph.h>
8.10 -#include <lemon/dimacs.h>
8.11 -#include <lemon/time_measure.h>
8.12 -//#include <graph_wrapper.h>
8.13 -#include <lemon/preflow.h>
8.14 -#include <augmenting_flow.h>
8.15 -//#include <preflow_res.h>
8.16 -//#include <lp_solver_wrapper_2.h>
8.17 -#include <min_cost_gen_flow.h>
8.18 -
8.19 -// Use a DIMACS max flow file as stdin.
8.20 -// max_flow_demo < dimacs_max_flow_file
8.21 -
8.22 -using namespace lemon;
8.23 -
8.24 -int main(int, char **) {
8.25 -
8.26 - typedef ListGraph MutableGraph;
8.27 - typedef ListGraph Graph;
8.28 - typedef Graph::Node Node;
8.29 - typedef Graph::Edge Edge;
8.30 - typedef Graph::EdgeIt EdgeIt;
8.31 -
8.32 - Graph g;
8.33 -
8.34 - Node s, t;
8.35 - Graph::EdgeMap<int> cap(g);
8.36 - //readDimacsMaxFlow(std::cin, g, s, t, cap);
8.37 - readDimacs(std::cin, g, cap, s, t);
8.38 - Timer ts;
8.39 - Graph::EdgeMap<int> flow(g); //0 flow
8.40 - Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
8.41 - max_flow_test(g, s, t, cap, flow);
8.42 - AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
8.43 - augmenting_flow_test(g, s, t, cap, flow);
8.44 -
8.45 - Graph::NodeMap<bool> cut(g);
8.46 -
8.47 - {
8.48 - std::cout << "preflow ..." << std::endl;
8.49 - ts.reset();
8.50 - max_flow_test.run();
8.51 - std::cout << "elapsed time: " << ts << std::endl;
8.52 - std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
8.53 - max_flow_test.minCut(cut);
8.54 -
8.55 - for (EdgeIt e(g); e!=INVALID; ++e) {
8.56 - if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
8.57 - std::cout << "Slackness does not hold!" << std::endl;
8.58 - if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
8.59 - std::cout << "Slackness does not hold!" << std::endl;
8.60 - }
8.61 - }
8.62 -
8.63 -// {
8.64 -// std::cout << "preflow ..." << std::endl;
8.65 -// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
8.66 -// ts.reset();
8.67 -// max_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
8.68 -// std::cout << "elapsed time: " << ts << std::endl;
8.69 -// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
8.70 -
8.71 -// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
8.72 -// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
8.73 -// std::cout << "Slackness does not hold!" << std::endl;
8.74 -// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
8.75 -// std::cout << "Slackness does not hold!" << std::endl;
8.76 -// }
8.77 -// }
8.78 -
8.79 -// {
8.80 -// std::cout << "wrapped preflow ..." << std::endl;
8.81 -// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
8.82 -// ts.reset();
8.83 -// pre_flow_res.run();
8.84 -// std::cout << "elapsed time: " << ts << std::endl;
8.85 -// std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
8.86 -// }
8.87 -
8.88 - {
8.89 - std::cout << "physical blocking flow augmentation ..." << std::endl;
8.90 - for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
8.91 - ts.reset();
8.92 - int i=0;
8.93 - while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
8.94 - std::cout << "elapsed time: " << ts << std::endl;
8.95 - std::cout << "number of augmentation phases: " << i << std::endl;
8.96 - std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
8.97 -
8.98 - for (EdgeIt e(g); e!=INVALID; ++e) {
8.99 - if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
8.100 - std::cout << "Slackness does not hold!" << std::endl;
8.101 - if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
8.102 - std::cout << "Slackness does not hold!" << std::endl;
8.103 - }
8.104 - }
8.105 -
8.106 -// {
8.107 -// std::cout << "faster physical blocking flow augmentation ..." << std::endl;
8.108 -// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
8.109 -// ts.reset();
8.110 -// int i=0;
8.111 -// while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
8.112 -// std::cout << "elapsed time: " << ts << std::endl;
8.113 -// std::cout << "number of augmentation phases: " << i << std::endl;
8.114 -// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
8.115 -// }
8.116 -
8.117 - {
8.118 - std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
8.119 - for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
8.120 - ts.reset();
8.121 - int i=0;
8.122 - while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
8.123 - std::cout << "elapsed time: " << ts << std::endl;
8.124 - std::cout << "number of augmentation phases: " << i << std::endl;
8.125 - std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
8.126 -
8.127 - for (EdgeIt e(g); e!=INVALID; ++e) {
8.128 - if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
8.129 - std::cout << "Slackness does not hold!" << std::endl;
8.130 - if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
8.131 - std::cout << "Slackness does not hold!" << std::endl;
8.132 - }
8.133 - }
8.134 -
8.135 -// {
8.136 -// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
8.137 -// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
8.138 -// ts.reset();
8.139 -// int i=0;
8.140 -// while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
8.141 -// std::cout << "elapsed time: " << ts << std::endl;
8.142 -// std::cout << "number of augmentation phases: " << i << std::endl;
8.143 -// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
8.144 -
8.145 -// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
8.146 -// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
8.147 -// std::cout << "Slackness does not hold!" << std::endl;
8.148 -// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
8.149 -// std::cout << "Slackness does not hold!" << std::endl;
8.150 -// }
8.151 -// }
8.152 -
8.153 -// {
8.154 -// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
8.155 -// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
8.156 -// ts.reset();
8.157 -// int i=0;
8.158 -// while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
8.159 -// std::cout << "elapsed time: " << ts << std::endl;
8.160 -// std::cout << "number of augmentation phases: " << i << std::endl;
8.161 -// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
8.162 -
8.163 -// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
8.164 -// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
8.165 -// std::cout << "Slackness does not hold!" << std::endl;
8.166 -// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
8.167 -// std::cout << "Slackness does not hold!" << std::endl;
8.168 -// }
8.169 -// }
8.170 -
8.171 - ts.reset();
8.172 -
8.173 - Edge e=g.addEdge(t, s);
8.174 - Graph::EdgeMap<int> cost(g, 0);
8.175 - cost.set(e, -1);
8.176 - cap.set(e, 10000);
8.177 - typedef ConstMap<Node, int> Excess;
8.178 - Excess excess(0);
8.179 - typedef ConstMap<Edge, int> LCap;
8.180 - LCap lcap(0);
8.181 -
8.182 - MinCostGenFlow<Graph, int, Excess, LCap>
8.183 - min_cost(g, excess, lcap, cap, flow, cost);
8.184 - min_cost.feasible();
8.185 - min_cost.runByLP();
8.186 -
8.187 - std::cout << "elapsed time: " << ts << std::endl;
8.188 - std::cout << "flow value: "<< flow[e] << std::endl;
8.189 -}
9.1 --- a/src/work/athos/lp/max_flow_expression.cc Tue Mar 22 16:49:30 2005 +0000
9.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
9.3 @@ -1,109 +0,0 @@
9.4 -// -*- c++ -*-
9.5 -#include <iostream>
9.6 -#include <fstream>
9.7 -
9.8 -#include <lemon/graph_utils.h>
9.9 -#include <lemon/smart_graph.h>
9.10 -#include <lemon/list_graph.h>
9.11 -#include <lemon/dimacs.h>
9.12 -#include <lemon/time_measure.h>
9.13 -#include <lp_solver_glpk.h>
9.14 -
9.15 -using std::cout;
9.16 -using std::endl;
9.17 -using namespace lemon;
9.18 -
9.19 -template<typename Edge, typename EdgeIndexMap>
9.20 -class PrimalMap {
9.21 -protected:
9.22 - LpGlpk* lp;
9.23 - EdgeIndexMap* edge_index_map;
9.24 -public:
9.25 - PrimalMap(LpGlpk& _lp, EdgeIndexMap& _edge_index_map) :
9.26 - lp(&_lp), edge_index_map(&_edge_index_map) { }
9.27 - double operator[](Edge e) const {
9.28 - return lp->getPrimal((*edge_index_map)[e]);
9.29 - }
9.30 -};
9.31 -
9.32 -// Use a DIMACS max flow file as stdin.
9.33 -// max_flow_expresion < dimacs_max_flow_file
9.34 -
9.35 -int main(int, char **) {
9.36 -
9.37 - typedef ListGraph Graph;
9.38 - typedef Graph::Node Node;
9.39 - typedef Graph::Edge Edge;
9.40 - typedef Graph::EdgeIt EdgeIt;
9.41 -
9.42 - Graph g;
9.43 -
9.44 - Node s, t;
9.45 - Graph::EdgeMap<int> cap(g);
9.46 - readDimacs(std::cin, g, cap, s, t);
9.47 - Timer ts;
9.48 -
9.49 - typedef LpGlpk LPSolver;
9.50 - LPSolver lp;
9.51 - lp.setMaximize();
9.52 - typedef LPSolver::Col Col;
9.53 - typedef LPSolver::Row Row;
9.54 - typedef Graph::EdgeMap<Col> EdgeIndexMap;
9.55 - typedef Graph::NodeMap<Row> NodeIndexMap;
9.56 - EdgeIndexMap edge_index_map(g);
9.57 - NodeIndexMap node_index_map(g);
9.58 - PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
9.59 -
9.60 - // nonnegativity of flow and capacity function
9.61 - for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
9.62 - Col col=lp.addCol();
9.63 - edge_index_map.set(e, col);
9.64 - // interesting property in GLPK:
9.65 - // if you change the order of the following two lines, the
9.66 - // two runs of GLPK are extremely different
9.67 - lp.setColLowerBound(col, 0);
9.68 - lp.setColUpperBound(col, cap[e]);
9.69 - }
9.70 -
9.71 - for (Graph::NodeIt n(g); n!=INVALID; ++n) {
9.72 - LPSolver::Expression expr;
9.73 - for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
9.74 - expr+=edge_index_map[e];
9.75 - for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
9.76 - expr-=edge_index_map[e];
9.77 - // cost function
9.78 - if (n==s) {
9.79 - lp.setObjCoeffs(expr);
9.80 - }
9.81 - // flow conservation constraints
9.82 - if ((n!=s) && (n!=t)) {
9.83 - node_index_map.set(n, lp.addRow(expr == 0.0));
9.84 - }
9.85 - }
9.86 - lp.solveSimplex();
9.87 - cout << "elapsed time: " << ts << endl;
9.88 -// cout << "rows:" << endl;
9.89 -// for (
9.90 -// LPSolver::Rows::ClassIt i(lp.row_iter_map, 0);
9.91 -// i!=INVALID;
9.92 -// ++i) {
9.93 -// cout << i << " ";
9.94 -// }
9.95 -// cout << endl;
9.96 -// cout << "cols:" << endl;
9.97 -// for (
9.98 -// LPSolver::Cols::ClassIt i(lp.col_iter_map, 0);
9.99 -// i!=INVALID;
9.100 -// ++i) {
9.101 -// cout << i << " ";
9.102 -// }
9.103 -// cout << endl;
9.104 - lp.setMIP();
9.105 - cout << "elapsed time: " << ts << endl;
9.106 - for (LPSolver::Cols::ClassIt it(lp.col_iter_map ,1); it!=INVALID; ++it) {
9.107 - lp.setColInt(it);
9.108 - }
9.109 - cout << "elapsed time: " << ts << endl;
9.110 - lp.solveBandB();
9.111 - cout << "elapsed time: " << ts << endl;
9.112 -}
10.1 --- a/src/work/athos/lp/min_cost_gen_flow.h Tue Mar 22 16:49:30 2005 +0000
10.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
10.3 @@ -1,268 +0,0 @@
10.4 -// -*- c++ -*-
10.5 -#ifndef LEMON_MIN_COST_GEN_FLOW_H
10.6 -#define LEMON_MIN_COST_GEN_FLOW_H
10.7 -#include <iostream>
10.8 -//#include <fstream>
10.9 -
10.10 -#include <lemon/smart_graph.h>
10.11 -#include <lemon/list_graph.h>
10.12 -//#include <lemon/dimacs.h>
10.13 -//#include <lemon/time_measure.h>
10.14 -//#include <graph_wrapper.h>
10.15 -#include <lemon/preflow.h>
10.16 -#include <lemon/min_cost_flow.h>
10.17 -//#include <augmenting_flow.h>
10.18 -//#include <preflow_res.h>
10.19 -#include <work/marci/merge_node_graph_wrapper.h>
10.20 -#include <work/marci/lp/lp_solver_wrapper_3.h>
10.21 -
10.22 -namespace lemon {
10.23 -
10.24 - template<typename Edge, typename EdgeIndexMap>
10.25 - class PrimalMap {
10.26 - protected:
10.27 - LPGLPK* lp;
10.28 - EdgeIndexMap* edge_index_map;
10.29 - public:
10.30 - PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) :
10.31 - lp(&_lp), edge_index_map(&_edge_index_map) { }
10.32 - double operator[](Edge e) const {
10.33 - return lp->getPrimal((*edge_index_map)[e]);
10.34 - }
10.35 - };
10.36 -
10.37 - // excess: rho-delta egyelore csak =0-ra.
10.38 - template <typename Graph, typename Num,
10.39 - typename Excess=typename Graph::template NodeMap<Num>,
10.40 - typename LCapMap=typename Graph::template EdgeMap<Num>,
10.41 - typename CapMap=typename Graph::template EdgeMap<Num>,
10.42 - typename FlowMap=typename Graph::template EdgeMap<Num>,
10.43 - typename CostMap=typename Graph::template EdgeMap<Num> >
10.44 - class MinCostGenFlow {
10.45 - protected:
10.46 - const Graph& g;
10.47 - const Excess& excess;
10.48 - const LCapMap& lcapacity;
10.49 - const CapMap& capacity;
10.50 - FlowMap& flow;
10.51 - const CostMap& cost;
10.52 - public:
10.53 - MinCostGenFlow(const Graph& _g, const Excess& _excess,
10.54 - const LCapMap& _lcapacity, const CapMap& _capacity,
10.55 - FlowMap& _flow,
10.56 - const CostMap& _cost) :
10.57 - g(_g), excess(_excess), lcapacity(_lcapacity),
10.58 - capacity(_capacity), flow(_flow), cost(_cost) { }
10.59 - bool feasible() {
10.60 - // std::cout << "making new vertices..." << std::endl;
10.61 - typedef ListGraph Graph2;
10.62 - Graph2 g2;
10.63 - typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
10.64 - // std::cout << "merging..." << std::endl;
10.65 - GW gw(g, g2);
10.66 - typename GW::Node s(INVALID, g2.addNode(), true);
10.67 - typename GW::Node t(INVALID, g2.addNode(), true);
10.68 - typedef SmartGraph Graph3;
10.69 - // std::cout << "making extender graph..." << std::endl;
10.70 - typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
10.71 -// {
10.72 -// checkConcept<StaticGraph, GWW>();
10.73 -// }
10.74 - GWW gww(gw);
10.75 - typedef AugmentingGraphWrapper<GW, GWW> GWWW;
10.76 - GWWW gwww(gw, gww);
10.77 -
10.78 - // std::cout << "making new edges..." << std::endl;
10.79 - typename GWWW::template EdgeMap<Num> translated_cap(gwww);
10.80 -
10.81 - for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
10.82 - translated_cap.set(typename GWWW::Edge(e,INVALID,false),
10.83 - capacity[e]-lcapacity[e]);
10.84 - // cout << "t_cap " << gw.id(e) << " "
10.85 - // << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
10.86 - }
10.87 -
10.88 - Num expected=0;
10.89 -
10.90 - // std::cout << "making new edges 2..." << std::endl;
10.91 - for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
10.92 - Num a=0;
10.93 - for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
10.94 - a+=lcapacity[e];
10.95 - for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
10.96 - a-=lcapacity[e];
10.97 - if (excess[n]>a) {
10.98 - typename GWW::Edge e=
10.99 - gww.addEdge(typename GW::Node(n,INVALID,false), t);
10.100 - translated_cap.set(typename GWWW::Edge(INVALID, e, true),
10.101 - excess[n]-a);
10.102 - // std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
10.103 - }
10.104 - if (excess[n]<a) {
10.105 - typename GWW::Edge e=
10.106 - gww.addEdge(s, typename GW::Node(n,INVALID,false));
10.107 - translated_cap.set(typename GWWW::Edge(INVALID, e, true),
10.108 - a-excess[n]);
10.109 - expected+=a-excess[n];
10.110 - // std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
10.111 - }
10.112 - }
10.113 -
10.114 - // std::cout << "preflow..." << std::endl;
10.115 - typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
10.116 - Preflow<GWWW, Num> preflow(gwww, s, t,
10.117 - translated_cap, translated_flow);
10.118 - preflow.run();
10.119 - // std::cout << "fv: " << preflow.flowValue() << std::endl;
10.120 - // std::cout << "expected: " << expected << std::endl;
10.121 -
10.122 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
10.123 - typename GW::Edge ew(e, INVALID, false);
10.124 - typename GWWW::Edge ewww(ew, INVALID, false);
10.125 - flow.set(e, translated_flow[ewww]+lcapacity[e]);
10.126 - }
10.127 - return (preflow.flowValue()>=expected);
10.128 - }
10.129 - // for nonnegative costs
10.130 - bool run() {
10.131 - // std::cout << "making new vertices..." << std::endl;
10.132 - typedef ListGraph Graph2;
10.133 - Graph2 g2;
10.134 - typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
10.135 - // std::cout << "merging..." << std::endl;
10.136 - GW gw(g, g2);
10.137 - typename GW::Node s(INVALID, g2.addNode(), true);
10.138 - typename GW::Node t(INVALID, g2.addNode(), true);
10.139 - typedef SmartGraph Graph3;
10.140 - // std::cout << "making extender graph..." << std::endl;
10.141 - typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
10.142 -// {
10.143 -// checkConcept<StaticGraph, GWW>();
10.144 -// }
10.145 - GWW gww(gw);
10.146 - typedef AugmentingGraphWrapper<GW, GWW> GWWW;
10.147 - GWWW gwww(gw, gww);
10.148 -
10.149 - // std::cout << "making new edges..." << std::endl;
10.150 - typename GWWW::template EdgeMap<Num> translated_cap(gwww);
10.151 -
10.152 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
10.153 - typename GW::Edge ew(e, INVALID, false);
10.154 - typename GWWW::Edge ewww(ew, INVALID, false);
10.155 - translated_cap.set(ewww, capacity[e]-lcapacity[e]);
10.156 - // cout << "t_cap " << g.id(e) << " "
10.157 - // << translated_cap[ewww] << endl;
10.158 - }
10.159 -
10.160 - Num expected=0;
10.161 -
10.162 - // std::cout << "making new edges 2..." << std::endl;
10.163 - for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
10.164 - // std::cout << "node: " << g.id(n) << std::endl;
10.165 - Num a=0;
10.166 - for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) {
10.167 - a+=lcapacity[e];
10.168 - // std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl;
10.169 - }
10.170 - for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) {
10.171 - a-=lcapacity[e];
10.172 - // std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl;
10.173 - }
10.174 - // std::cout << "excess " << g.id(n) << ": " << a << std::endl;
10.175 - if (0>a) {
10.176 - typename GWW::Edge e=
10.177 - gww.addEdge(typename GW::Node(n,INVALID,false), t);
10.178 - translated_cap.set(typename GWWW::Edge(INVALID, e, true),
10.179 - -a);
10.180 - // std::cout << g.id(n) << "->t " << -a << std::endl;
10.181 - }
10.182 - if (0<a) {
10.183 - typename GWW::Edge e=
10.184 - gww.addEdge(s, typename GW::Node(n,INVALID,false));
10.185 - translated_cap.set(typename GWWW::Edge(INVALID, e, true),
10.186 - a);
10.187 - expected+=a;
10.188 - // std::cout << "s->" << g.id(n) << " "<< a <<std:: endl;
10.189 - }
10.190 - }
10.191 -
10.192 - // std::cout << "preflow..." << std::endl;
10.193 - typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0);
10.194 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
10.195 - translated_cost.set(typename GWWW::Edge(
10.196 - typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]);
10.197 - }
10.198 - // typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
10.199 - MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>,
10.200 - typename GWWW::template EdgeMap<Num> >
10.201 - min_cost_flow(gwww, translated_cost, translated_cap,
10.202 - s, t);
10.203 - while (min_cost_flow.augment()) { }
10.204 - std::cout << "fv: " << min_cost_flow.flowValue() << std::endl;
10.205 - std::cout << "expected: " << expected << std::endl;
10.206 -
10.207 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
10.208 - typename GW::Edge ew(e, INVALID, false);
10.209 - typename GWWW::Edge ewww(ew, INVALID, false);
10.210 - // std::cout << g.id(e) << " " << flow[e] << std::endl;
10.211 - flow.set(e, lcapacity[e]+
10.212 - min_cost_flow.getFlow()[ewww]);
10.213 - }
10.214 - return (min_cost_flow.flowValue()>=expected);
10.215 - }
10.216 - void runByLP() {
10.217 - typedef LPGLPK LPSolver;
10.218 - LPSolver lp;
10.219 - lp.setMinimize();
10.220 - typedef LPSolver::ColIt ColIt;
10.221 - typedef LPSolver::RowIt RowIt;
10.222 - typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
10.223 - EdgeIndexMap edge_index_map(g);
10.224 - PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
10.225 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
10.226 - ColIt col_it=lp.addCol();
10.227 - edge_index_map.set(e, col_it);
10.228 - if (lcapacity[e]==capacity[e])
10.229 - lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]);
10.230 - else
10.231 - lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]);
10.232 - lp.setObjCoef(col_it, cost[e]);
10.233 - }
10.234 - LPSolver::ColIt col_it;
10.235 - for (lp.col_iter_map.first(col_it, lp.VALID_CLASS);
10.236 - lp.col_iter_map.valid(col_it);
10.237 - lp.col_iter_map.next(col_it)) {
10.238 -// std::cout << "ize " << lp.col_iter_map[col_it] << std::endl;
10.239 - }
10.240 - for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
10.241 - typename Graph::template EdgeMap<Num> coeffs(g, 0);
10.242 - for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
10.243 - coeffs.set(e, coeffs[e]+1);
10.244 - for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
10.245 - coeffs.set(e, coeffs[e]-1);
10.246 - RowIt row_it=lp.addRow();
10.247 - typename std::vector< std::pair<ColIt, double> > row;
10.248 - //std::cout << "node:" <<g.id(n)<<std::endl;
10.249 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
10.250 - if (coeffs[e]!=0) {
10.251 - //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
10.252 - row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
10.253 - }
10.254 - }
10.255 - //std::cout << std::endl;
10.256 - //std::cout << " " << g.id(n) << " " << row.size() << std::endl;
10.257 - lp.setRowCoeffs(row_it, row.begin(), row.end());
10.258 - lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
10.259 - }
10.260 - lp.solveSimplex();
10.261 - //std::cout << lp.colNum() << std::endl;
10.262 - //std::cout << lp.rowNum() << std::endl;
10.263 - //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
10.264 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e)
10.265 - flow.set(e, lp_flow[e]);
10.266 - }
10.267 - };
10.268 -
10.269 -} // namespace lemon
10.270 -
10.271 -#endif //LEMON_MIN_COST_GEN_FLOW_H
11.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
11.2 +++ b/src/work/athos/lp_old/expression.h Wed Mar 23 09:49:41 2005 +0000
11.3 @@ -0,0 +1,197 @@
11.4 +// -*- c++ -*-
11.5 +#ifndef LEMON_EXPRESSION_H
11.6 +#define LEMON_EXPRESSION_H
11.7 +
11.8 +#include <iostream>
11.9 +#include <map>
11.10 +#include <limits>
11.11 +
11.12 +namespace lemon {
11.13 +
11.14 + /*! \brief Linear expression
11.15 +
11.16 + \c Expr<_Col,_Value> implements a class of linear expressions with the
11.17 + operations of addition and multiplication with scalar.
11.18 +
11.19 + \author Marton Makai
11.20 + */
11.21 + template <typename _Col, typename _Value>
11.22 + class Expr;
11.23 +
11.24 + template <typename _Col, typename _Value>
11.25 + class Expr {
11.26 +// protected:
11.27 + public:
11.28 + typedef
11.29 + typename std::map<_Col, _Value> Data;
11.30 + Data data;
11.31 + public:
11.32 + void simplify() {
11.33 + for (typename Data::iterator i=data.begin();
11.34 + i!=data.end(); ++i) {
11.35 + if ((*i).second==0) data.erase(i);
11.36 + }
11.37 + }
11.38 + Expr() { }
11.39 + Expr(_Col _col) {
11.40 + data.insert(std::make_pair(_col, 1));
11.41 + }
11.42 + Expr& operator*=(_Value _value) {
11.43 + for (typename Data::iterator i=data.begin();
11.44 + i!=data.end(); ++i) {
11.45 + (*i).second *= _value;
11.46 + }
11.47 + simplify();
11.48 + return *this;
11.49 + }
11.50 + Expr& operator+=(const Expr<_Col, _Value>& expr) {
11.51 + for (typename Data::const_iterator j=expr.data.begin();
11.52 + j!=expr.data.end(); ++j) {
11.53 + typename Data::iterator i=data.find((*j).first);
11.54 + if (i==data.end()) {
11.55 + data.insert(std::make_pair((*j).first, (*j).second));
11.56 + } else {
11.57 + (*i).second+=(*j).second;
11.58 + }
11.59 + }
11.60 + simplify();
11.61 + return *this;
11.62 + }
11.63 + Expr& operator-=(const Expr<_Col, _Value>& expr) {
11.64 + for (typename Data::const_iterator j=expr.data.begin();
11.65 + j!=expr.data.end(); ++j) {
11.66 + typename Data::iterator i=data.find((*j).first);
11.67 + if (i==data.end()) {
11.68 + data.insert(std::make_pair((*j).first, -(*j).second));
11.69 + } else {
11.70 + (*i).second+=-(*j).second;
11.71 + }
11.72 + }
11.73 + simplify();
11.74 + return *this;
11.75 + }
11.76 + template <typename _C, typename _V>
11.77 + friend std::ostream& operator<<(std::ostream& os,
11.78 + const Expr<_C, _V>& expr);
11.79 + };
11.80 +
11.81 + template <typename _Col, typename _Value>
11.82 + Expr<_Col, _Value> operator*(_Value _value, _Col _col) {
11.83 + Expr<_Col, _Value> tmp(_col);
11.84 + tmp*=_value;
11.85 + tmp.simplify();
11.86 + return tmp;
11.87 + }
11.88 +
11.89 + template <typename _Col, typename _Value>
11.90 + Expr<_Col, _Value> operator*(_Value _value,
11.91 + const Expr<_Col, _Value>& expr) {
11.92 + Expr<_Col, _Value> tmp(expr);
11.93 + tmp*=_value;
11.94 + tmp.simplify();
11.95 + return tmp;
11.96 + }
11.97 +
11.98 + template <typename _Col, typename _Value>
11.99 + Expr<_Col, _Value> operator+(const Expr<_Col, _Value>& expr1,
11.100 + const Expr<_Col, _Value>& expr2) {
11.101 + Expr<_Col, _Value> tmp(expr1);
11.102 + tmp+=expr2;
11.103 + tmp.simplify();
11.104 + return tmp;
11.105 + }
11.106 +
11.107 + template <typename _Col, typename _Value>
11.108 + Expr<_Col, _Value> operator-(const Expr<_Col, _Value>& expr1,
11.109 + const Expr<_Col, _Value>& expr2) {
11.110 + Expr<_Col, _Value> tmp(expr1);
11.111 + tmp-=expr2;
11.112 + tmp.simplify();
11.113 + return tmp;
11.114 + }
11.115 +
11.116 + template <typename _Col, typename _Value>
11.117 + std::ostream& operator<<(std::ostream& os,
11.118 + const Expr<_Col, _Value>& expr) {
11.119 + for (typename Expr<_Col, _Value>::Data::const_iterator i=
11.120 + expr.data.begin();
11.121 + i!=expr.data.end(); ++i) {
11.122 + os << (*i).second << "*" << (*i).first << " ";
11.123 + }
11.124 + return os;
11.125 + }
11.126 +
11.127 + template <typename _Col, typename _Value>
11.128 + class LConstr {
11.129 + // protected:
11.130 + public:
11.131 + Expr<_Col, _Value> expr;
11.132 + _Value lo;
11.133 + public:
11.134 + LConstr(const Expr<_Col, _Value>& _expr, _Value _lo) :
11.135 + expr(_expr), lo(_lo) { }
11.136 + };
11.137 +
11.138 + template <typename _Col, typename _Value>
11.139 + LConstr<_Col, _Value>
11.140 + operator<=(_Value lo, const Expr<_Col, _Value>& expr) {
11.141 + return LConstr<_Col, _Value>(expr, lo);
11.142 + }
11.143 +
11.144 + template <typename _Col, typename _Value>
11.145 + class UConstr {
11.146 + // protected:
11.147 + public:
11.148 + Expr<_Col, _Value> expr;
11.149 + _Value up;
11.150 + public:
11.151 + UConstr(const Expr<_Col, _Value>& _expr, _Value _up) :
11.152 + expr(_expr), up(_up) { }
11.153 + };
11.154 +
11.155 + template <typename _Col, typename _Value>
11.156 + UConstr<_Col, _Value>
11.157 + operator<=(const Expr<_Col, _Value>& expr, _Value up) {
11.158 + return UConstr<_Col, _Value>(expr, up);
11.159 + }
11.160 +
11.161 + template <typename _Col, typename _Value>
11.162 + class Constr {
11.163 + // protected:
11.164 + public:
11.165 + Expr<_Col, _Value> expr;
11.166 + _Value lo, up;
11.167 + public:
11.168 + Constr(const Expr<_Col, _Value>& _expr, _Value _lo, _Value _up) :
11.169 + expr(_expr), lo(_lo), up(_up) { }
11.170 + Constr(const LConstr<_Col, _Value>& _lconstr) :
11.171 + expr(_lconstr.expr),
11.172 + lo(_lconstr.lo),
11.173 + up(std::numeric_limits<_Value>::infinity()) { }
11.174 + Constr(const UConstr<_Col, _Value>& _uconstr) :
11.175 + expr(_uconstr.expr),
11.176 + lo(-std::numeric_limits<_Value>::infinity()),
11.177 + up(_uconstr.up) { }
11.178 + };
11.179 +
11.180 + template <typename _Col, typename _Value>
11.181 + Constr<_Col, _Value>
11.182 + operator<=(const LConstr<_Col, _Value>& lconstr, _Value up) {
11.183 + return Constr<_Col, _Value>(lconstr.expr, lconstr.lo, up);
11.184 + }
11.185 +
11.186 + template <typename _Col, typename _Value>
11.187 + Constr<_Col, _Value>
11.188 + operator<=(_Value lo, const UConstr<_Col, _Value>& uconstr) {
11.189 + return Constr<_Col, _Value>(uconstr.expr, lo, uconstr.up);
11.190 + }
11.191 +
11.192 + template <typename _Col, typename _Value>
11.193 + Constr<_Col, _Value>
11.194 + operator==(const Expr<_Col, _Value>& expr, _Value value) {
11.195 + return Constr<_Col, _Value>(expr, value, value);
11.196 + }
11.197 +
11.198 +} //namespace lemon
11.199 +
11.200 +#endif //LEMON_EXPRESSION_H
12.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
12.2 +++ b/src/work/athos/lp_old/expression_test.cc Wed Mar 23 09:49:41 2005 +0000
12.3 @@ -0,0 +1,23 @@
12.4 +#include <expression.h>
12.5 +#include <iostream>
12.6 +#include <string>
12.7 +
12.8 +using std::cout;
12.9 +using std::endl;
12.10 +using std::string;
12.11 +using namespace lemon;
12.12 +
12.13 +int main() {
12.14 + Expr<string, double> b;
12.15 + cout << b << endl;
12.16 + Expr<string, double> c("f");
12.17 + cout << c << endl;
12.18 + Expr<string, double> d=8.0*string("g");
12.19 + cout << d << endl;
12.20 + c*=5;
12.21 + cout << c << endl;
12.22 + Expr<string, double> e=c;
12.23 + e+=8.9*9.0*string("l");
12.24 + cout << e << endl;
12.25 + cout << c+d << endl;
12.26 +}
13.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
13.2 +++ b/src/work/athos/lp_old/lp_solver_base.h Wed Mar 23 09:49:41 2005 +0000
13.3 @@ -0,0 +1,639 @@
13.4 +// -*- c++ -*-
13.5 +#ifndef LEMON_LP_SOLVER_BASE_H
13.6 +#define LEMON_LP_SOLVER_BASE_H
13.7 +
13.8 +///\ingroup misc
13.9 +///\file
13.10 +
13.11 +// #include <stdio.h>
13.12 +#include <stdlib.h>
13.13 +#include <iostream>
13.14 +#include <map>
13.15 +#include <limits>
13.16 +// #include <stdio>
13.17 +//#include <stdlib>
13.18 +
13.19 +#include <iostream>
13.20 +#include <vector>
13.21 +#include <string>
13.22 +#include <list>
13.23 +#include <memory>
13.24 +#include <utility>
13.25 +
13.26 +#include <lemon/invalid.h>
13.27 +#include <expression.h>
13.28 +//#include <stp.h>
13.29 +//#include <lemon/max_flow.h>
13.30 +//#include <augmenting_flow.h>
13.31 +//#include <iter_map.h>
13.32 +
13.33 +using std::cout;
13.34 +using std::cin;
13.35 +using std::endl;
13.36 +
13.37 +namespace lemon {
13.38 +
13.39 + /// \addtogroup misc
13.40 + /// @{
13.41 +
13.42 + /// \brief A partitioned vector with iterable classes.
13.43 + ///
13.44 + /// This class implements a container in which the data is stored in an
13.45 + /// stl vector, the range is partitioned into sets and each set is
13.46 + /// doubly linked in a list.
13.47 + /// That is, each class is iterable by lemon iterators, and any member of
13.48 + /// the vector can bo moved to an other class.
13.49 + template <typename T>
13.50 + class IterablePartition {
13.51 + protected:
13.52 + struct Node {
13.53 + T data;
13.54 + int prev; //invalid az -1
13.55 + int next;
13.56 + };
13.57 + std::vector<Node> nodes;
13.58 + struct Tip {
13.59 + int first;
13.60 + int last;
13.61 + };
13.62 + std::vector<Tip> tips;
13.63 + public:
13.64 + /// The classes are indexed by integers from \c 0 to \c classNum()-1.
13.65 + int classNum() const { return tips.size(); }
13.66 + /// This lemon style iterator iterates through a class.
13.67 + class Class;
13.68 + /// Constructor. The number of classes is to be given which is fixed
13.69 + /// over the life of the container.
13.70 + /// The partition classes are indexed from 0 to class_num-1.
13.71 + IterablePartition(int class_num) {
13.72 + for (int i=0; i<class_num; ++i) {
13.73 + Tip t;
13.74 + t.first=t.last=-1;
13.75 + tips.push_back(t);
13.76 + }
13.77 + }
13.78 + protected:
13.79 + void befuz(Class it, int class_id) {
13.80 + if (tips[class_id].first==-1) {
13.81 + if (tips[class_id].last==-1) {
13.82 + nodes[it.i].prev=nodes[it.i].next=-1;
13.83 + tips[class_id].first=tips[class_id].last=it.i;
13.84 + }
13.85 + } else {
13.86 + nodes[it.i].prev=tips[class_id].last;
13.87 + nodes[it.i].next=-1;
13.88 + nodes[tips[class_id].last].next=it.i;
13.89 + tips[class_id].last=it.i;
13.90 + }
13.91 + }
13.92 + void kifuz(Class it, int class_id) {
13.93 + if (tips[class_id].first==it.i) {
13.94 + if (tips[class_id].last==it.i) {
13.95 + tips[class_id].first=tips[class_id].last=-1;
13.96 + } else {
13.97 + tips[class_id].first=nodes[it.i].next;
13.98 + nodes[nodes[it.i].next].prev=-1;
13.99 + }
13.100 + } else {
13.101 + if (tips[class_id].last==it.i) {
13.102 + tips[class_id].last=nodes[it.i].prev;
13.103 + nodes[nodes[it.i].prev].next=-1;
13.104 + } else {
13.105 + nodes[nodes[it.i].next].prev=nodes[it.i].prev;
13.106 + nodes[nodes[it.i].prev].next=nodes[it.i].next;
13.107 + }
13.108 + }
13.109 + }
13.110 + public:
13.111 + /// A new element with data \c t is pushed into the vector and into class
13.112 + /// \c class_id.
13.113 + Class push_back(const T& t, int class_id) {
13.114 + Node n;
13.115 + n.data=t;
13.116 + nodes.push_back(n);
13.117 + int i=nodes.size()-1;
13.118 + befuz(i, class_id);
13.119 + return i;
13.120 + }
13.121 + /// A member is moved to an other class.
13.122 + void set(Class it, int old_class_id, int new_class_id) {
13.123 + kifuz(it.i, old_class_id);
13.124 + befuz(it.i, new_class_id);
13.125 + }
13.126 + /// Returns the data pointed by \c it.
13.127 + T& operator[](Class it) { return nodes[it.i].data; }
13.128 + /// Returns the data pointed by \c it.
13.129 + const T& operator[](Class it) const { return nodes[it.i].data; }
13.130 + ///.
13.131 + class Class {
13.132 + friend class IterablePartition;
13.133 + protected:
13.134 + int i;
13.135 + public:
13.136 + /// Default constructor.
13.137 + Class() { }
13.138 + /// This constructor constructs an iterator which points
13.139 + /// to the member of th container indexed by the integer _i.
13.140 + Class(const int& _i) : i(_i) { }
13.141 + /// Invalid constructor.
13.142 + Class(const Invalid&) : i(-1) { }
13.143 + friend bool operator<(const Class& x, const Class& y);
13.144 + friend std::ostream& operator<<(std::ostream& os,
13.145 + const Class& it);
13.146 + bool operator==(const Class& node) const {return i == node.i;}
13.147 + bool operator!=(const Class& node) const {return i != node.i;}
13.148 + };
13.149 + friend bool operator<(const Class& x, const Class& y) {
13.150 + return (x.i < y.i);
13.151 + }
13.152 + friend std::ostream& operator<<(std::ostream& os,
13.153 + const Class& it) {
13.154 + os << it.i;
13.155 + return os;
13.156 + }
13.157 + /// First member of class \c class_id.
13.158 + Class& first(Class& it, int class_id) const {
13.159 + it.i=tips[class_id].first;
13.160 + return it;
13.161 + }
13.162 + /// Next member.
13.163 + Class& next(Class& it) const {
13.164 + it.i=nodes[it.i].next;
13.165 + return it;
13.166 + }
13.167 + /// True iff the iterator is valid.
13.168 + bool valid(const Class& it) const { return it.i!=-1; }
13.169 +
13.170 + class ClassIt : public Class {
13.171 + const IterablePartition* iterable_partition;
13.172 + public:
13.173 + ClassIt() { }
13.174 + ClassIt(Invalid i) : Class(i) { }
13.175 + ClassIt(const IterablePartition& _iterable_partition,
13.176 + const int& i) : iterable_partition(&_iterable_partition) {
13.177 + _iterable_partition.first(*this, i);
13.178 + }
13.179 + ClassIt(const IterablePartition& _iterable_partition,
13.180 + const Class& _class) :
13.181 + Class(_class), iterable_partition(&_iterable_partition) { }
13.182 + ClassIt& operator++() {
13.183 + iterable_partition->next(*this);
13.184 + return *this;
13.185 + }
13.186 + };
13.187 +
13.188 + };
13.189 +
13.190 +
13.191 + /*! \e
13.192 + \todo kellenene uj iterable structure bele, mert ez nem az igazi
13.193 + \todo A[x,y]-t cserel. Jobboldal, baloldal csere.
13.194 + \todo LEKERDEZESEK!!!
13.195 + \todo DOKSI!!!! Doxygen group!!!
13.196 + The aim of this class is to give a general surface to different
13.197 + solvers, i.e. it makes possible to write algorithms using LP's,
13.198 + in which the solver can be changed to an other one easily.
13.199 + \nosubgrouping
13.200 + */
13.201 + template <typename _Value>
13.202 + class LpSolverBase {
13.203 +
13.204 + /*! @name Uncategorized functions and types (public members)
13.205 + */
13.206 + //@{
13.207 + public:
13.208 +
13.209 + //UNCATEGORIZED
13.210 +
13.211 + /// \e
13.212 + typedef IterablePartition<int> Rows;
13.213 + /// \e
13.214 + typedef IterablePartition<int> Cols;
13.215 + /// \e
13.216 + typedef _Value Value;
13.217 + /// \e
13.218 + typedef Rows::Class Row;
13.219 + /// \e
13.220 + typedef Cols::Class Col;
13.221 + public:
13.222 + /// \e
13.223 + IterablePartition<int> row_iter_map;
13.224 + /// \e
13.225 + IterablePartition<int> col_iter_map;
13.226 + /// \e
13.227 + std::vector<Row> int_row_map;
13.228 + /// \e
13.229 + std::vector<Col> int_col_map;
13.230 + /// \e
13.231 + const int VALID_CLASS;
13.232 + /// \e
13.233 + const int INVALID_CLASS;
13.234 + /// \e
13.235 + static const Value INF;
13.236 + public:
13.237 + /// \e
13.238 + LpSolverBase() : row_iter_map(2),
13.239 + col_iter_map(2),
13.240 + VALID_CLASS(0), INVALID_CLASS(1) { }
13.241 + /// \e
13.242 + virtual ~LpSolverBase() { }
13.243 + //@}
13.244 +
13.245 + /*! @name Medium level interface (public members)
13.246 + These functions appear in the low level and also in the high level
13.247 + interfaces thus these each of these functions have to be implemented
13.248 + only once in the different interfaces.
13.249 + This means that these functions have to be reimplemented for all of the
13.250 + different lp solvers. These are basic functions, and have the same
13.251 + parameter lists in the low and high level interfaces.
13.252 + */
13.253 + //@{
13.254 + public:
13.255 +
13.256 + //UNCATEGORIZED FUNCTIONS
13.257 +
13.258 + /// \e
13.259 + virtual void setMinimize() = 0;
13.260 + /// \e
13.261 + virtual void setMaximize() = 0;
13.262 +
13.263 + //SOLVER FUNCTIONS
13.264 +
13.265 + /// \e
13.266 + virtual void solveSimplex() = 0;
13.267 + /// \e
13.268 + virtual void solvePrimalSimplex() = 0;
13.269 + /// \e
13.270 + virtual void solveDualSimplex() = 0;
13.271 +
13.272 + //SOLUTION RETRIEVING
13.273 +
13.274 + /// \e
13.275 + virtual Value getObjVal() = 0;
13.276 +
13.277 + //OTHER FUNCTIONS
13.278 +
13.279 + /// \e
13.280 + virtual int rowNum() const = 0;
13.281 + /// \e
13.282 + virtual int colNum() const = 0;
13.283 + /// \e
13.284 + virtual int warmUp() = 0;
13.285 + /// \e
13.286 + virtual void printWarmUpStatus(int i) = 0;
13.287 + /// \e
13.288 + virtual int getPrimalStatus() = 0;
13.289 + /// \e
13.290 + virtual void printPrimalStatus(int i) = 0;
13.291 + /// \e
13.292 + virtual int getDualStatus() = 0;
13.293 + /// \e
13.294 + virtual void printDualStatus(int i) = 0;
13.295 + /// Returns the status of the slack variable assigned to row \c row.
13.296 + virtual int getRowStat(const Row& row) = 0;
13.297 + /// \e
13.298 + virtual void printRowStatus(int i) = 0;
13.299 + /// Returns the status of the variable assigned to column \c col.
13.300 + virtual int getColStat(const Col& col) = 0;
13.301 + /// \e
13.302 + virtual void printColStatus(int i) = 0;
13.303 +
13.304 + //@}
13.305 +
13.306 + /*! @name Low level interface (protected members)
13.307 + Problem manipulating functions in the low level interface
13.308 + */
13.309 + //@{
13.310 + protected:
13.311 +
13.312 + //MATRIX MANIPULATING FUNCTIONS
13.313 +
13.314 + /// \e
13.315 + virtual int _addCol() = 0;
13.316 + /// \e
13.317 + virtual int _addRow() = 0;
13.318 + /// \e
13.319 + virtual void _eraseCol(int i) = 0;
13.320 + /// \e
13.321 + virtual void _eraseRow(int i) = 0;
13.322 + /// \e
13.323 + virtual void _setRowCoeffs(int i,
13.324 + const std::vector<std::pair<int, Value> >& coeffs) = 0;
13.325 + /// \e
13.326 + /// This routine modifies \c coeffs only by the \c push_back method.
13.327 + virtual void _getRowCoeffs(int i,
13.328 + std::vector<std::pair<int, Value> >& coeffs) = 0;
13.329 + /// \e
13.330 + virtual void _setColCoeffs(int i,
13.331 + const std::vector<std::pair<int, Value> >& coeffs) = 0;
13.332 + /// \e
13.333 + /// This routine modifies \c coeffs only by the \c push_back method.
13.334 + virtual void _getColCoeffs(int i,
13.335 + std::vector<std::pair<int, Value> >& coeffs) = 0;
13.336 + /// \e
13.337 + virtual void _setCoeff(int col, int row, Value value) = 0;
13.338 + /// \e
13.339 + virtual Value _getCoeff(int col, int row) = 0;
13.340 + // public:
13.341 + // /// \e
13.342 + // enum Bound { FREE, LOWER, UPPER, DOUBLE, FIXED };
13.343 + protected:
13.344 + /// \e
13.345 + /// The lower bound of a variable (column) have to be given by an
13.346 + /// extended number of type Value, i.e. a finite number of type
13.347 + /// Value or -INF.
13.348 + virtual void _setColLowerBound(int i, Value value) = 0;
13.349 + /// \e
13.350 + /// The lower bound of a variable (column) is an
13.351 + /// extended number of type Value, i.e. a finite number of type
13.352 + /// Value or -INF.
13.353 + virtual Value _getColLowerBound(int i) = 0;
13.354 + /// \e
13.355 + /// The upper bound of a variable (column) have to be given by an
13.356 + /// extended number of type Value, i.e. a finite number of type
13.357 + /// Value or INF.
13.358 + virtual void _setColUpperBound(int i, Value value) = 0;
13.359 + /// \e
13.360 + /// The upper bound of a variable (column) is an
13.361 + /// extended number of type Value, i.e. a finite number of type
13.362 + /// Value or INF.
13.363 + virtual Value _getColUpperBound(int i) = 0;
13.364 + /// \e
13.365 + /// The lower bound of a linear expression (row) have to be given by an
13.366 + /// extended number of type Value, i.e. a finite number of type
13.367 + /// Value or -INF.
13.368 + virtual void _setRowLowerBound(int i, Value value) = 0;
13.369 + /// \e
13.370 + /// The lower bound of a linear expression (row) is an
13.371 + /// extended number of type Value, i.e. a finite number of type
13.372 + /// Value or -INF.
13.373 + virtual Value _getRowLowerBound(int i) = 0;
13.374 + /// \e
13.375 + /// The upper bound of a linear expression (row) have to be given by an
13.376 + /// extended number of type Value, i.e. a finite number of type
13.377 + /// Value or INF.
13.378 + virtual void _setRowUpperBound(int i, Value value) = 0;
13.379 + /// \e
13.380 + /// The upper bound of a linear expression (row) is an
13.381 + /// extended number of type Value, i.e. a finite number of type
13.382 + /// Value or INF.
13.383 + virtual Value _getRowUpperBound(int i) = 0;
13.384 + /// \e
13.385 + virtual void _setObjCoeff(int i, Value obj_coef) = 0;
13.386 + /// \e
13.387 + virtual Value _getObjCoeff(int i) = 0;
13.388 +
13.389 + //SOLUTION RETRIEVING
13.390 +
13.391 + /// \e
13.392 + virtual Value _getPrimal(int i) = 0;
13.393 + //@}
13.394 +
13.395 + /*! @name High level interface (public members)
13.396 + Problem manipulating functions in the high level interface
13.397 + */
13.398 + //@{
13.399 + public:
13.400 +
13.401 + //MATRIX MANIPULATING FUNCTIONS
13.402 +
13.403 + /// \e
13.404 + Col addCol() {
13.405 + int i=_addCol();
13.406 + Col col;
13.407 + col_iter_map.first(col, INVALID_CLASS);
13.408 + if (col_iter_map.valid(col)) { //van hasznalhato hely
13.409 + col_iter_map.set(col, INVALID_CLASS, VALID_CLASS);
13.410 + col_iter_map[col]=i;
13.411 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
13.412 + col=col_iter_map.push_back(i, VALID_CLASS);
13.413 + }
13.414 + int_col_map.push_back(col);
13.415 + return col;
13.416 + }
13.417 + /// \e
13.418 + Row addRow() {
13.419 + int i=_addRow();
13.420 + Row row;
13.421 + row_iter_map.first(row, INVALID_CLASS);
13.422 + if (row_iter_map.valid(row)) { //van hasznalhato hely
13.423 + row_iter_map.set(row, INVALID_CLASS, VALID_CLASS);
13.424 + row_iter_map[row]=i;
13.425 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
13.426 + row=row_iter_map.push_back(i, VALID_CLASS);
13.427 + }
13.428 + int_row_map.push_back(row);
13.429 + return row;
13.430 + }
13.431 + /// \e
13.432 + void eraseCol(const Col& col) {
13.433 + col_iter_map.set(col, VALID_CLASS, INVALID_CLASS);
13.434 + int cols[2];
13.435 + cols[1]=col_iter_map[col];
13.436 + _eraseCol(cols[1]);
13.437 + col_iter_map[col]=0; //glpk specifikus, de kell ez??
13.438 + Col it;
13.439 + for (col_iter_map.first(it, VALID_CLASS);
13.440 + col_iter_map.valid(it); col_iter_map.next(it)) {
13.441 + if (col_iter_map[it]>cols[1]) --col_iter_map[it];
13.442 + }
13.443 + int_col_map.erase(int_col_map.begin()+cols[1]);
13.444 + }
13.445 + /// \e
13.446 + void eraseRow(const Row& row) {
13.447 + row_iter_map.set(row, VALID_CLASS, INVALID_CLASS);
13.448 + int rows[2];
13.449 + rows[1]=row_iter_map[row];
13.450 + _eraseRow(rows[1]);
13.451 + row_iter_map[row]=0; //glpk specifikus, de kell ez??
13.452 + Row it;
13.453 + for (row_iter_map.first(it, VALID_CLASS);
13.454 + row_iter_map.valid(it); row_iter_map.next(it)) {
13.455 + if (row_iter_map[it]>rows[1]) --row_iter_map[it];
13.456 + }
13.457 + int_row_map.erase(int_row_map.begin()+rows[1]);
13.458 + }
13.459 + /// \e
13.460 + void setCoeff(Col col, Row row, Value value) {
13.461 + _setCoeff(col_iter_map[col], row_iter_map[row], value);
13.462 + }
13.463 + /// \e
13.464 + Value getCoeff(Col col, Row row) {
13.465 + return _getCoeff(col_iter_map[col], row_iter_map[row], value);
13.466 + }
13.467 + /// \e
13.468 + void setColLowerBound(Col col, Value lo) {
13.469 + _setColLowerBound(col_iter_map[col], lo);
13.470 + }
13.471 + /// \e
13.472 + Value getColLowerBound(Col col) {
13.473 + return _getColLowerBound(col_iter_map[col]);
13.474 + }
13.475 + /// \e
13.476 + void setColUpperBound(Col col, Value up) {
13.477 + _setColUpperBound(col_iter_map[col], up);
13.478 + }
13.479 + /// \e
13.480 + Value getColUpperBound(Col col) {
13.481 + return _getColUpperBound(col_iter_map[col]);
13.482 + }
13.483 + /// \e
13.484 + void setRowLowerBound(Row row, Value lo) {
13.485 + _setRowLowerBound(row_iter_map[row], lo);
13.486 + }
13.487 + /// \e
13.488 + Value getRowLowerBound(Row row) {
13.489 + return _getRowLowerBound(row_iter_map[row]);
13.490 + }
13.491 + /// \e
13.492 + void setRowUpperBound(Row row, Value up) {
13.493 + _setRowUpperBound(row_iter_map[row], up);
13.494 + }
13.495 + /// \e
13.496 + Value getRowUpperBound(Row row) {
13.497 + return _getRowUpperBound(row_iter_map[row]);
13.498 + }
13.499 + /// \e
13.500 + void setObjCoeff(const Col& col, Value obj_coef) {
13.501 + _setObjCoeff(col_iter_map[col], obj_coef);
13.502 + }
13.503 + /// \e
13.504 + Value getObjCoeff(const Col& col) {
13.505 + return _getObjCoeff(col_iter_map[col]);
13.506 + }
13.507 +
13.508 + //SOLUTION RETRIEVING FUNCTIONS
13.509 +
13.510 + /// \e
13.511 + Value getPrimal(const Col& col) {
13.512 + return _getPrimal(col_iter_map[col]);
13.513 + }
13.514 +
13.515 + //@}
13.516 +
13.517 + /*! @name User friend interface
13.518 + Problem manipulating functions in the user friend interface
13.519 + */
13.520 + //@{
13.521 +
13.522 + //EXPRESSION TYPES
13.523 +
13.524 + /// \e
13.525 + typedef Expr<Col, Value> Expression;
13.526 + /// \e
13.527 + typedef Expr<Row, Value> DualExpression;
13.528 + /// \e
13.529 + typedef Constr<Col, Value> Constraint;
13.530 +
13.531 + //MATRIX MANIPULATING FUNCTIONS
13.532 +
13.533 + /// \e
13.534 + void setRowCoeffs(Row row, const Expression& expr) {
13.535 + std::vector<std::pair<int, Value> > row_coeffs;
13.536 + for(typename Expression::Data::const_iterator i=expr.data.begin();
13.537 + i!=expr.data.end(); ++i) {
13.538 + row_coeffs.push_back(std::make_pair
13.539 + (col_iter_map[(*i).first], (*i).second));
13.540 + }
13.541 + _setRowCoeffs(row_iter_map[row], row_coeffs);
13.542 + }
13.543 + /// \e
13.544 + void setRow(Row row, const Constraint& constr) {
13.545 + setRowCoeffs(row, constr.expr);
13.546 + setRowLowerBound(row, constr.lo);
13.547 + setRowUpperBound(row, constr.up);
13.548 + }
13.549 + /// \e
13.550 + Row addRow(const Constraint& constr) {
13.551 + Row row=addRow();
13.552 + setRowCoeffs(row, constr.expr);
13.553 + setRowLowerBound(row, constr.lo);
13.554 + setRowUpperBound(row, constr.up);
13.555 + return row;
13.556 + }
13.557 + /// \e
13.558 + /// This routine modifies \c expr by only adding to it.
13.559 + void getRowCoeffs(Row row, Expression& expr) {
13.560 + std::vector<std::pair<int, Value> > row_coeffs;
13.561 + _getRowCoeffs(row_iter_map[row], row_coeffs);
13.562 + for(typename std::vector<std::pair<int, Value> >::const_iterator
13.563 + i=row_coeffs.begin(); i!=row_coeffs.end(); ++i) {
13.564 + expr+= (*i).second*int_col_map[(*i).first];
13.565 + }
13.566 + }
13.567 + /// \e
13.568 + void setColCoeffs(Col col, const DualExpression& expr) {
13.569 + std::vector<std::pair<int, Value> > col_coeffs;
13.570 + for(typename DualExpression::Data::const_iterator i=expr.data.begin();
13.571 + i!=expr.data.end(); ++i) {
13.572 + col_coeffs.push_back(std::make_pair
13.573 + (row_iter_map[(*i).first], (*i).second));
13.574 + }
13.575 + _setColCoeffs(col_iter_map[col], col_coeffs);
13.576 + }
13.577 + /// \e
13.578 + /// This routine modifies \c expr by only adding to it.
13.579 + void getColCoeffs(Col col, DualExpression& expr) {
13.580 + std::vector<std::pair<int, Value> > col_coeffs;
13.581 + _getColCoeffs(col_iter_map[col], col_coeffs);
13.582 + for(typename std::vector<std::pair<int, Value> >::const_iterator
13.583 + i=col_coeffs.begin(); i!=col_coeffs.end(); ++i) {
13.584 + expr+= (*i).second*int_row_map[(*i).first];
13.585 + }
13.586 + }
13.587 + /// \e
13.588 + void setObjCoeffs(const Expression& expr) {
13.589 + // writing zero everywhere
13.590 + for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
13.591 + setObjCoeff(it, 0.0);
13.592 + // writing the data needed
13.593 + for(typename Expression::Data::const_iterator i=expr.data.begin();
13.594 + i!=expr.data.end(); ++i) {
13.595 + setObjCoeff((*i).first, (*i).second);
13.596 + }
13.597 + }
13.598 + /// \e
13.599 + /// This routine modifies \c expr by only adding to it.
13.600 + void getObjCoeffs(Expression& expr) {
13.601 + for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
13.602 + expr+=getObjCoeff(it)*it;
13.603 + }
13.604 + //@}
13.605 +
13.606 +
13.607 + /*! @name MIP functions and types (public members)
13.608 + */
13.609 + //@{
13.610 + public:
13.611 + /// \e
13.612 + virtual void solveBandB() = 0;
13.613 + /// \e
13.614 + virtual void setLP() = 0;
13.615 + /// \e
13.616 + virtual void setMIP() = 0;
13.617 + protected:
13.618 + /// \e
13.619 + virtual void _setColCont(int i) = 0;
13.620 + /// \e
13.621 + virtual void _setColInt(int i) = 0;
13.622 + /// \e
13.623 + virtual Value _getMIPPrimal(int i) = 0;
13.624 + public:
13.625 + /// \e
13.626 + void setColCont(Col col) {
13.627 + _setColCont(col_iter_map[col]);
13.628 + }
13.629 + /// \e
13.630 + void setColInt(Col col) {
13.631 + _setColInt(col_iter_map[col]);
13.632 + }
13.633 + /// \e
13.634 + Value getMIPPrimal(Col col) {
13.635 + return _getMIPPrimal(col_iter_map[col]);
13.636 + }
13.637 + //@}
13.638 + };
13.639 +
13.640 +} //namespace lemon
13.641 +
13.642 +#endif //LEMON_LP_SOLVER_BASE_H
14.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
14.2 +++ b/src/work/athos/lp_old/lp_solver_glpk.h Wed Mar 23 09:49:41 2005 +0000
14.3 @@ -0,0 +1,545 @@
14.4 +// -*- c++ -*-
14.5 +#ifndef LEMON_LP_SOLVER_GLPK_H
14.6 +#define LEMON_LP_SOLVER_GLPK_H
14.7 +
14.8 +///\ingroup misc
14.9 +///\file
14.10 +
14.11 +// #include <stdio.h>
14.12 +/* #include <stdlib.h> */
14.13 +/* #include <iostream> */
14.14 +/* #include <map> */
14.15 +/* #include <limits> */
14.16 +// #include <stdio>
14.17 +//#include <stdlib>
14.18 +extern "C" {
14.19 +#include "glpk.h"
14.20 +}
14.21 +
14.22 +/* #include <iostream> */
14.23 +/* #include <vector> */
14.24 +/* #include <string> */
14.25 +/* #include <list> */
14.26 +/* #include <memory> */
14.27 +/* #include <utility> */
14.28 +
14.29 +//#include <lemon/invalid.h>
14.30 +//#include <expression.h>
14.31 +#include <lp_solver_base.h>
14.32 +//#include <stp.h>
14.33 +//#include <lemon/max_flow.h>
14.34 +//#include <augmenting_flow.h>
14.35 +//#include <iter_map.h>
14.36 +
14.37 +using std::cout;
14.38 +using std::cin;
14.39 +using std::endl;
14.40 +
14.41 +namespace lemon {
14.42 +
14.43 +
14.44 + template <typename Value>
14.45 + const Value LpSolverBase<Value>::INF=std::numeric_limits<Value>::infinity();
14.46 +
14.47 +
14.48 + /// \brief Wrapper for GLPK solver
14.49 + ///
14.50 + /// This class implements a lemon wrapper for GLPK.
14.51 + class LpGlpk : public LpSolverBase<double> {
14.52 + public:
14.53 + typedef LpSolverBase<double> Parent;
14.54 +
14.55 + public:
14.56 + /// \e
14.57 + LPX* lp;
14.58 +
14.59 + public:
14.60 + /// \e
14.61 + LpGlpk() : Parent(),
14.62 + lp(lpx_create_prob()) {
14.63 + int_row_map.push_back(Row());
14.64 + int_col_map.push_back(Col());
14.65 + lpx_set_int_parm(lp, LPX_K_DUAL, 1);
14.66 + }
14.67 + /// \e
14.68 + ~LpGlpk() {
14.69 + lpx_delete_prob(lp);
14.70 + }
14.71 +
14.72 + //MATRIX INDEPEDENT MANIPULATING FUNCTIONS
14.73 +
14.74 + /// \e
14.75 + void setMinimize() {
14.76 + lpx_set_obj_dir(lp, LPX_MIN);
14.77 + }
14.78 + /// \e
14.79 + void setMaximize() {
14.80 + lpx_set_obj_dir(lp, LPX_MAX);
14.81 + }
14.82 +
14.83 + //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS
14.84 +
14.85 + protected:
14.86 + /// \e
14.87 + int _addCol() {
14.88 + int i=lpx_add_cols(lp, 1);
14.89 + _setColLowerBound(i, -INF);
14.90 + _setColUpperBound(i, INF);
14.91 + return i;
14.92 + }
14.93 + /// \e
14.94 + int _addRow() {
14.95 + int i=lpx_add_rows(lp, 1);
14.96 + return i;
14.97 + }
14.98 + /// \e
14.99 + virtual void _setRowCoeffs(int i,
14.100 + const std::vector<std::pair<int, double> >& coeffs) {
14.101 + int mem_length=1+colNum();
14.102 + int* indices = new int[mem_length];
14.103 + double* doubles = new double[mem_length];
14.104 + int length=0;
14.105 + for (std::vector<std::pair<int, double> >::
14.106 + const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
14.107 + ++length;
14.108 + indices[length]=it->first;
14.109 + doubles[length]=it->second;
14.110 + }
14.111 + lpx_set_mat_row(lp, i, length, indices, doubles);
14.112 + delete [] indices;
14.113 + delete [] doubles;
14.114 + }
14.115 + /// \e
14.116 + virtual void _getRowCoeffs(int i,
14.117 + std::vector<std::pair<int, double> >& coeffs) {
14.118 + int mem_length=1+colNum();
14.119 + int* indices = new int[mem_length];
14.120 + double* doubles = new double[mem_length];
14.121 + int length=lpx_get_mat_row(lp, i, indices, doubles);
14.122 + for (int i=1; i<=length; ++i) {
14.123 + coeffs.push_back(std::make_pair(indices[i], doubles[i]));
14.124 + }
14.125 + delete [] indices;
14.126 + delete [] doubles;
14.127 + }
14.128 + /// \e
14.129 + virtual void _setColCoeffs(int i,
14.130 + const std::vector<std::pair<int, double> >& coeffs) {
14.131 + int mem_length=1+rowNum();
14.132 + int* indices = new int[mem_length];
14.133 + double* doubles = new double[mem_length];
14.134 + int length=0;
14.135 + for (std::vector<std::pair<int, double> >::
14.136 + const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
14.137 + ++length;
14.138 + indices[length]=it->first;
14.139 + doubles[length]=it->second;
14.140 + }
14.141 + lpx_set_mat_col(lp, i, length, indices, doubles);
14.142 + delete [] indices;
14.143 + delete [] doubles;
14.144 + }
14.145 + /// \e
14.146 + virtual void _getColCoeffs(int i,
14.147 + std::vector<std::pair<int, double> >& coeffs) {
14.148 + int mem_length=1+rowNum();
14.149 + int* indices = new int[mem_length];
14.150 + double* doubles = new double[mem_length];
14.151 + int length=lpx_get_mat_col(lp, i, indices, doubles);
14.152 + for (int i=1; i<=length; ++i) {
14.153 + coeffs.push_back(std::make_pair(indices[i], doubles[i]));
14.154 + }
14.155 + delete [] indices;
14.156 + delete [] doubles;
14.157 + }
14.158 + /// \e
14.159 + virtual void _eraseCol(int i) {
14.160 + int cols[2];
14.161 + cols[1]=i;
14.162 + lpx_del_cols(lp, 1, cols);
14.163 + }
14.164 + virtual void _eraseRow(int i) {
14.165 + int rows[2];
14.166 + rows[1]=i;
14.167 + lpx_del_rows(lp, 1, rows);
14.168 + }
14.169 + void _setCoeff(int col, int row, double value) {
14.170 + ///FIXME Of course this is not efficient at all, but GLPK knows not more.
14.171 + int change_this;
14.172 + bool get_set_row;
14.173 + //The only thing we can do is optimize on whether working with a row
14.174 + //or a coloumn
14.175 + int row_num = rowNum();
14.176 + int col_num = colNum();
14.177 + if (col_num<row_num){
14.178 + //There are more rows than coloumns
14.179 + get_set_row=true;
14.180 + int mem_length=1+row_num;
14.181 + int* indices = new int[mem_length];
14.182 + double* doubles = new double[mem_length];
14.183 + int length=lpx_get_mat_col(lp, i, indices, doubles);
14.184 + }else{
14.185 + get_set_row=false;
14.186 + int mem_length=1+col_num;
14.187 + int* indices = new int[mem_length];
14.188 + double* doubles = new double[mem_length];
14.189 + int length=lpx_get_mat_row(lp, i, indices, doubles);
14.190 + }
14.191 + //Itten
14.192 +int* indices = new int[mem_length];
14.193 + double* doubles = new double[mem_length];
14.194 + int length=lpx_get_mat_col(lp, i, indices, doubles);
14.195 +
14.196 + delete [] indices;
14.197 + delete [] doubles;
14.198 +
14.199 + }
14.200 + double _getCoeff(int col, int row) {
14.201 + /// FIXME not yet implemented
14.202 + return 0.0;
14.203 + }
14.204 + virtual void _setColLowerBound(int i, double lo) {
14.205 + if (lo==INF) {
14.206 + //FIXME error
14.207 + }
14.208 + int b=lpx_get_col_type(lp, i);
14.209 + double up=lpx_get_col_ub(lp, i);
14.210 + if (lo==-INF) {
14.211 + switch (b) {
14.212 + case LPX_FR:
14.213 + case LPX_LO:
14.214 + lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
14.215 + break;
14.216 + case LPX_UP:
14.217 + break;
14.218 + case LPX_DB:
14.219 + case LPX_FX:
14.220 + lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
14.221 + break;
14.222 + default: ;
14.223 + //FIXME error
14.224 + }
14.225 + } else {
14.226 + switch (b) {
14.227 + case LPX_FR:
14.228 + case LPX_LO:
14.229 + lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
14.230 + break;
14.231 + case LPX_UP:
14.232 + case LPX_DB:
14.233 + case LPX_FX:
14.234 + if (lo==up)
14.235 + lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
14.236 + else
14.237 + lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
14.238 + break;
14.239 + default: ;
14.240 + //FIXME error
14.241 + }
14.242 + }
14.243 + }
14.244 + virtual double _getColLowerBound(int i) {
14.245 + int b=lpx_get_col_type(lp, i);
14.246 + switch (b) {
14.247 + case LPX_FR:
14.248 + return -INF;
14.249 + case LPX_LO:
14.250 + return lpx_get_col_lb(lp, i);
14.251 + case LPX_UP:
14.252 + return -INF;
14.253 + case LPX_DB:
14.254 + case LPX_FX:
14.255 + return lpx_get_col_lb(lp, i);
14.256 + default: ;
14.257 + //FIXME error
14.258 + return 0.0;
14.259 + }
14.260 + }
14.261 + virtual void _setColUpperBound(int i, double up) {
14.262 + if (up==-INF) {
14.263 + //FIXME error
14.264 + }
14.265 + int b=lpx_get_col_type(lp, i);
14.266 + double lo=lpx_get_col_lb(lp, i);
14.267 + if (up==INF) {
14.268 + switch (b) {
14.269 + case LPX_FR:
14.270 + case LPX_LO:
14.271 + break;
14.272 + case LPX_UP:
14.273 + lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
14.274 + break;
14.275 + case LPX_DB:
14.276 + case LPX_FX:
14.277 + lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
14.278 + break;
14.279 + default: ;
14.280 + //FIXME error
14.281 + }
14.282 + } else {
14.283 + switch (b) {
14.284 + case LPX_FR:
14.285 + lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
14.286 + case LPX_LO:
14.287 + if (lo==up)
14.288 + lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
14.289 + else
14.290 + lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
14.291 + break;
14.292 + case LPX_UP:
14.293 + lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
14.294 + break;
14.295 + case LPX_DB:
14.296 + case LPX_FX:
14.297 + if (lo==up)
14.298 + lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
14.299 + else
14.300 + lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
14.301 + break;
14.302 + default: ;
14.303 + //FIXME error
14.304 + }
14.305 + }
14.306 + }
14.307 + virtual double _getColUpperBound(int i) {
14.308 + int b=lpx_get_col_type(lp, i);
14.309 + switch (b) {
14.310 + case LPX_FR:
14.311 + case LPX_LO:
14.312 + return INF;
14.313 + case LPX_UP:
14.314 + case LPX_DB:
14.315 + case LPX_FX:
14.316 + return lpx_get_col_ub(lp, i);
14.317 + default: ;
14.318 + //FIXME error
14.319 + return 0.0;
14.320 + }
14.321 + }
14.322 + virtual void _setRowLowerBound(int i, double lo) {
14.323 + if (lo==INF) {
14.324 + //FIXME error
14.325 + }
14.326 + int b=lpx_get_row_type(lp, i);
14.327 + double up=lpx_get_row_ub(lp, i);
14.328 + if (lo==-INF) {
14.329 + switch (b) {
14.330 + case LPX_FR:
14.331 + case LPX_LO:
14.332 + lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
14.333 + break;
14.334 + case LPX_UP:
14.335 + break;
14.336 + case LPX_DB:
14.337 + case LPX_FX:
14.338 + lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
14.339 + break;
14.340 + default: ;
14.341 + //FIXME error
14.342 + }
14.343 + } else {
14.344 + switch (b) {
14.345 + case LPX_FR:
14.346 + case LPX_LO:
14.347 + lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
14.348 + break;
14.349 + case LPX_UP:
14.350 + case LPX_DB:
14.351 + case LPX_FX:
14.352 + if (lo==up)
14.353 + lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
14.354 + else
14.355 + lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
14.356 + break;
14.357 + default: ;
14.358 + //FIXME error
14.359 + }
14.360 + }
14.361 + }
14.362 + virtual double _getRowLowerBound(int i) {
14.363 + int b=lpx_get_row_type(lp, i);
14.364 + switch (b) {
14.365 + case LPX_FR:
14.366 + return -INF;
14.367 + case LPX_LO:
14.368 + return lpx_get_row_lb(lp, i);
14.369 + case LPX_UP:
14.370 + return -INF;
14.371 + case LPX_DB:
14.372 + case LPX_FX:
14.373 + return lpx_get_row_lb(lp, i);
14.374 + default: ;
14.375 + //FIXME error
14.376 + return 0.0;
14.377 + }
14.378 + }
14.379 + virtual void _setRowUpperBound(int i, double up) {
14.380 + if (up==-INF) {
14.381 + //FIXME error
14.382 + }
14.383 + int b=lpx_get_row_type(lp, i);
14.384 + double lo=lpx_get_row_lb(lp, i);
14.385 + if (up==INF) {
14.386 + switch (b) {
14.387 + case LPX_FR:
14.388 + case LPX_LO:
14.389 + break;
14.390 + case LPX_UP:
14.391 + lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
14.392 + break;
14.393 + case LPX_DB:
14.394 + case LPX_FX:
14.395 + lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
14.396 + break;
14.397 + default: ;
14.398 + //FIXME error
14.399 + }
14.400 + } else {
14.401 + switch (b) {
14.402 + case LPX_FR:
14.403 + lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
14.404 + case LPX_LO:
14.405 + if (lo==up)
14.406 + lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
14.407 + else
14.408 + lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
14.409 + break;
14.410 + case LPX_UP:
14.411 + lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
14.412 + break;
14.413 + case LPX_DB:
14.414 + case LPX_FX:
14.415 + if (lo==up)
14.416 + lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
14.417 + else
14.418 + lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
14.419 + break;
14.420 + default: ;
14.421 + //FIXME error
14.422 + }
14.423 + }
14.424 + }
14.425 + virtual double _getRowUpperBound(int i) {
14.426 + int b=lpx_get_row_type(lp, i);
14.427 + switch (b) {
14.428 + case LPX_FR:
14.429 + case LPX_LO:
14.430 + return INF;
14.431 + case LPX_UP:
14.432 + case LPX_DB:
14.433 + case LPX_FX:
14.434 + return lpx_get_row_ub(lp, i);
14.435 + default: ;
14.436 + //FIXME error
14.437 + return 0.0;
14.438 + }
14.439 + }
14.440 + /// \e
14.441 + virtual double _getObjCoeff(int i) {
14.442 + return lpx_get_obj_coef(lp, i);
14.443 + }
14.444 + /// \e
14.445 + virtual void _setObjCoeff(int i, double obj_coef) {
14.446 + lpx_set_obj_coef(lp, i, obj_coef);
14.447 + }
14.448 + public:
14.449 + /// \e
14.450 + void solveSimplex() { lpx_simplex(lp); }
14.451 + /// \e
14.452 + void solvePrimalSimplex() { lpx_simplex(lp); }
14.453 + /// \e
14.454 + void solveDualSimplex() { lpx_simplex(lp); }
14.455 + protected:
14.456 + virtual double _getPrimal(int i) {
14.457 + return lpx_get_col_prim(lp, i);
14.458 + }
14.459 + public:
14.460 + /// \e
14.461 + double getObjVal() { return lpx_get_obj_val(lp); }
14.462 + /// \e
14.463 + int rowNum() const { return lpx_get_num_rows(lp); }
14.464 + /// \e
14.465 + int colNum() const { return lpx_get_num_cols(lp); }
14.466 + /// \e
14.467 + int warmUp() { return lpx_warm_up(lp); }
14.468 + /// \e
14.469 + void printWarmUpStatus(int i) {
14.470 + switch (i) {
14.471 + case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
14.472 + case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;
14.473 + case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
14.474 + case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
14.475 + }
14.476 + }
14.477 + /// \e
14.478 + int getPrimalStatus() { return lpx_get_prim_stat(lp); }
14.479 + /// \e
14.480 + void printPrimalStatus(int i) {
14.481 + switch (i) {
14.482 + case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
14.483 + case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;
14.484 + case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
14.485 + case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
14.486 + }
14.487 + }
14.488 + /// \e
14.489 + int getDualStatus() { return lpx_get_dual_stat(lp); }
14.490 + /// \e
14.491 + void printDualStatus(int i) {
14.492 + switch (i) {
14.493 + case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
14.494 + case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;
14.495 + case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
14.496 + case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
14.497 + }
14.498 + }
14.499 + /// Returns the status of the slack variable assigned to row \c row.
14.500 + int getRowStat(const Row& row) {
14.501 + return lpx_get_row_stat(lp, row_iter_map[row]);
14.502 + }
14.503 + /// \e
14.504 + void printRowStatus(int i) {
14.505 + switch (i) {
14.506 + case LPX_BS: cout << "LPX_BS" << endl; break;
14.507 + case LPX_NL: cout << "LPX_NL" << endl; break;
14.508 + case LPX_NU: cout << "LPX_NU" << endl; break;
14.509 + case LPX_NF: cout << "LPX_NF" << endl; break;
14.510 + case LPX_NS: cout << "LPX_NS" << endl; break;
14.511 + }
14.512 + }
14.513 + /// Returns the status of the variable assigned to column \c col.
14.514 + int getColStat(const Col& col) {
14.515 + return lpx_get_col_stat(lp, col_iter_map[col]);
14.516 + }
14.517 + /// \e
14.518 + void printColStatus(int i) {
14.519 + switch (i) {
14.520 + case LPX_BS: cout << "LPX_BS" << endl; break;
14.521 + case LPX_NL: cout << "LPX_NL" << endl; break;
14.522 + case LPX_NU: cout << "LPX_NU" << endl; break;
14.523 + case LPX_NF: cout << "LPX_NF" << endl; break;
14.524 + case LPX_NS: cout << "LPX_NS" << endl; break;
14.525 + }
14.526 + }
14.527 +
14.528 + // MIP
14.529 + /// \e
14.530 + void solveBandB() { lpx_integer(lp); }
14.531 + /// \e
14.532 + void setLP() { lpx_set_class(lp, LPX_LP); }
14.533 + /// \e
14.534 + void setMIP() { lpx_set_class(lp, LPX_MIP); }
14.535 + protected:
14.536 + /// \e
14.537 + void _setColCont(int i) { lpx_set_col_kind(lp, i, LPX_CV); }
14.538 + /// \e
14.539 + void _setColInt(int i) { lpx_set_col_kind(lp, i, LPX_IV); }
14.540 + /// \e
14.541 + double _getMIPPrimal(int i) { return lpx_mip_col_val(lp, i); }
14.542 + };
14.543 +
14.544 + /// @}
14.545 +
14.546 +} //namespace lemon
14.547 +
14.548 +#endif //LEMON_LP_SOLVER_GLPK_H
15.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
15.2 +++ b/src/work/athos/lp_old/lp_solver_wrapper.h Wed Mar 23 09:49:41 2005 +0000
15.3 @@ -0,0 +1,431 @@
15.4 +// -*- c++ -*-
15.5 +#ifndef LEMON_LP_SOLVER_WRAPPER_H
15.6 +#define LEMON_LP_SOLVER_WRAPPER_H
15.7 +
15.8 +///\ingroup misc
15.9 +///\file
15.10 +///\brief Dijkstra algorithm.
15.11 +
15.12 +// #include <stdio.h>
15.13 +#include <stdlib.h>
15.14 +// #include <stdio>
15.15 +//#include <stdlib>
15.16 +extern "C" {
15.17 +#include "glpk.h"
15.18 +}
15.19 +
15.20 +#include <iostream>
15.21 +#include <vector>
15.22 +#include <string>
15.23 +#include <list>
15.24 +#include <memory>
15.25 +#include <utility>
15.26 +
15.27 +//#include <sage_graph.h>
15.28 +//#include <lemon/list_graph.h>
15.29 +//#include <lemon/graph_wrapper.h>
15.30 +#include <lemon/invalid.h>
15.31 +//#include <bfs_dfs.h>
15.32 +//#include <stp.h>
15.33 +//#include <lemon/max_flow.h>
15.34 +//#include <augmenting_flow.h>
15.35 +//#include <iter_map.h>
15.36 +
15.37 +using std::cout;
15.38 +using std::cin;
15.39 +using std::endl;
15.40 +
15.41 +namespace lemon {
15.42 +
15.43 +
15.44 + /// \addtogroup misc
15.45 + /// @{
15.46 +
15.47 + /// \brief A partitioned vector with iterable classes.
15.48 + ///
15.49 + /// This class implements a container in which the data is stored in an
15.50 + /// stl vector, the range is partitioned into sets and each set is
15.51 + /// doubly linked in a list.
15.52 + /// That is, each class is iterable by lemon iterators, and any member of
15.53 + /// the vector can bo moved to an other class.
15.54 + template <typename T>
15.55 + class IterablePartition {
15.56 + protected:
15.57 + struct Node {
15.58 + T data;
15.59 + int prev; //invalid az -1
15.60 + int next;
15.61 + };
15.62 + std::vector<Node> nodes;
15.63 + struct Tip {
15.64 + int first;
15.65 + int last;
15.66 + };
15.67 + std::vector<Tip> tips;
15.68 + public:
15.69 + /// The classes are indexed by integers from \c 0 to \c classNum()-1.
15.70 + int classNum() const { return tips.size(); }
15.71 + /// This lemon style iterator iterates through a class.
15.72 + class ClassIt;
15.73 + /// Constructor. The number of classes is to be given which is fixed
15.74 + /// over the life of the container.
15.75 + /// The partition classes are indexed from 0 to class_num-1.
15.76 + IterablePartition(int class_num) {
15.77 + for (int i=0; i<class_num; ++i) {
15.78 + Tip t;
15.79 + t.first=t.last=-1;
15.80 + tips.push_back(t);
15.81 + }
15.82 + }
15.83 + protected:
15.84 + void befuz(ClassIt it, int class_id) {
15.85 + if (tips[class_id].first==-1) {
15.86 + if (tips[class_id].last==-1) {
15.87 + nodes[it.i].prev=nodes[it.i].next=-1;
15.88 + tips[class_id].first=tips[class_id].last=it.i;
15.89 + }
15.90 + } else {
15.91 + nodes[it.i].prev=tips[class_id].last;
15.92 + nodes[it.i].next=-1;
15.93 + nodes[tips[class_id].last].next=it.i;
15.94 + tips[class_id].last=it.i;
15.95 + }
15.96 + }
15.97 + void kifuz(ClassIt it, int class_id) {
15.98 + if (tips[class_id].first==it.i) {
15.99 + if (tips[class_id].last==it.i) {
15.100 + tips[class_id].first=tips[class_id].last=-1;
15.101 + } else {
15.102 + tips[class_id].first=nodes[it.i].next;
15.103 + nodes[nodes[it.i].next].prev=-1;
15.104 + }
15.105 + } else {
15.106 + if (tips[class_id].last==it.i) {
15.107 + tips[class_id].last=nodes[it.i].prev;
15.108 + nodes[nodes[it.i].prev].next=-1;
15.109 + } else {
15.110 + nodes[nodes[it.i].next].prev=nodes[it.i].prev;
15.111 + nodes[nodes[it.i].prev].next=nodes[it.i].next;
15.112 + }
15.113 + }
15.114 + }
15.115 + public:
15.116 + /// A new element with data \c t is pushed into the vector and into class
15.117 + /// \c class_id.
15.118 + ClassIt push_back(const T& t, int class_id) {
15.119 + Node n;
15.120 + n.data=t;
15.121 + nodes.push_back(n);
15.122 + int i=nodes.size()-1;
15.123 + befuz(i, class_id);
15.124 + return i;
15.125 + }
15.126 + /// A member is moved to an other class.
15.127 + void set(ClassIt it, int old_class_id, int new_class_id) {
15.128 + kifuz(it.i, old_class_id);
15.129 + befuz(it.i, new_class_id);
15.130 + }
15.131 + /// Returns the data pointed by \c it.
15.132 + T& operator[](ClassIt it) { return nodes[it.i].data; }
15.133 + /// Returns the data pointed by \c it.
15.134 + const T& operator[](ClassIt it) const { return nodes[it.i].data; }
15.135 + ///.
15.136 + class ClassIt {
15.137 + friend class IterablePartition;
15.138 + protected:
15.139 + int i;
15.140 + public:
15.141 + /// Default constructor.
15.142 + ClassIt() { }
15.143 + /// This constructor constructs an iterator which points
15.144 + /// to the member of th container indexed by the integer _i.
15.145 + ClassIt(const int& _i) : i(_i) { }
15.146 + /// Invalid constructor.
15.147 + ClassIt(const Invalid&) : i(-1) { }
15.148 + };
15.149 + /// First member of class \c class_id.
15.150 + ClassIt& first(ClassIt& it, int class_id) const {
15.151 + it.i=tips[class_id].first;
15.152 + return it;
15.153 + }
15.154 + /// Next member.
15.155 + ClassIt& next(ClassIt& it) const {
15.156 + it.i=nodes[it.i].next;
15.157 + return it;
15.158 + }
15.159 + /// True iff the iterator is valid.
15.160 + bool valid(const ClassIt& it) const { return it.i!=-1; }
15.161 + };
15.162 +
15.163 + /// \brief Wrappers for LP solvers
15.164 + ///
15.165 + /// This class implements a lemon wrapper for glpk.
15.166 + /// Later other LP-solvers will be wrapped into lemon.
15.167 + /// The aim of this class is to give a general surface to different
15.168 + /// solvers, i.e. it makes possible to write algorithms using LP's,
15.169 + /// in which the solver can be changed to an other one easily.
15.170 + class LPSolverWrapper {
15.171 + public:
15.172 +
15.173 +// class Row {
15.174 +// protected:
15.175 +// int i;
15.176 +// public:
15.177 +// Row() { }
15.178 +// Row(const Invalid&) : i(0) { }
15.179 +// Row(const int& _i) : i(_i) { }
15.180 +// operator int() const { return i; }
15.181 +// };
15.182 +// class RowIt : public Row {
15.183 +// public:
15.184 +// RowIt(const Row& row) : Row(row) { }
15.185 +// };
15.186 +
15.187 +// class Col {
15.188 +// protected:
15.189 +// int i;
15.190 +// public:
15.191 +// Col() { }
15.192 +// Col(const Invalid&) : i(0) { }
15.193 +// Col(const int& _i) : i(_i) { }
15.194 +// operator int() const { return i; }
15.195 +// };
15.196 +// class ColIt : public Col {
15.197 +// ColIt(const Col& col) : Col(col) { }
15.198 +// };
15.199 +
15.200 + public:
15.201 + ///.
15.202 + LPX* lp;
15.203 + ///.
15.204 + typedef IterablePartition<int>::ClassIt RowIt;
15.205 + ///.
15.206 + IterablePartition<int> row_iter_map;
15.207 + ///.
15.208 + typedef IterablePartition<int>::ClassIt ColIt;
15.209 + ///.
15.210 + IterablePartition<int> col_iter_map;
15.211 + //std::vector<int> row_id_to_lp_row_id;
15.212 + //std::vector<int> col_id_to_lp_col_id;
15.213 + ///.
15.214 + const int VALID_ID;
15.215 + ///.
15.216 + const int INVALID_ID;
15.217 +
15.218 + public:
15.219 + ///.
15.220 + LPSolverWrapper() : lp(lpx_create_prob()),
15.221 + row_iter_map(2),
15.222 + col_iter_map(2),
15.223 + //row_id_to_lp_row_id(), col_id_to_lp_col_id(),
15.224 + VALID_ID(0), INVALID_ID(1) {
15.225 + lpx_set_int_parm(lp, LPX_K_DUAL, 1);
15.226 + }
15.227 + ///.
15.228 + ~LPSolverWrapper() {
15.229 + lpx_delete_prob(lp);
15.230 + }
15.231 + ///.
15.232 + void setMinimize() {
15.233 + lpx_set_obj_dir(lp, LPX_MIN);
15.234 + }
15.235 + ///.
15.236 + void setMaximize() {
15.237 + lpx_set_obj_dir(lp, LPX_MAX);
15.238 + }
15.239 + ///.
15.240 + ColIt addCol() {
15.241 + int i=lpx_add_cols(lp, 1);
15.242 + ColIt col_it;
15.243 + col_iter_map.first(col_it, INVALID_ID);
15.244 + if (col_iter_map.valid(col_it)) { //van hasznalhato hely
15.245 + col_iter_map.set(col_it, INVALID_ID, VALID_ID);
15.246 + col_iter_map[col_it]=i;
15.247 + //col_id_to_lp_col_id[col_iter_map[col_it]]=i;
15.248 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
15.249 + //col_id_to_lp_col_id.push_back(i);
15.250 + //int j=col_id_to_lp_col_id.size()-1;
15.251 + col_it=col_iter_map.push_back(i, VALID_ID);
15.252 + }
15.253 +// edge_index_map.set(e, i);
15.254 +// lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
15.255 +// lpx_set_obj_coef(lp, i, cost[e]);
15.256 + return col_it;
15.257 + }
15.258 + ///.
15.259 + RowIt addRow() {
15.260 + int i=lpx_add_rows(lp, 1);
15.261 + RowIt row_it;
15.262 + row_iter_map.first(row_it, INVALID_ID);
15.263 + if (row_iter_map.valid(row_it)) { //van hasznalhato hely
15.264 + row_iter_map.set(row_it, INVALID_ID, VALID_ID);
15.265 + row_iter_map[row_it]=i;
15.266 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
15.267 + row_it=row_iter_map.push_back(i, VALID_ID);
15.268 + }
15.269 + return row_it;
15.270 + }
15.271 + //pair<RowIt, double>-bol kell megadni egy std range-et
15.272 + ///.
15.273 + template <typename Begin, typename End>
15.274 + void setColCoeffs(const ColIt& col_it,
15.275 + Begin begin, End end) {
15.276 + int mem_length=1+lpx_get_num_rows(lp);
15.277 + int* indices = new int[mem_length];
15.278 + double* doubles = new double[mem_length];
15.279 + int length=0;
15.280 + for ( ; begin!=end; ++begin) {
15.281 + ++length;
15.282 + indices[length]=row_iter_map[begin->first];
15.283 + doubles[length]=begin->second;
15.284 + }
15.285 + lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
15.286 + delete [] indices;
15.287 + delete [] doubles;
15.288 + }
15.289 + //pair<ColIt, double>-bol kell megadni egy std range-et
15.290 + ///.
15.291 + template <typename Begin, typename End>
15.292 + void setRowCoeffs(const RowIt& row_it,
15.293 + Begin begin, End end) {
15.294 + int mem_length=1+lpx_get_num_cols(lp);
15.295 + int* indices = new int[mem_length];
15.296 + double* doubles = new double[mem_length];
15.297 + int length=0;
15.298 + for ( ; begin!=end; ++begin) {
15.299 + ++length;
15.300 + indices[length]=col_iter_map[begin->first];
15.301 + doubles[length]=begin->second;
15.302 + }
15.303 + lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
15.304 + delete [] indices;
15.305 + delete [] doubles;
15.306 + }
15.307 + ///.
15.308 + void eraseCol(const ColIt& col_it) {
15.309 + col_iter_map.set(col_it, VALID_ID, INVALID_ID);
15.310 + int cols[2];
15.311 + cols[1]=col_iter_map[col_it];
15.312 + lpx_del_cols(lp, 1, cols);
15.313 + col_iter_map[col_it]=0; //glpk specifikus
15.314 + ColIt it;
15.315 + for (col_iter_map.first(it, VALID_ID);
15.316 + col_iter_map.valid(it); col_iter_map.next(it)) {
15.317 + if (col_iter_map[it]>cols[1]) --col_iter_map[it];
15.318 + }
15.319 + }
15.320 + ///.
15.321 + void eraseRow(const RowIt& row_it) {
15.322 + row_iter_map.set(row_it, VALID_ID, INVALID_ID);
15.323 + int rows[2];
15.324 + rows[1]=row_iter_map[row_it];
15.325 + lpx_del_rows(lp, 1, rows);
15.326 + row_iter_map[row_it]=0; //glpk specifikus
15.327 + RowIt it;
15.328 + for (row_iter_map.first(it, VALID_ID);
15.329 + row_iter_map.valid(it); row_iter_map.next(it)) {
15.330 + if (row_iter_map[it]>rows[1]) --row_iter_map[it];
15.331 + }
15.332 + }
15.333 + ///.
15.334 + void setColBounds(const ColIt& col_it, int bound_type,
15.335 + double lo, double up) {
15.336 + lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
15.337 + }
15.338 + ///.
15.339 + double getObjCoef(const ColIt& col_it) {
15.340 + return lpx_get_obj_coef(lp, col_iter_map[col_it]);
15.341 + }
15.342 + ///.
15.343 + void setRowBounds(const RowIt& row_it, int bound_type,
15.344 + double lo, double up) {
15.345 + lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
15.346 + }
15.347 + ///.
15.348 + void setObjCoef(const ColIt& col_it, double obj_coef) {
15.349 + lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
15.350 + }
15.351 + ///.
15.352 + void solveSimplex() { lpx_simplex(lp); }
15.353 + ///.
15.354 + void solvePrimalSimplex() { lpx_simplex(lp); }
15.355 + ///.
15.356 + void solveDualSimplex() { lpx_simplex(lp); }
15.357 + ///.
15.358 + double getPrimal(const ColIt& col_it) {
15.359 + return lpx_get_col_prim(lp, col_iter_map[col_it]);
15.360 + }
15.361 + ///.
15.362 + double getObjVal() { return lpx_get_obj_val(lp); }
15.363 + ///.
15.364 + int rowNum() const { return lpx_get_num_rows(lp); }
15.365 + ///.
15.366 + int colNum() const { return lpx_get_num_cols(lp); }
15.367 + ///.
15.368 + int warmUp() { return lpx_warm_up(lp); }
15.369 + ///.
15.370 + void printWarmUpStatus(int i) {
15.371 + switch (i) {
15.372 + case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
15.373 + case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;
15.374 + case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
15.375 + case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
15.376 + }
15.377 + }
15.378 + ///.
15.379 + int getPrimalStatus() { return lpx_get_prim_stat(lp); }
15.380 + ///.
15.381 + void printPrimalStatus(int i) {
15.382 + switch (i) {
15.383 + case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
15.384 + case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;
15.385 + case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
15.386 + case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
15.387 + }
15.388 + }
15.389 + ///.
15.390 + int getDualStatus() { return lpx_get_dual_stat(lp); }
15.391 + ///.
15.392 + void printDualStatus(int i) {
15.393 + switch (i) {
15.394 + case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
15.395 + case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;
15.396 + case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
15.397 + case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
15.398 + }
15.399 + }
15.400 + /// Returns the status of the slack variable assigned to row \c row_it.
15.401 + int getRowStat(const RowIt& row_it) {
15.402 + return lpx_get_row_stat(lp, row_iter_map[row_it]);
15.403 + }
15.404 + ///.
15.405 + void printRowStatus(int i) {
15.406 + switch (i) {
15.407 + case LPX_BS: cout << "LPX_BS" << endl; break;
15.408 + case LPX_NL: cout << "LPX_NL" << endl; break;
15.409 + case LPX_NU: cout << "LPX_NU" << endl; break;
15.410 + case LPX_NF: cout << "LPX_NF" << endl; break;
15.411 + case LPX_NS: cout << "LPX_NS" << endl; break;
15.412 + }
15.413 + }
15.414 + /// Returns the status of the variable assigned to column \c col_it.
15.415 + int getColStat(const ColIt& col_it) {
15.416 + return lpx_get_col_stat(lp, col_iter_map[col_it]);
15.417 + }
15.418 + ///.
15.419 + void printColStatus(int i) {
15.420 + switch (i) {
15.421 + case LPX_BS: cout << "LPX_BS" << endl; break;
15.422 + case LPX_NL: cout << "LPX_NL" << endl; break;
15.423 + case LPX_NU: cout << "LPX_NU" << endl; break;
15.424 + case LPX_NF: cout << "LPX_NF" << endl; break;
15.425 + case LPX_NS: cout << "LPX_NS" << endl; break;
15.426 + }
15.427 + }
15.428 + };
15.429 +
15.430 + /// @}
15.431 +
15.432 +} //namespace lemon
15.433 +
15.434 +#endif //LEMON_LP_SOLVER_WRAPPER_H
16.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
16.2 +++ b/src/work/athos/lp_old/magic_square.cc Wed Mar 23 09:49:41 2005 +0000
16.3 @@ -0,0 +1,99 @@
16.4 +// -*- c++ -*-
16.5 +#include <iostream>
16.6 +#include <fstream>
16.7 +
16.8 +#include <lemon/time_measure.h>
16.9 +#include <lp_solver_glpk.h>
16.10 +
16.11 +using std::cout;
16.12 +using std::endl;
16.13 +using namespace lemon;
16.14 +
16.15 +/*
16.16 + On an 1537Mhz PC, the run times with
16.17 + glpk are the following.
16.18 + for n=3,4, some secondes
16.19 + for n=5, 25 hours
16.20 + */
16.21 +
16.22 +int main(int, char **) {
16.23 + const int n=3;
16.24 + const double row_sum=(1.0+n*n)*n/2;
16.25 + Timer ts;
16.26 + ts.reset();
16.27 + typedef LpGlpk LPSolver;
16.28 + typedef LPSolver::Col Col;
16.29 + LPSolver lp;
16.30 + typedef std::map<std::pair<int, int>, Col> Coords;
16.31 + Coords x;
16.32 + // we create a new variable for each entry
16.33 + // of the magic square
16.34 + for (int i=1; i<=n; ++i) {
16.35 + for (int j=1; j<=n; ++j) {
16.36 + Col col=lp.addCol();
16.37 + x[std::make_pair(i,j)]=col;
16.38 + lp.setColLowerBound(col, 1.0);
16.39 + lp.setColUpperBound(col, double(n*n));
16.40 + }
16.41 + }
16.42 + LPSolver::Expression expr3, expr4;
16.43 + for (int i=1; i<=n; ++i) {
16.44 + LPSolver::Expression expr1, expr2;
16.45 + for (int j=1; j<=n; ++j) {
16.46 + expr1+=x[std::make_pair(i, j)];
16.47 + expr2+=x[std::make_pair(j, i)];
16.48 + }
16.49 +
16.50 + // sum of rows and columns
16.51 + lp.addRow(expr1==row_sum);
16.52 + lp.addRow(expr2==row_sum);
16.53 + cout <<"Expr1: "<<expr1<<endl;
16.54 + cout <<"Expr2: "<<expr2<<endl;
16.55 +
16.56 + expr3+=x[std::make_pair(i, i)];
16.57 + expr4+=x[std::make_pair(i, (n+1)-i)];
16.58 + }
16.59 + cout <<"Expr3: "<<expr3<<endl;
16.60 + cout <<"Expr4: "<<expr4<<endl;
16.61 + // sum of the diagonal entries
16.62 + lp.addRow(expr3==row_sum);
16.63 + lp.addRow(expr4==row_sum);
16.64 + lp.solveSimplex();
16.65 + cout << "elapsed time: " << ts << endl;
16.66 + for (int i=1; i<=n; ++i) {
16.67 + for (int j=1; j<=n; ++j) {
16.68 + cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)])
16.69 + << endl;
16.70 + }
16.71 + }
16.72 +
16.73 +
16.74 +
16.75 +// // we make new binary variables for each pair of
16.76 +// // entries of the square to achieve that
16.77 +// // the values of different entries are different
16.78 +// lp.setMIP();
16.79 +// for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) {
16.80 +// Coords::const_iterator jt=it; ++jt;
16.81 +// for(; jt!=x.end(); ++jt) {
16.82 +// Col col1=(*it).second;
16.83 +// Col col2=(*jt).second;
16.84 +// Col col=lp.addCol();
16.85 +// lp.setColLowerBound(col, 0.0);
16.86 +// lp.setColUpperBound(col, 1.0);
16.87 +// lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0);
16.88 +// lp.setColInt(col);
16.89 +// }
16.90 +// }
16.91 +// cout << "elapsed time: " << ts << endl;
16.92 +// lp.solveSimplex();
16.93 +// // let's solve the integer problem
16.94 +// lp.solveBandB();
16.95 +// cout << "elapsed time: " << ts << endl;
16.96 +// for (int i=1; i<=n; ++i) {
16.97 +// for (int j=1; j<=n; ++j) {
16.98 +// cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)])
16.99 +// << endl;
16.100 +// }
16.101 +// }
16.102 +}
17.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
17.2 +++ b/src/work/athos/lp_old/makefile Wed Mar 23 09:49:41 2005 +0000
17.3 @@ -0,0 +1,73 @@
17.4 +#INCLUDEDIRS ?= -I.. -I. -I./{marci,jacint,alpar,klao,akos}
17.5 +#GLPKROOT = /usr/local/glpk-4.4
17.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
17.7 +#INCLUDEDIRS ?= -I../.. -I../.. -I../../{marci,jacint,alpar,klao,akos} -I/usr/local/glpk-4.4/include
17.8 +CXXFLAGS = -g -O2 -W -Wall $(INCLUDEDIRS) -ansi -pedantic
17.9 +LDFLAGS = -lglpk#-lcplex -lm -lpthread -lilocplex -L/usr/local/cplex/cplex75/lib/i86_linux2_glibc2.2_gcc3.0/static_mt# -L$(GLPKROOT)/lib
17.10 +
17.11 +BINARIES = magic_square max_flow_expression #expression_test max_flow_by_lp# sample sample2 sample11 sample15
17.12 +
17.13 +#include ../makefile
17.14 +
17.15 +# Hat, ez elismerem, hogy nagyon ronda, de mukodik minden altalam
17.16 +# ismert rendszeren :-) (Misi)
17.17 +ifdef GCCVER
17.18 +CXX := g++-$(GCCVER)
17.19 +else
17.20 +CXX := $(shell type -p g++-3.3 || type -p g++-3.2 || type -p g++-3.0 || type -p g++-3 || echo g++)
17.21 +endif
17.22 +
17.23 +ifdef DEBUG
17.24 +CXXFLAGS += -DDEBUG
17.25 +endif
17.26 +
17.27 +CC := $(CXX)
17.28 +
17.29 +
17.30 +all: $(BINARIES)
17.31 +
17.32 +################
17.33 +# Minden binarishoz egy sor, felsorolva, hogy mely object file-okbol
17.34 +# all elo.
17.35 +# Kiveve ha siman file.cc -> file esetrol van szo, amikor is nem kell
17.36 +# irni semmit.
17.37 +
17.38 +#proba: proba.o seged.o
17.39 +
17.40 +################
17.41 +
17.42 +
17.43 +# .depend dep depend:
17.44 +# -$(CXX) $(CXXFLAGS) -M $(BINARIES:=.cc) > .depend
17.45 +
17.46 +#makefile: .depend
17.47 +#sinclude .depend
17.48 +#moert nem megy az eredeti /usr/bin/ld-vel?
17.49 +
17.50 +# %: %.o
17.51 +# $(CXX) -o $@ $< $(LDFLAGS)
17.52 +
17.53 +# %.o: %.cc
17.54 +# $(CXX) $(CXXFLAGS) -c $<
17.55 +
17.56 +%: %.cc
17.57 + $(CXX) $(CXXFLAGS) -o $@ $< $(LDFLAGS)
17.58 +
17.59 +sample11prof: sample11prof.o
17.60 + $(CXX) -pg -o sample11prof sample11prof.o -L$(GLPKROOT)/lib -lglpk
17.61 +sample11prof.o: sample11.cc
17.62 + $(CXX) -pg $(CXXFLAGS) -c -o sample11prof.o sample11.cc
17.63 +
17.64 +# sample.o: sample.cc
17.65 +# $(CXX) $(CXXFLAGS) -c -o sample.o sample.cc
17.66 +
17.67 +# sample2: sample2.o
17.68 +# $(CXX) -o sample2 sample2.o -L/usr/local/glpk-4.4/lib -lglpk
17.69 +# sample2.o: sample2.cc
17.70 +# $(CXX) $(CXXFLAGS) -c -o sample2.o sample2.cc
17.71 +
17.72 +
17.73 +clean:
17.74 + $(RM) *.o $(BINARIES) .depend
17.75 +
17.76 +.PHONY: all clean dep depend
18.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
18.2 +++ b/src/work/athos/lp_old/max_flow_by_lp.cc Wed Mar 23 09:49:41 2005 +0000
18.3 @@ -0,0 +1,186 @@
18.4 +// -*- c++ -*-
18.5 +#include <iostream>
18.6 +#include <fstream>
18.7 +
18.8 +#include <lemon/smart_graph.h>
18.9 +#include <lemon/list_graph.h>
18.10 +#include <lemon/dimacs.h>
18.11 +#include <lemon/time_measure.h>
18.12 +//#include <graph_wrapper.h>
18.13 +#include <lemon/preflow.h>
18.14 +#include <augmenting_flow.h>
18.15 +//#include <preflow_res.h>
18.16 +//#include <lp_solver_wrapper_2.h>
18.17 +#include <min_cost_gen_flow.h>
18.18 +
18.19 +// Use a DIMACS max flow file as stdin.
18.20 +// max_flow_demo < dimacs_max_flow_file
18.21 +
18.22 +using namespace lemon;
18.23 +
18.24 +int main(int, char **) {
18.25 +
18.26 + typedef ListGraph MutableGraph;
18.27 + typedef ListGraph Graph;
18.28 + typedef Graph::Node Node;
18.29 + typedef Graph::Edge Edge;
18.30 + typedef Graph::EdgeIt EdgeIt;
18.31 +
18.32 + Graph g;
18.33 +
18.34 + Node s, t;
18.35 + Graph::EdgeMap<int> cap(g);
18.36 + //readDimacsMaxFlow(std::cin, g, s, t, cap);
18.37 + readDimacs(std::cin, g, cap, s, t);
18.38 + Timer ts;
18.39 + Graph::EdgeMap<int> flow(g); //0 flow
18.40 + Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
18.41 + max_flow_test(g, s, t, cap, flow);
18.42 + AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
18.43 + augmenting_flow_test(g, s, t, cap, flow);
18.44 +
18.45 + Graph::NodeMap<bool> cut(g);
18.46 +
18.47 + {
18.48 + std::cout << "preflow ..." << std::endl;
18.49 + ts.reset();
18.50 + max_flow_test.run();
18.51 + std::cout << "elapsed time: " << ts << std::endl;
18.52 + std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
18.53 + max_flow_test.minCut(cut);
18.54 +
18.55 + for (EdgeIt e(g); e!=INVALID; ++e) {
18.56 + if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
18.57 + std::cout << "Slackness does not hold!" << std::endl;
18.58 + if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
18.59 + std::cout << "Slackness does not hold!" << std::endl;
18.60 + }
18.61 + }
18.62 +
18.63 +// {
18.64 +// std::cout << "preflow ..." << std::endl;
18.65 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
18.66 +// ts.reset();
18.67 +// max_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
18.68 +// std::cout << "elapsed time: " << ts << std::endl;
18.69 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
18.70 +
18.71 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
18.72 +// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
18.73 +// std::cout << "Slackness does not hold!" << std::endl;
18.74 +// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
18.75 +// std::cout << "Slackness does not hold!" << std::endl;
18.76 +// }
18.77 +// }
18.78 +
18.79 +// {
18.80 +// std::cout << "wrapped preflow ..." << std::endl;
18.81 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
18.82 +// ts.reset();
18.83 +// pre_flow_res.run();
18.84 +// std::cout << "elapsed time: " << ts << std::endl;
18.85 +// std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
18.86 +// }
18.87 +
18.88 + {
18.89 + std::cout << "physical blocking flow augmentation ..." << std::endl;
18.90 + for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
18.91 + ts.reset();
18.92 + int i=0;
18.93 + while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
18.94 + std::cout << "elapsed time: " << ts << std::endl;
18.95 + std::cout << "number of augmentation phases: " << i << std::endl;
18.96 + std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
18.97 +
18.98 + for (EdgeIt e(g); e!=INVALID; ++e) {
18.99 + if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
18.100 + std::cout << "Slackness does not hold!" << std::endl;
18.101 + if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
18.102 + std::cout << "Slackness does not hold!" << std::endl;
18.103 + }
18.104 + }
18.105 +
18.106 +// {
18.107 +// std::cout << "faster physical blocking flow augmentation ..." << std::endl;
18.108 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
18.109 +// ts.reset();
18.110 +// int i=0;
18.111 +// while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
18.112 +// std::cout << "elapsed time: " << ts << std::endl;
18.113 +// std::cout << "number of augmentation phases: " << i << std::endl;
18.114 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
18.115 +// }
18.116 +
18.117 + {
18.118 + std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
18.119 + for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
18.120 + ts.reset();
18.121 + int i=0;
18.122 + while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
18.123 + std::cout << "elapsed time: " << ts << std::endl;
18.124 + std::cout << "number of augmentation phases: " << i << std::endl;
18.125 + std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
18.126 +
18.127 + for (EdgeIt e(g); e!=INVALID; ++e) {
18.128 + if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
18.129 + std::cout << "Slackness does not hold!" << std::endl;
18.130 + if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
18.131 + std::cout << "Slackness does not hold!" << std::endl;
18.132 + }
18.133 + }
18.134 +
18.135 +// {
18.136 +// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
18.137 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
18.138 +// ts.reset();
18.139 +// int i=0;
18.140 +// while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
18.141 +// std::cout << "elapsed time: " << ts << std::endl;
18.142 +// std::cout << "number of augmentation phases: " << i << std::endl;
18.143 +// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
18.144 +
18.145 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
18.146 +// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
18.147 +// std::cout << "Slackness does not hold!" << std::endl;
18.148 +// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
18.149 +// std::cout << "Slackness does not hold!" << std::endl;
18.150 +// }
18.151 +// }
18.152 +
18.153 +// {
18.154 +// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
18.155 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
18.156 +// ts.reset();
18.157 +// int i=0;
18.158 +// while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
18.159 +// std::cout << "elapsed time: " << ts << std::endl;
18.160 +// std::cout << "number of augmentation phases: " << i << std::endl;
18.161 +// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
18.162 +
18.163 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
18.164 +// if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
18.165 +// std::cout << "Slackness does not hold!" << std::endl;
18.166 +// if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
18.167 +// std::cout << "Slackness does not hold!" << std::endl;
18.168 +// }
18.169 +// }
18.170 +
18.171 + ts.reset();
18.172 +
18.173 + Edge e=g.addEdge(t, s);
18.174 + Graph::EdgeMap<int> cost(g, 0);
18.175 + cost.set(e, -1);
18.176 + cap.set(e, 10000);
18.177 + typedef ConstMap<Node, int> Excess;
18.178 + Excess excess(0);
18.179 + typedef ConstMap<Edge, int> LCap;
18.180 + LCap lcap(0);
18.181 +
18.182 + MinCostGenFlow<Graph, int, Excess, LCap>
18.183 + min_cost(g, excess, lcap, cap, flow, cost);
18.184 + min_cost.feasible();
18.185 + min_cost.runByLP();
18.186 +
18.187 + std::cout << "elapsed time: " << ts << std::endl;
18.188 + std::cout << "flow value: "<< flow[e] << std::endl;
18.189 +}
19.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
19.2 +++ b/src/work/athos/lp_old/max_flow_expression.cc Wed Mar 23 09:49:41 2005 +0000
19.3 @@ -0,0 +1,109 @@
19.4 +// -*- c++ -*-
19.5 +#include <iostream>
19.6 +#include <fstream>
19.7 +
19.8 +#include <lemon/graph_utils.h>
19.9 +#include <lemon/smart_graph.h>
19.10 +#include <lemon/list_graph.h>
19.11 +#include <lemon/dimacs.h>
19.12 +#include <lemon/time_measure.h>
19.13 +#include <lp_solver_glpk.h>
19.14 +
19.15 +using std::cout;
19.16 +using std::endl;
19.17 +using namespace lemon;
19.18 +
19.19 +template<typename Edge, typename EdgeIndexMap>
19.20 +class PrimalMap {
19.21 +protected:
19.22 + LpGlpk* lp;
19.23 + EdgeIndexMap* edge_index_map;
19.24 +public:
19.25 + PrimalMap(LpGlpk& _lp, EdgeIndexMap& _edge_index_map) :
19.26 + lp(&_lp), edge_index_map(&_edge_index_map) { }
19.27 + double operator[](Edge e) const {
19.28 + return lp->getPrimal((*edge_index_map)[e]);
19.29 + }
19.30 +};
19.31 +
19.32 +// Use a DIMACS max flow file as stdin.
19.33 +// max_flow_expresion < dimacs_max_flow_file
19.34 +
19.35 +int main(int, char **) {
19.36 +
19.37 + typedef ListGraph Graph;
19.38 + typedef Graph::Node Node;
19.39 + typedef Graph::Edge Edge;
19.40 + typedef Graph::EdgeIt EdgeIt;
19.41 +
19.42 + Graph g;
19.43 +
19.44 + Node s, t;
19.45 + Graph::EdgeMap<int> cap(g);
19.46 + readDimacs(std::cin, g, cap, s, t);
19.47 + Timer ts;
19.48 +
19.49 + typedef LpGlpk LPSolver;
19.50 + LPSolver lp;
19.51 + lp.setMaximize();
19.52 + typedef LPSolver::Col Col;
19.53 + typedef LPSolver::Row Row;
19.54 + typedef Graph::EdgeMap<Col> EdgeIndexMap;
19.55 + typedef Graph::NodeMap<Row> NodeIndexMap;
19.56 + EdgeIndexMap edge_index_map(g);
19.57 + NodeIndexMap node_index_map(g);
19.58 + PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
19.59 +
19.60 + // nonnegativity of flow and capacity function
19.61 + for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
19.62 + Col col=lp.addCol();
19.63 + edge_index_map.set(e, col);
19.64 + // interesting property in GLPK:
19.65 + // if you change the order of the following two lines, the
19.66 + // two runs of GLPK are extremely different
19.67 + lp.setColLowerBound(col, 0);
19.68 + lp.setColUpperBound(col, cap[e]);
19.69 + }
19.70 +
19.71 + for (Graph::NodeIt n(g); n!=INVALID; ++n) {
19.72 + LPSolver::Expression expr;
19.73 + for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
19.74 + expr+=edge_index_map[e];
19.75 + for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
19.76 + expr-=edge_index_map[e];
19.77 + // cost function
19.78 + if (n==s) {
19.79 + lp.setObjCoeffs(expr);
19.80 + }
19.81 + // flow conservation constraints
19.82 + if ((n!=s) && (n!=t)) {
19.83 + node_index_map.set(n, lp.addRow(expr == 0.0));
19.84 + }
19.85 + }
19.86 + lp.solveSimplex();
19.87 + cout << "elapsed time: " << ts << endl;
19.88 +// cout << "rows:" << endl;
19.89 +// for (
19.90 +// LPSolver::Rows::ClassIt i(lp.row_iter_map, 0);
19.91 +// i!=INVALID;
19.92 +// ++i) {
19.93 +// cout << i << " ";
19.94 +// }
19.95 +// cout << endl;
19.96 +// cout << "cols:" << endl;
19.97 +// for (
19.98 +// LPSolver::Cols::ClassIt i(lp.col_iter_map, 0);
19.99 +// i!=INVALID;
19.100 +// ++i) {
19.101 +// cout << i << " ";
19.102 +// }
19.103 +// cout << endl;
19.104 + lp.setMIP();
19.105 + cout << "elapsed time: " << ts << endl;
19.106 + for (LPSolver::Cols::ClassIt it(lp.col_iter_map ,1); it!=INVALID; ++it) {
19.107 + lp.setColInt(it);
19.108 + }
19.109 + cout << "elapsed time: " << ts << endl;
19.110 + lp.solveBandB();
19.111 + cout << "elapsed time: " << ts << endl;
19.112 +}
20.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
20.2 +++ b/src/work/athos/lp_old/min_cost_gen_flow.h Wed Mar 23 09:49:41 2005 +0000
20.3 @@ -0,0 +1,268 @@
20.4 +// -*- c++ -*-
20.5 +#ifndef LEMON_MIN_COST_GEN_FLOW_H
20.6 +#define LEMON_MIN_COST_GEN_FLOW_H
20.7 +#include <iostream>
20.8 +//#include <fstream>
20.9 +
20.10 +#include <lemon/smart_graph.h>
20.11 +#include <lemon/list_graph.h>
20.12 +//#include <lemon/dimacs.h>
20.13 +//#include <lemon/time_measure.h>
20.14 +//#include <graph_wrapper.h>
20.15 +#include <lemon/preflow.h>
20.16 +#include <lemon/min_cost_flow.h>
20.17 +//#include <augmenting_flow.h>
20.18 +//#include <preflow_res.h>
20.19 +#include <work/marci/merge_node_graph_wrapper.h>
20.20 +#include <work/marci/lp/lp_solver_wrapper_3.h>
20.21 +
20.22 +namespace lemon {
20.23 +
20.24 + template<typename Edge, typename EdgeIndexMap>
20.25 + class PrimalMap {
20.26 + protected:
20.27 + LPGLPK* lp;
20.28 + EdgeIndexMap* edge_index_map;
20.29 + public:
20.30 + PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) :
20.31 + lp(&_lp), edge_index_map(&_edge_index_map) { }
20.32 + double operator[](Edge e) const {
20.33 + return lp->getPrimal((*edge_index_map)[e]);
20.34 + }
20.35 + };
20.36 +
20.37 + // excess: rho-delta egyelore csak =0-ra.
20.38 + template <typename Graph, typename Num,
20.39 + typename Excess=typename Graph::template NodeMap<Num>,
20.40 + typename LCapMap=typename Graph::template EdgeMap<Num>,
20.41 + typename CapMap=typename Graph::template EdgeMap<Num>,
20.42 + typename FlowMap=typename Graph::template EdgeMap<Num>,
20.43 + typename CostMap=typename Graph::template EdgeMap<Num> >
20.44 + class MinCostGenFlow {
20.45 + protected:
20.46 + const Graph& g;
20.47 + const Excess& excess;
20.48 + const LCapMap& lcapacity;
20.49 + const CapMap& capacity;
20.50 + FlowMap& flow;
20.51 + const CostMap& cost;
20.52 + public:
20.53 + MinCostGenFlow(const Graph& _g, const Excess& _excess,
20.54 + const LCapMap& _lcapacity, const CapMap& _capacity,
20.55 + FlowMap& _flow,
20.56 + const CostMap& _cost) :
20.57 + g(_g), excess(_excess), lcapacity(_lcapacity),
20.58 + capacity(_capacity), flow(_flow), cost(_cost) { }
20.59 + bool feasible() {
20.60 + // std::cout << "making new vertices..." << std::endl;
20.61 + typedef ListGraph Graph2;
20.62 + Graph2 g2;
20.63 + typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
20.64 + // std::cout << "merging..." << std::endl;
20.65 + GW gw(g, g2);
20.66 + typename GW::Node s(INVALID, g2.addNode(), true);
20.67 + typename GW::Node t(INVALID, g2.addNode(), true);
20.68 + typedef SmartGraph Graph3;
20.69 + // std::cout << "making extender graph..." << std::endl;
20.70 + typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
20.71 +// {
20.72 +// checkConcept<StaticGraph, GWW>();
20.73 +// }
20.74 + GWW gww(gw);
20.75 + typedef AugmentingGraphWrapper<GW, GWW> GWWW;
20.76 + GWWW gwww(gw, gww);
20.77 +
20.78 + // std::cout << "making new edges..." << std::endl;
20.79 + typename GWWW::template EdgeMap<Num> translated_cap(gwww);
20.80 +
20.81 + for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
20.82 + translated_cap.set(typename GWWW::Edge(e,INVALID,false),
20.83 + capacity[e]-lcapacity[e]);
20.84 + // cout << "t_cap " << gw.id(e) << " "
20.85 + // << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
20.86 + }
20.87 +
20.88 + Num expected=0;
20.89 +
20.90 + // std::cout << "making new edges 2..." << std::endl;
20.91 + for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
20.92 + Num a=0;
20.93 + for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
20.94 + a+=lcapacity[e];
20.95 + for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
20.96 + a-=lcapacity[e];
20.97 + if (excess[n]>a) {
20.98 + typename GWW::Edge e=
20.99 + gww.addEdge(typename GW::Node(n,INVALID,false), t);
20.100 + translated_cap.set(typename GWWW::Edge(INVALID, e, true),
20.101 + excess[n]-a);
20.102 + // std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
20.103 + }
20.104 + if (excess[n]<a) {
20.105 + typename GWW::Edge e=
20.106 + gww.addEdge(s, typename GW::Node(n,INVALID,false));
20.107 + translated_cap.set(typename GWWW::Edge(INVALID, e, true),
20.108 + a-excess[n]);
20.109 + expected+=a-excess[n];
20.110 + // std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
20.111 + }
20.112 + }
20.113 +
20.114 + // std::cout << "preflow..." << std::endl;
20.115 + typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
20.116 + Preflow<GWWW, Num> preflow(gwww, s, t,
20.117 + translated_cap, translated_flow);
20.118 + preflow.run();
20.119 + // std::cout << "fv: " << preflow.flowValue() << std::endl;
20.120 + // std::cout << "expected: " << expected << std::endl;
20.121 +
20.122 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
20.123 + typename GW::Edge ew(e, INVALID, false);
20.124 + typename GWWW::Edge ewww(ew, INVALID, false);
20.125 + flow.set(e, translated_flow[ewww]+lcapacity[e]);
20.126 + }
20.127 + return (preflow.flowValue()>=expected);
20.128 + }
20.129 + // for nonnegative costs
20.130 + bool run() {
20.131 + // std::cout << "making new vertices..." << std::endl;
20.132 + typedef ListGraph Graph2;
20.133 + Graph2 g2;
20.134 + typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
20.135 + // std::cout << "merging..." << std::endl;
20.136 + GW gw(g, g2);
20.137 + typename GW::Node s(INVALID, g2.addNode(), true);
20.138 + typename GW::Node t(INVALID, g2.addNode(), true);
20.139 + typedef SmartGraph Graph3;
20.140 + // std::cout << "making extender graph..." << std::endl;
20.141 + typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
20.142 +// {
20.143 +// checkConcept<StaticGraph, GWW>();
20.144 +// }
20.145 + GWW gww(gw);
20.146 + typedef AugmentingGraphWrapper<GW, GWW> GWWW;
20.147 + GWWW gwww(gw, gww);
20.148 +
20.149 + // std::cout << "making new edges..." << std::endl;
20.150 + typename GWWW::template EdgeMap<Num> translated_cap(gwww);
20.151 +
20.152 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
20.153 + typename GW::Edge ew(e, INVALID, false);
20.154 + typename GWWW::Edge ewww(ew, INVALID, false);
20.155 + translated_cap.set(ewww, capacity[e]-lcapacity[e]);
20.156 + // cout << "t_cap " << g.id(e) << " "
20.157 + // << translated_cap[ewww] << endl;
20.158 + }
20.159 +
20.160 + Num expected=0;
20.161 +
20.162 + // std::cout << "making new edges 2..." << std::endl;
20.163 + for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
20.164 + // std::cout << "node: " << g.id(n) << std::endl;
20.165 + Num a=0;
20.166 + for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) {
20.167 + a+=lcapacity[e];
20.168 + // std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl;
20.169 + }
20.170 + for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) {
20.171 + a-=lcapacity[e];
20.172 + // std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl;
20.173 + }
20.174 + // std::cout << "excess " << g.id(n) << ": " << a << std::endl;
20.175 + if (0>a) {
20.176 + typename GWW::Edge e=
20.177 + gww.addEdge(typename GW::Node(n,INVALID,false), t);
20.178 + translated_cap.set(typename GWWW::Edge(INVALID, e, true),
20.179 + -a);
20.180 + // std::cout << g.id(n) << "->t " << -a << std::endl;
20.181 + }
20.182 + if (0<a) {
20.183 + typename GWW::Edge e=
20.184 + gww.addEdge(s, typename GW::Node(n,INVALID,false));
20.185 + translated_cap.set(typename GWWW::Edge(INVALID, e, true),
20.186 + a);
20.187 + expected+=a;
20.188 + // std::cout << "s->" << g.id(n) << " "<< a <<std:: endl;
20.189 + }
20.190 + }
20.191 +
20.192 + // std::cout << "preflow..." << std::endl;
20.193 + typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0);
20.194 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
20.195 + translated_cost.set(typename GWWW::Edge(
20.196 + typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]);
20.197 + }
20.198 + // typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
20.199 + MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>,
20.200 + typename GWWW::template EdgeMap<Num> >
20.201 + min_cost_flow(gwww, translated_cost, translated_cap,
20.202 + s, t);
20.203 + while (min_cost_flow.augment()) { }
20.204 + std::cout << "fv: " << min_cost_flow.flowValue() << std::endl;
20.205 + std::cout << "expected: " << expected << std::endl;
20.206 +
20.207 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
20.208 + typename GW::Edge ew(e, INVALID, false);
20.209 + typename GWWW::Edge ewww(ew, INVALID, false);
20.210 + // std::cout << g.id(e) << " " << flow[e] << std::endl;
20.211 + flow.set(e, lcapacity[e]+
20.212 + min_cost_flow.getFlow()[ewww]);
20.213 + }
20.214 + return (min_cost_flow.flowValue()>=expected);
20.215 + }
20.216 + void runByLP() {
20.217 + typedef LPGLPK LPSolver;
20.218 + LPSolver lp;
20.219 + lp.setMinimize();
20.220 + typedef LPSolver::ColIt ColIt;
20.221 + typedef LPSolver::RowIt RowIt;
20.222 + typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
20.223 + EdgeIndexMap edge_index_map(g);
20.224 + PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
20.225 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
20.226 + ColIt col_it=lp.addCol();
20.227 + edge_index_map.set(e, col_it);
20.228 + if (lcapacity[e]==capacity[e])
20.229 + lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]);
20.230 + else
20.231 + lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]);
20.232 + lp.setObjCoef(col_it, cost[e]);
20.233 + }
20.234 + LPSolver::ColIt col_it;
20.235 + for (lp.col_iter_map.first(col_it, lp.VALID_CLASS);
20.236 + lp.col_iter_map.valid(col_it);
20.237 + lp.col_iter_map.next(col_it)) {
20.238 +// std::cout << "ize " << lp.col_iter_map[col_it] << std::endl;
20.239 + }
20.240 + for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
20.241 + typename Graph::template EdgeMap<Num> coeffs(g, 0);
20.242 + for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
20.243 + coeffs.set(e, coeffs[e]+1);
20.244 + for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
20.245 + coeffs.set(e, coeffs[e]-1);
20.246 + RowIt row_it=lp.addRow();
20.247 + typename std::vector< std::pair<ColIt, double> > row;
20.248 + //std::cout << "node:" <<g.id(n)<<std::endl;
20.249 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
20.250 + if (coeffs[e]!=0) {
20.251 + //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
20.252 + row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
20.253 + }
20.254 + }
20.255 + //std::cout << std::endl;
20.256 + //std::cout << " " << g.id(n) << " " << row.size() << std::endl;
20.257 + lp.setRowCoeffs(row_it, row.begin(), row.end());
20.258 + lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
20.259 + }
20.260 + lp.solveSimplex();
20.261 + //std::cout << lp.colNum() << std::endl;
20.262 + //std::cout << lp.rowNum() << std::endl;
20.263 + //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
20.264 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e)
20.265 + flow.set(e, lp_flow[e]);
20.266 + }
20.267 + };
20.268 +
20.269 +} // namespace lemon
20.270 +
20.271 +#endif //LEMON_MIN_COST_GEN_FLOW_H