1 files changed, 82 insertions, 0 deletions
diff --git a/impl/Operator.hpp b/impl/Operator.hpp
index 64ef096..4a573b8 100644
--- a/impl/Operator.hpp
+++ b/impl/Operator.hpp
@@ -135,6 +135,88 @@ struct Guard : public Operator<Domain> {
   }
 };
 
+#include "MCFSimplex.h"
+
+template<typename Domain>
+struct MinCostFlow : public Operator<Domain> {
+  MinCostFlow(const std::vector<Domain>& supplies, const std::vector<std::pair<int,int> > arcs)
+    : _supplies(supplies), _arcs(arcs), _solver(0,0) {
+    MCFClass::FNumber* deficits = new MCFClass::FNumber[supplies.size()];
+    MCFClass::Index* starts = new MCFClass::Index[arcs.size()];
+    MCFClass::Index* ends = new MCFClass::Index[arcs.size()];
+
+    for (int i = 0, size = supplies.size(); i < size; ++i) {
+      deficits[i] = -supplies[i].template as<MCFClass::FNumber>();
+    }
+    for (int i = 0, size = arcs.size(); i < size; ++i) {
+      starts[i] = arcs[i].first;
+      ends[i] = arcs[i].second;
+    }
+
+    _solver.LoadNet(supplies.size(), supplies.size(), // max nodes/arcs
+                    supplies.size(), supplies.size(), // current nodes/arcs
+                    NULL, NULL, // arcs have inf cap, zero cost (to begin)
+                    deficits, // deficits for each node
+                    starts, ends); // start/end of each edge
+
+    delete[] deficits;
+    delete[] starts;
+    delete[] ends;
+  }
+  Domain eval (const std::vector<Domain>& costs) const {
+    assert(costs.size() == _arcs.size());
+    if (costs.size() == 0)
+      return 0;
+    for (int i = 0, size = costs.size(); i < size; ++i) {
+      _solver.ChgCost(i, costs[i].template as<MCFClass::CNumber>());
+    }
+    _solver.SolveMCF();
+    if (_solver.MCFGetStatus() == MCFClass::kUnfeasible){
+      return -infinity<Domain>();
+    } else if (_solver.MCFGetFO() == MCFClass::Inf<MCFClass::FONumber>()) {
+      return infinity<Domain>();
+    } else if (_solver.MCFGetFO() == -MCFClass::Inf<MCFClass::FONumber>()) {
+      return -infinity<Domain>();        
+    } else {
+      return _solver.MCFGetFO();
+    }
+  }
+  void print(std::ostream& cout) const {
+    std::string supplyString = "[";
+    for (int i = 0, size = _supplies.size(); i < size; ++i) {
+      if (i > 0)
+        supplyString += ",";
+      std::stringstream stream;
+      stream << _supplies[i];
+      supplyString += stream.str();
+    }
+    supplyString += ']';
+
+    std::string arcString = "[";
+    for (int i = 0, size = _arcs.size(); i < size; ++i) {
+      if (i > 0)
+        arcString += ",";
+      {
+        std::stringstream stream;
+        stream << _arcs[i].first;
+        arcString += stream.str() + ":";
+      }
+      {
+        std::stringstream stream;
+        stream << _arcs[i].second;
+        arcString += stream.str();
+      }
+    }
+    arcString += ']';
+
+    cout << "MCF<" << supplyString << "," << arcString << ">";
+  }
+private:
+  std::vector<Domain> _supplies;
+  std::vector<std::pair<int,int> > _arcs;
+  mutable MCFSimplex _solver;
+};
+
 /*#include "TemplateConstraintMatrix.hpp"
 
 template<typename Domain>