Add an MCF operator to the separate solver

For the solver utility it'd be good to have MCF problems, so here they
are!

Format is:
  MCF<supplies, arcs>(cost*)

Supplies is a [int,int,int,...], where each int represents a new node
Arcs is [int:int, int:int, int:int, ...] where each int:int pair
represents an edge from the first to the second (1 indexed from the
"supplies" array).

Costs is the argument to the function. There must be as many costs as
arcs, and they are set from left to right, in order.

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>