# HG changeset patch
# User marci
# Date 1106922812 0
# Node ID 23a54f88927263f550c5c94a53d909475615c907
# Parent  f196dc4f1b313b45fb8f027574915b118f972eb3
small changes, a try for max flow using expression

diff -r f196dc4f1b31 -r 23a54f889272 src/work/marci/lp/lp_solver_wrapper_3.h
--- a/src/work/marci/lp/lp_solver_wrapper_3.h	Fri Jan 28 09:09:59 2005 +0000
+++ b/src/work/marci/lp/lp_solver_wrapper_3.h	Fri Jan 28 14:33:32 2005 +0000
@@ -10,6 +10,7 @@
 #include <stdlib.h>
 #include <iostream>
 #include <map>
+#include <limits>
 // #include <stdio>
 //#include <stdlib>
 extern "C" {
@@ -196,6 +197,8 @@
     const int VALID_CLASS;
     /// \e
     const int INVALID_CLASS;
+    /// \e 
+    static const _Value INF;
   public:
     /// \e
     LPSolverBase() : row_iter_map(2), 
@@ -221,10 +224,10 @@
     virtual int _addCol() = 0;
     /// \e
     virtual void _setRowCoeffs(int i, 
-			       std::vector<std::pair<int, double> > coeffs) = 0;
+			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;
     /// \e
     virtual void _setColCoeffs(int i, 
-			       std::vector<std::pair<int, double> > coeffs) = 0;
+			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;
     /// \e
     virtual void _eraseCol(int i) = 0;
     /// \e
@@ -423,6 +426,9 @@
     virtual void printColStatus(int i) = 0;
   };
   
+  template <typename _Value>
+  const _Value LPSolverBase<_Value>::INF=std::numeric_limits<_Value>::infinity();
+
 
   /// \brief Wrappers for LP solvers
   /// 
@@ -474,7 +480,7 @@
     }
     /// \e
     virtual void _setRowCoeffs(int i, 
-			       std::vector<std::pair<int, double> > coeffs) {
+			       const std::vector<std::pair<int, double> >& coeffs) {
       int mem_length=1+colNum();
       int* indices = new int[mem_length];
       double* doubles = new double[mem_length];
@@ -494,7 +500,7 @@
     }
     /// \e
     virtual void _setColCoeffs(int i, 
-			       std::vector<std::pair<int, double> > coeffs) {
+			       const std::vector<std::pair<int, double> >& coeffs) {
       int mem_length=1+rowNum();
       int* indices = new int[mem_length];
       double* doubles = new double[mem_length];
diff -r f196dc4f1b31 -r 23a54f889272 src/work/marci/lp/makefile
--- a/src/work/marci/lp/makefile	Fri Jan 28 09:09:59 2005 +0000
+++ b/src/work/marci/lp/makefile	Fri Jan 28 14:33:32 2005 +0000
@@ -5,7 +5,7 @@
 CXXFLAGS = -g -O2 -W -Wall $(INCLUDEDIRS) -ansi -pedantic
 LDFLAGS  =  -lglpk#-lcplex -lm -lpthread -lilocplex -L/usr/local/cplex/cplex75/lib/i86_linux2_glibc2.2_gcc3.0/static_mt# -L$(GLPKROOT)/lib
 
-BINARIES = expression_test max_flow_by_lp# sample sample2 sample11 sample15
+BINARIES = max_flow_expression expression_test max_flow_by_lp# sample sample2 sample11 sample15
 
 #include ../makefile
 
diff -r f196dc4f1b31 -r 23a54f889272 src/work/marci/lp/max_flow_expression.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/marci/lp/max_flow_expression.cc	Fri Jan 28 14:33:32 2005 +0000
@@ -0,0 +1,87 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+#include <lemon/dimacs.h>
+#include <lemon/time_measure.h>
+#include <lp_solver_wrapper_3.h>
+
+using std::cout;
+using std::endl;
+using namespace lemon;
+
+template<typename Edge, typename EdgeIndexMap> 
+class PrimalMap {
+protected:
+  LPGLPK* lp;
+  EdgeIndexMap* edge_index_map;
+public:
+  PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
+    lp(&_lp), edge_index_map(&_edge_index_map) { }
+  double operator[](Edge e) const { 
+    return lp->getPrimal((*edge_index_map)[e]);
+  }
+};
+
+// Use a DIMACS max flow file as stdin.
+// max_flow_expresion < dimacs_max_flow_file
+
+int main(int, char **) {
+
+  typedef ListGraph Graph;
+  typedef Graph::Node Node;
+  typedef Graph::Edge Edge;
+  typedef Graph::EdgeIt EdgeIt;
+
+  Graph g;
+  
+  Node s, t;
+  Graph::EdgeMap<int> cap(g);
+  readDimacs(std::cin, g, cap, s, t);
+  Timer ts;
+  
+  typedef LPGLPK LPSolver;
+  LPSolver lp;
+  lp.setMaximize();
+  typedef LPSolver::ColIt ColIt;
+  typedef LPSolver::RowIt RowIt;
+  typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
+  EdgeIndexMap edge_index_map(g);
+  PrimalMap<Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
+
+  // capacity function
+  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
+    ColIt col_it=lp.addCol();
+    edge_index_map.set(e, col_it);
+    if (cap[e]==0)
+      lp.setColBounds(col_it, LPSolver::FIXED, 0, cap[e]);
+    else 
+      lp.setColBounds(col_it, LPSolver::DOUBLE, 0, cap[e]);
+  }
+  
+  for (Graph::NodeIt n(g); n!=INVALID; ++n) {
+    LPSolver::Expression expr;
+    for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
+      expr+=edge_index_map[e];
+    for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
+      expr-=edge_index_map[e];
+    // cost function
+    if (n==s) {
+      lp.setObjCoeffs(expr);      
+    }
+    // flow conservation
+    if ((n!=s) && (n!=t)) {
+      RowIt row_it=lp.addRow();
+      lp.setRowCoeffs(row_it, expr);
+      lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
+    }
+  }
+  lp.solveSimplex();
+  //std::cout << lp.colNum() << std::endl;
+  //std::cout << lp.rowNum() << std::endl;
+  //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
+  for (Graph::EdgeIt e(g); e!=INVALID; ++e) 
+    flow.set(e, lp_flow[e]);
+}