From b19bd8d8a41664328f33c9b459b2b0100f0b303f Mon Sep 17 00:00:00 2001
From: Carlo Zancanaro <carlo@carlo-laptop>
Date: Fri, 9 Nov 2012 12:17:46 +1100
Subject: 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.
---
 impl/Operator.hpp | 82 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 82 insertions(+)

(limited to 'impl/Operator.hpp')

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>
-- 
cgit v1.2.3