1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/athos/lp/lp_solver_wrapper.h Tue Mar 22 11:45:47 2005 +0000
1.3 @@ -0,0 +1,431 @@
1.4 +// -*- c++ -*-
1.5 +#ifndef LEMON_LP_SOLVER_WRAPPER_H
1.6 +#define LEMON_LP_SOLVER_WRAPPER_H
1.7 +
1.8 +///\ingroup misc
1.9 +///\file
1.10 +///\brief Dijkstra algorithm.
1.11 +
1.12 +// #include <stdio.h>
1.13 +#include <stdlib.h>
1.14 +// #include <stdio>
1.15 +//#include <stdlib>
1.16 +extern "C" {
1.17 +#include "glpk.h"
1.18 +}
1.19 +
1.20 +#include <iostream>
1.21 +#include <vector>
1.22 +#include <string>
1.23 +#include <list>
1.24 +#include <memory>
1.25 +#include <utility>
1.26 +
1.27 +//#include <sage_graph.h>
1.28 +//#include <lemon/list_graph.h>
1.29 +//#include <lemon/graph_wrapper.h>
1.30 +#include <lemon/invalid.h>
1.31 +//#include <bfs_dfs.h>
1.32 +//#include <stp.h>
1.33 +//#include <lemon/max_flow.h>
1.34 +//#include <augmenting_flow.h>
1.35 +//#include <iter_map.h>
1.36 +
1.37 +using std::cout;
1.38 +using std::cin;
1.39 +using std::endl;
1.40 +
1.41 +namespace lemon {
1.42 +
1.43 +
1.44 + /// \addtogroup misc
1.45 + /// @{
1.46 +
1.47 + /// \brief A partitioned vector with iterable classes.
1.48 + ///
1.49 + /// This class implements a container in which the data is stored in an
1.50 + /// stl vector, the range is partitioned into sets and each set is
1.51 + /// doubly linked in a list.
1.52 + /// That is, each class is iterable by lemon iterators, and any member of
1.53 + /// the vector can bo moved to an other class.
1.54 + template <typename T>
1.55 + class IterablePartition {
1.56 + protected:
1.57 + struct Node {
1.58 + T data;
1.59 + int prev; //invalid az -1
1.60 + int next;
1.61 + };
1.62 + std::vector<Node> nodes;
1.63 + struct Tip {
1.64 + int first;
1.65 + int last;
1.66 + };
1.67 + std::vector<Tip> tips;
1.68 + public:
1.69 + /// The classes are indexed by integers from \c 0 to \c classNum()-1.
1.70 + int classNum() const { return tips.size(); }
1.71 + /// This lemon style iterator iterates through a class.
1.72 + class ClassIt;
1.73 + /// Constructor. The number of classes is to be given which is fixed
1.74 + /// over the life of the container.
1.75 + /// The partition classes are indexed from 0 to class_num-1.
1.76 + IterablePartition(int class_num) {
1.77 + for (int i=0; i<class_num; ++i) {
1.78 + Tip t;
1.79 + t.first=t.last=-1;
1.80 + tips.push_back(t);
1.81 + }
1.82 + }
1.83 + protected:
1.84 + void befuz(ClassIt it, int class_id) {
1.85 + if (tips[class_id].first==-1) {
1.86 + if (tips[class_id].last==-1) {
1.87 + nodes[it.i].prev=nodes[it.i].next=-1;
1.88 + tips[class_id].first=tips[class_id].last=it.i;
1.89 + }
1.90 + } else {
1.91 + nodes[it.i].prev=tips[class_id].last;
1.92 + nodes[it.i].next=-1;
1.93 + nodes[tips[class_id].last].next=it.i;
1.94 + tips[class_id].last=it.i;
1.95 + }
1.96 + }
1.97 + void kifuz(ClassIt it, int class_id) {
1.98 + if (tips[class_id].first==it.i) {
1.99 + if (tips[class_id].last==it.i) {
1.100 + tips[class_id].first=tips[class_id].last=-1;
1.101 + } else {
1.102 + tips[class_id].first=nodes[it.i].next;
1.103 + nodes[nodes[it.i].next].prev=-1;
1.104 + }
1.105 + } else {
1.106 + if (tips[class_id].last==it.i) {
1.107 + tips[class_id].last=nodes[it.i].prev;
1.108 + nodes[nodes[it.i].prev].next=-1;
1.109 + } else {
1.110 + nodes[nodes[it.i].next].prev=nodes[it.i].prev;
1.111 + nodes[nodes[it.i].prev].next=nodes[it.i].next;
1.112 + }
1.113 + }
1.114 + }
1.115 + public:
1.116 + /// A new element with data \c t is pushed into the vector and into class
1.117 + /// \c class_id.
1.118 + ClassIt push_back(const T& t, int class_id) {
1.119 + Node n;
1.120 + n.data=t;
1.121 + nodes.push_back(n);
1.122 + int i=nodes.size()-1;
1.123 + befuz(i, class_id);
1.124 + return i;
1.125 + }
1.126 + /// A member is moved to an other class.
1.127 + void set(ClassIt it, int old_class_id, int new_class_id) {
1.128 + kifuz(it.i, old_class_id);
1.129 + befuz(it.i, new_class_id);
1.130 + }
1.131 + /// Returns the data pointed by \c it.
1.132 + T& operator[](ClassIt it) { return nodes[it.i].data; }
1.133 + /// Returns the data pointed by \c it.
1.134 + const T& operator[](ClassIt it) const { return nodes[it.i].data; }
1.135 + ///.
1.136 + class ClassIt {
1.137 + friend class IterablePartition;
1.138 + protected:
1.139 + int i;
1.140 + public:
1.141 + /// Default constructor.
1.142 + ClassIt() { }
1.143 + /// This constructor constructs an iterator which points
1.144 + /// to the member of th container indexed by the integer _i.
1.145 + ClassIt(const int& _i) : i(_i) { }
1.146 + /// Invalid constructor.
1.147 + ClassIt(const Invalid&) : i(-1) { }
1.148 + };
1.149 + /// First member of class \c class_id.
1.150 + ClassIt& first(ClassIt& it, int class_id) const {
1.151 + it.i=tips[class_id].first;
1.152 + return it;
1.153 + }
1.154 + /// Next member.
1.155 + ClassIt& next(ClassIt& it) const {
1.156 + it.i=nodes[it.i].next;
1.157 + return it;
1.158 + }
1.159 + /// True iff the iterator is valid.
1.160 + bool valid(const ClassIt& it) const { return it.i!=-1; }
1.161 + };
1.162 +
1.163 + /// \brief Wrappers for LP solvers
1.164 + ///
1.165 + /// This class implements a lemon wrapper for glpk.
1.166 + /// Later other LP-solvers will be wrapped into lemon.
1.167 + /// The aim of this class is to give a general surface to different
1.168 + /// solvers, i.e. it makes possible to write algorithms using LP's,
1.169 + /// in which the solver can be changed to an other one easily.
1.170 + class LPSolverWrapper {
1.171 + public:
1.172 +
1.173 +// class Row {
1.174 +// protected:
1.175 +// int i;
1.176 +// public:
1.177 +// Row() { }
1.178 +// Row(const Invalid&) : i(0) { }
1.179 +// Row(const int& _i) : i(_i) { }
1.180 +// operator int() const { return i; }
1.181 +// };
1.182 +// class RowIt : public Row {
1.183 +// public:
1.184 +// RowIt(const Row& row) : Row(row) { }
1.185 +// };
1.186 +
1.187 +// class Col {
1.188 +// protected:
1.189 +// int i;
1.190 +// public:
1.191 +// Col() { }
1.192 +// Col(const Invalid&) : i(0) { }
1.193 +// Col(const int& _i) : i(_i) { }
1.194 +// operator int() const { return i; }
1.195 +// };
1.196 +// class ColIt : public Col {
1.197 +// ColIt(const Col& col) : Col(col) { }
1.198 +// };
1.199 +
1.200 + public:
1.201 + ///.
1.202 + LPX* lp;
1.203 + ///.
1.204 + typedef IterablePartition<int>::ClassIt RowIt;
1.205 + ///.
1.206 + IterablePartition<int> row_iter_map;
1.207 + ///.
1.208 + typedef IterablePartition<int>::ClassIt ColIt;
1.209 + ///.
1.210 + IterablePartition<int> col_iter_map;
1.211 + //std::vector<int> row_id_to_lp_row_id;
1.212 + //std::vector<int> col_id_to_lp_col_id;
1.213 + ///.
1.214 + const int VALID_ID;
1.215 + ///.
1.216 + const int INVALID_ID;
1.217 +
1.218 + public:
1.219 + ///.
1.220 + LPSolverWrapper() : lp(lpx_create_prob()),
1.221 + row_iter_map(2),
1.222 + col_iter_map(2),
1.223 + //row_id_to_lp_row_id(), col_id_to_lp_col_id(),
1.224 + VALID_ID(0), INVALID_ID(1) {
1.225 + lpx_set_int_parm(lp, LPX_K_DUAL, 1);
1.226 + }
1.227 + ///.
1.228 + ~LPSolverWrapper() {
1.229 + lpx_delete_prob(lp);
1.230 + }
1.231 + ///.
1.232 + void setMinimize() {
1.233 + lpx_set_obj_dir(lp, LPX_MIN);
1.234 + }
1.235 + ///.
1.236 + void setMaximize() {
1.237 + lpx_set_obj_dir(lp, LPX_MAX);
1.238 + }
1.239 + ///.
1.240 + ColIt addCol() {
1.241 + int i=lpx_add_cols(lp, 1);
1.242 + ColIt col_it;
1.243 + col_iter_map.first(col_it, INVALID_ID);
1.244 + if (col_iter_map.valid(col_it)) { //van hasznalhato hely
1.245 + col_iter_map.set(col_it, INVALID_ID, VALID_ID);
1.246 + col_iter_map[col_it]=i;
1.247 + //col_id_to_lp_col_id[col_iter_map[col_it]]=i;
1.248 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
1.249 + //col_id_to_lp_col_id.push_back(i);
1.250 + //int j=col_id_to_lp_col_id.size()-1;
1.251 + col_it=col_iter_map.push_back(i, VALID_ID);
1.252 + }
1.253 +// edge_index_map.set(e, i);
1.254 +// lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
1.255 +// lpx_set_obj_coef(lp, i, cost[e]);
1.256 + return col_it;
1.257 + }
1.258 + ///.
1.259 + RowIt addRow() {
1.260 + int i=lpx_add_rows(lp, 1);
1.261 + RowIt row_it;
1.262 + row_iter_map.first(row_it, INVALID_ID);
1.263 + if (row_iter_map.valid(row_it)) { //van hasznalhato hely
1.264 + row_iter_map.set(row_it, INVALID_ID, VALID_ID);
1.265 + row_iter_map[row_it]=i;
1.266 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
1.267 + row_it=row_iter_map.push_back(i, VALID_ID);
1.268 + }
1.269 + return row_it;
1.270 + }
1.271 + //pair<RowIt, double>-bol kell megadni egy std range-et
1.272 + ///.
1.273 + template <typename Begin, typename End>
1.274 + void setColCoeffs(const ColIt& col_it,
1.275 + Begin begin, End end) {
1.276 + int mem_length=1+lpx_get_num_rows(lp);
1.277 + int* indices = new int[mem_length];
1.278 + double* doubles = new double[mem_length];
1.279 + int length=0;
1.280 + for ( ; begin!=end; ++begin) {
1.281 + ++length;
1.282 + indices[length]=row_iter_map[begin->first];
1.283 + doubles[length]=begin->second;
1.284 + }
1.285 + lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
1.286 + delete [] indices;
1.287 + delete [] doubles;
1.288 + }
1.289 + //pair<ColIt, double>-bol kell megadni egy std range-et
1.290 + ///.
1.291 + template <typename Begin, typename End>
1.292 + void setRowCoeffs(const RowIt& row_it,
1.293 + Begin begin, End end) {
1.294 + int mem_length=1+lpx_get_num_cols(lp);
1.295 + int* indices = new int[mem_length];
1.296 + double* doubles = new double[mem_length];
1.297 + int length=0;
1.298 + for ( ; begin!=end; ++begin) {
1.299 + ++length;
1.300 + indices[length]=col_iter_map[begin->first];
1.301 + doubles[length]=begin->second;
1.302 + }
1.303 + lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
1.304 + delete [] indices;
1.305 + delete [] doubles;
1.306 + }
1.307 + ///.
1.308 + void eraseCol(const ColIt& col_it) {
1.309 + col_iter_map.set(col_it, VALID_ID, INVALID_ID);
1.310 + int cols[2];
1.311 + cols[1]=col_iter_map[col_it];
1.312 + lpx_del_cols(lp, 1, cols);
1.313 + col_iter_map[col_it]=0; //glpk specifikus
1.314 + ColIt it;
1.315 + for (col_iter_map.first(it, VALID_ID);
1.316 + col_iter_map.valid(it); col_iter_map.next(it)) {
1.317 + if (col_iter_map[it]>cols[1]) --col_iter_map[it];
1.318 + }
1.319 + }
1.320 + ///.
1.321 + void eraseRow(const RowIt& row_it) {
1.322 + row_iter_map.set(row_it, VALID_ID, INVALID_ID);
1.323 + int rows[2];
1.324 + rows[1]=row_iter_map[row_it];
1.325 + lpx_del_rows(lp, 1, rows);
1.326 + row_iter_map[row_it]=0; //glpk specifikus
1.327 + RowIt it;
1.328 + for (row_iter_map.first(it, VALID_ID);
1.329 + row_iter_map.valid(it); row_iter_map.next(it)) {
1.330 + if (row_iter_map[it]>rows[1]) --row_iter_map[it];
1.331 + }
1.332 + }
1.333 + ///.
1.334 + void setColBounds(const ColIt& col_it, int bound_type,
1.335 + double lo, double up) {
1.336 + lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
1.337 + }
1.338 + ///.
1.339 + double getObjCoef(const ColIt& col_it) {
1.340 + return lpx_get_obj_coef(lp, col_iter_map[col_it]);
1.341 + }
1.342 + ///.
1.343 + void setRowBounds(const RowIt& row_it, int bound_type,
1.344 + double lo, double up) {
1.345 + lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
1.346 + }
1.347 + ///.
1.348 + void setObjCoef(const ColIt& col_it, double obj_coef) {
1.349 + lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
1.350 + }
1.351 + ///.
1.352 + void solveSimplex() { lpx_simplex(lp); }
1.353 + ///.
1.354 + void solvePrimalSimplex() { lpx_simplex(lp); }
1.355 + ///.
1.356 + void solveDualSimplex() { lpx_simplex(lp); }
1.357 + ///.
1.358 + double getPrimal(const ColIt& col_it) {
1.359 + return lpx_get_col_prim(lp, col_iter_map[col_it]);
1.360 + }
1.361 + ///.
1.362 + double getObjVal() { return lpx_get_obj_val(lp); }
1.363 + ///.
1.364 + int rowNum() const { return lpx_get_num_rows(lp); }
1.365 + ///.
1.366 + int colNum() const { return lpx_get_num_cols(lp); }
1.367 + ///.
1.368 + int warmUp() { return lpx_warm_up(lp); }
1.369 + ///.
1.370 + void printWarmUpStatus(int i) {
1.371 + switch (i) {
1.372 + case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
1.373 + case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;
1.374 + case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
1.375 + case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
1.376 + }
1.377 + }
1.378 + ///.
1.379 + int getPrimalStatus() { return lpx_get_prim_stat(lp); }
1.380 + ///.
1.381 + void printPrimalStatus(int i) {
1.382 + switch (i) {
1.383 + case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
1.384 + case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;
1.385 + case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
1.386 + case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
1.387 + }
1.388 + }
1.389 + ///.
1.390 + int getDualStatus() { return lpx_get_dual_stat(lp); }
1.391 + ///.
1.392 + void printDualStatus(int i) {
1.393 + switch (i) {
1.394 + case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
1.395 + case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;
1.396 + case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
1.397 + case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
1.398 + }
1.399 + }
1.400 + /// Returns the status of the slack variable assigned to row \c row_it.
1.401 + int getRowStat(const RowIt& row_it) {
1.402 + return lpx_get_row_stat(lp, row_iter_map[row_it]);
1.403 + }
1.404 + ///.
1.405 + void printRowStatus(int i) {
1.406 + switch (i) {
1.407 + case LPX_BS: cout << "LPX_BS" << endl; break;
1.408 + case LPX_NL: cout << "LPX_NL" << endl; break;
1.409 + case LPX_NU: cout << "LPX_NU" << endl; break;
1.410 + case LPX_NF: cout << "LPX_NF" << endl; break;
1.411 + case LPX_NS: cout << "LPX_NS" << endl; break;
1.412 + }
1.413 + }
1.414 + /// Returns the status of the variable assigned to column \c col_it.
1.415 + int getColStat(const ColIt& col_it) {
1.416 + return lpx_get_col_stat(lp, col_iter_map[col_it]);
1.417 + }
1.418 + ///.
1.419 + void printColStatus(int i) {
1.420 + switch (i) {
1.421 + case LPX_BS: cout << "LPX_BS" << endl; break;
1.422 + case LPX_NL: cout << "LPX_NL" << endl; break;
1.423 + case LPX_NU: cout << "LPX_NU" << endl; break;
1.424 + case LPX_NF: cout << "LPX_NF" << endl; break;
1.425 + case LPX_NS: cout << "LPX_NS" << endl; break;
1.426 + }
1.427 + }
1.428 + };
1.429 +
1.430 + /// @}
1.431 +
1.432 +} //namespace lemon
1.433 +
1.434 +#endif //LEMON_LP_SOLVER_WRAPPER_H