src/work/marci/lp/max_flow_expression.cc
author marci
Mon, 31 Jan 2005 17:00:12 +0000
changeset 1110 ba28dfbea5f2
parent 1104 23a54f889272
child 1111 88ade201ffc6
permissions -rw-r--r--
new functions for changing lower and upper bounds of variables
marci@1104
     1
// -*- c++ -*-
marci@1104
     2
#include <iostream>
marci@1104
     3
#include <fstream>
marci@1104
     4
marci@1104
     5
#include <lemon/smart_graph.h>
marci@1104
     6
#include <lemon/list_graph.h>
marci@1104
     7
#include <lemon/dimacs.h>
marci@1104
     8
#include <lemon/time_measure.h>
marci@1104
     9
#include <lp_solver_wrapper_3.h>
marci@1104
    10
marci@1104
    11
using std::cout;
marci@1104
    12
using std::endl;
marci@1104
    13
using namespace lemon;
marci@1104
    14
marci@1104
    15
template<typename Edge, typename EdgeIndexMap> 
marci@1104
    16
class PrimalMap {
marci@1104
    17
protected:
marci@1104
    18
  LPGLPK* lp;
marci@1104
    19
  EdgeIndexMap* edge_index_map;
marci@1104
    20
public:
marci@1104
    21
  PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
marci@1104
    22
    lp(&_lp), edge_index_map(&_edge_index_map) { }
marci@1104
    23
  double operator[](Edge e) const { 
marci@1104
    24
    return lp->getPrimal((*edge_index_map)[e]);
marci@1104
    25
  }
marci@1104
    26
};
marci@1104
    27
marci@1104
    28
// Use a DIMACS max flow file as stdin.
marci@1104
    29
// max_flow_expresion < dimacs_max_flow_file
marci@1104
    30
marci@1104
    31
int main(int, char **) {
marci@1104
    32
marci@1104
    33
  typedef ListGraph Graph;
marci@1104
    34
  typedef Graph::Node Node;
marci@1104
    35
  typedef Graph::Edge Edge;
marci@1104
    36
  typedef Graph::EdgeIt EdgeIt;
marci@1104
    37
marci@1104
    38
  Graph g;
marci@1104
    39
  
marci@1104
    40
  Node s, t;
marci@1104
    41
  Graph::EdgeMap<int> cap(g);
marci@1104
    42
  readDimacs(std::cin, g, cap, s, t);
marci@1104
    43
  Timer ts;
marci@1104
    44
  
marci@1104
    45
  typedef LPGLPK LPSolver;
marci@1104
    46
  LPSolver lp;
marci@1104
    47
  lp.setMaximize();
marci@1104
    48
  typedef LPSolver::ColIt ColIt;
marci@1104
    49
  typedef LPSolver::RowIt RowIt;
marci@1104
    50
  typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
marci@1104
    51
  EdgeIndexMap edge_index_map(g);
marci@1110
    52
  PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
marci@1104
    53
marci@1104
    54
  // capacity function
marci@1104
    55
  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
marci@1104
    56
    ColIt col_it=lp.addCol();
marci@1104
    57
    edge_index_map.set(e, col_it);
marci@1110
    58
    // interesting property in GLPK:
marci@1110
    59
    // if you change the order of the following two lines, the 
marci@1110
    60
    // two runs of GLPK are extremely different
marci@1110
    61
      lp.setColUpperBound(col_it, cap[e]);
marci@1110
    62
      lp.setColLowerBound(col_it, 0);
marci@1104
    63
  }
marci@1104
    64
  
marci@1104
    65
  for (Graph::NodeIt n(g); n!=INVALID; ++n) {
marci@1104
    66
    LPSolver::Expression expr;
marci@1104
    67
    for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
marci@1104
    68
      expr+=edge_index_map[e];
marci@1104
    69
    for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
marci@1104
    70
      expr-=edge_index_map[e];
marci@1104
    71
    // cost function
marci@1104
    72
    if (n==s) {
marci@1104
    73
      lp.setObjCoeffs(expr);      
marci@1104
    74
    }
marci@1104
    75
    // flow conservation
marci@1104
    76
    if ((n!=s) && (n!=t)) {
marci@1104
    77
      RowIt row_it=lp.addRow();
marci@1104
    78
      lp.setRowCoeffs(row_it, expr);
marci@1104
    79
      lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
marci@1104
    80
    }
marci@1104
    81
  }
marci@1104
    82
  lp.solveSimplex();
marci@1110
    83
  cout << "elapsed time: " << ts << endl;
marci@1104
    84
}