Differential D5594 Diff 18116 extern/quadriflow/src/flow.hpp

Changeset View

Standalone View

extern/quadriflow/src/flow.hpp

This file was added.

				#ifndef FLOW_H_
				#define FLOW_H_

				#include <Eigen/Core>
				#include <list>
				#include <map>
				#include <vector>

				#include "config.hpp"

				#include <boost/graph/adjacency_list.hpp>
				#include <boost/graph/boykov_kolmogorov_max_flow.hpp>
				#include <boost/graph/edmonds_karp_max_flow.hpp>
				#include <boost/graph/push_relabel_max_flow.hpp>

				#include <lemon/network_simplex.h>
				#include <lemon/preflow.h>
				#include <lemon/smart_graph.h>

				using namespace boost;
				using namespace Eigen;

				namespace qflow {

				class MaxFlowHelper {
				public:
				MaxFlowHelper() {}
				virtual ~MaxFlowHelper(){};
				virtual void resize(int n, int m) = 0;
				virtual void addEdge(int x, int y, int c, int rc, int v, int cost = 1) = 0;
				virtual int compute() = 0;
				virtual void applyTo(std::vector<Vector2i>& edge_diff) = 0;
				};

				class BoykovMaxFlowHelper : public MaxFlowHelper {
				public:
				typedef int EdgeWeightType;
				typedef adjacency_list_traits<vecS, vecS, directedS> Traits;
				// clang-format off
				typedef adjacency_list < vecS, vecS, directedS,
				property < vertex_name_t, std::string,
				property < vertex_index_t, long,
				property < vertex_color_t, boost::default_color_type,
				property < vertex_distance_t, long,
				property < vertex_predecessor_t, Traits::edge_descriptor > > > > >,

				property < edge_capacity_t, EdgeWeightType,
				property < edge_residual_capacity_t, EdgeWeightType,
				property < edge_reverse_t, Traits::edge_descriptor > > > > Graph;
				// clang-format on

				public:
				BoykovMaxFlowHelper() { rev = get(edge_reverse, g); }
				void resize(int n, int m) {
				vertex_descriptors.resize(n);
				for (int i = 0; i < n; ++i) vertex_descriptors[i] = add_vertex(g);
				}
				int compute() {
				EdgeWeightType flow =
				boykov_kolmogorov_max_flow(g, vertex_descriptors.front(), vertex_descriptors.back());
				return flow;
				}
				void addDirectEdge(Traits::vertex_descriptor& v1, Traits::vertex_descriptor& v2,
				property_map<Graph, edge_reverse_t>::type& rev, const int capacity,
				const int inv_capacity, Graph& g, Traits::edge_descriptor& e1,
				Traits::edge_descriptor& e2) {
				e1 = add_edge(v1, v2, g).first;
				e2 = add_edge(v2, v1, g).first;
				put(edge_capacity, g, e1, capacity);
				put(edge_capacity, g, e2, inv_capacity);

				rev[e1] = e2;
				rev[e2] = e1;
				}
				void addEdge(int x, int y, int c, int rc, int v, int cost = 1) {
				Traits::edge_descriptor e1, e2;
				addDirectEdge(vertex_descriptors[x], vertex_descriptors[y], rev, c, rc, g, e1, e2);
				if (v != -1) {
				edge_to_variables[e1] = std::make_pair(v, -1);
				edge_to_variables[e2] = std::make_pair(v, 1);
				}
				}
				void applyTo(std::vector<Vector2i>& edge_diff) {
				property_map<Graph, edge_capacity_t>::type capacity = get(edge_capacity, g);
				property_map<Graph, edge_residual_capacity_t>::type residual_capacity =
				get(edge_residual_capacity, g);

				graph_traits<Graph>::vertex_iterator u_iter, u_end;
				graph_traits<Graph>::out_edge_iterator ei, e_end;
				for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
				for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei)
				if (capacity[*ei] > 0) {
				int flow = (capacity[ei] - residual_capacity[ei]);
				if (flow > 0) {
				auto it = edge_to_variables.find(*ei);
				if (it != edge_to_variables.end()) {
				edge_diff[it->second.first / 2][it->second.first % 2] +=
				it->second.second * flow;
				}
				}
				}
				}

				private:
				Graph g;
				property_map<Graph, edge_reverse_t>::type rev;
				std::vector<Traits::vertex_descriptor> vertex_descriptors;
				std::map<Traits::edge_descriptor, std::pair<int, int>> edge_to_variables;
				};

				class NetworkSimplexFlowHelper : public MaxFlowHelper {
				public:
				using Weight = int;
				using Capacity = int;
				using Graph = lemon::SmartDigraph;
				using Node = Graph::Node;
				using Arc = Graph::Arc;
				template <typename ValueType>
				using ArcMap = lemon::SmartDigraph::ArcMap<ValueType>;
				using Preflow = lemon::Preflow<lemon::SmartDigraph, ArcMap<Capacity>>;
				using NetworkSimplex = lemon::NetworkSimplex<lemon::SmartDigraph, Capacity, Weight>;

				public:
				NetworkSimplexFlowHelper() : cost(graph), capacity(graph), flow(graph), variable(graph) {}
				~NetworkSimplexFlowHelper(){};
				void resize(int n, int m) {
				nodes.reserve(n);
				for (int i = 0; i < n; ++i) nodes.push_back(graph.addNode());
				}
				void addEdge(int x, int y, int c, int rc, int v, int cst = 1) {
				assert(x >= 0);
				assert(v >= -1);
				if (c) {
				auto e1 = graph.addArc(nodes[x], nodes[y]);
				cost[e1] = cst;
				capacity[e1] = c;
				variable[e1] = std::make_pair(v, 1);
				}

				if (rc) {
				auto e2 = graph.addArc(nodes[y], nodes[x]);
				cost[e2] = cst;
				capacity[e2] = rc;
				variable[e2] = std::make_pair(v, -1);
				}
				}
				int compute() {
				Preflow pf(graph, capacity, nodes.front(), nodes.back());
				NetworkSimplex ns(graph);

				// Run preflow to find maximum flow
				lprintf("push-relabel flow... ");
				pf.runMinCut();
				int maxflow = pf.flowValue();

				// Run network simplex to find minimum cost maximum flow
				ns.costMap(cost).upperMap(capacity).stSupply(nodes.front(), nodes.back(), maxflow);
				auto status = ns.run();
				switch (status) {
				case NetworkSimplex::OPTIMAL:
				ns.flowMap(flow);
				break;
				case NetworkSimplex::INFEASIBLE:
				lputs("NetworkSimplex::INFEASIBLE");
				assert(0);
				break;
				default:
				lputs("Unknown: NetworkSimplex::Default");
				assert(0);
				break;
				}

				return maxflow;
				}
				void applyTo(std::vector<Vector2i>& edge_diff) {
				for (Graph::ArcIt e(graph); e != lemon::INVALID; ++e) {
				int var = variable[e].first;
				if (var == -1) continue;
				int sgn = variable[e].second;
				edge_diff[var / 2][var % 2] -= sgn * flow[e];
				}
				}

				private:
				Graph graph;
				ArcMap<Weight> cost;
				ArcMap<Capacity> capacity;
				ArcMap<Capacity> flow;
				ArcMap<std::pair<int, int>> variable;
				std::vector<Node> nodes;
				std::vector<Arc> edges;
				};

				#ifdef WITH_GUROBI

				#include <gurobi_c++.h>

				class GurobiFlowHelper : public MaxFlowHelper {
				public:
				GurobiFlowHelper() {}
				virtual ~GurobiFlowHelper(){};
				virtual void resize(int n, int m) {
				nodes.resize(n * 2);
				edges.resize(m);
				}
				virtual void addEdge(int x, int y, int c, int rc, int v, int cost = 1) {
				nodes[x * 2 + 0].push_back(vars.size());
				nodes[y * 2 + 1].push_back(vars.size());
				vars.push_back(model.addVar(0, c, 0, GRB_INTEGER));
				edges.push_back(std::make_pair(v, 1));

				nodes[y * 2 + 0].push_back(vars.size());
				nodes[x * 2 + 1].push_back(vars.size());
				vars.push_back(model.addVar(0, rc, 0, GRB_INTEGER));
				edges.push_back(std::make_pair(v, -1));
				}
				virtual int compute() {
				std::cerr << "compute" << std::endl;
				int ns = nodes.size() / 2;

				int flow;
				for (int i = 1; i < ns - 1; ++i) {
				GRBLinExpr cons = 0;
				for (auto n : nodes[2 * i + 0]) cons += vars[n];
				for (auto n : nodes[2 * i + 1]) cons -= vars[n];
				model.addConstr(cons == 0);
				}

				// first pass, maximum flow
				GRBLinExpr outbound = 0;
				{
				lprintf("first pass\n");
				for (auto& n : nodes[0]) outbound += vars[n];
				for (auto& n : nodes[1]) outbound -= vars[n];
				model.setObjective(outbound, GRB_MAXIMIZE);
				model.optimize();

				flow = (int)model.get(GRB_DoubleAttr_ObjVal);
				lprintf("Gurobi result: %d\n", flow);
				}

				// second pass, minimum cost flow
				{
				lprintf("second pass\n");
				model.addConstr(outbound == flow);
				GRBLinExpr cost = 0;
				for (auto& v : vars) cost += v;
				model.setObjective(cost, GRB_MINIMIZE);
				model.optimize();

				double optimal_cost = (int)model.get(GRB_DoubleAttr_ObjVal);
				lprintf("Gurobi result: %.3f\n", optimal_cost);
				}
				return flow;
				}
				virtual void applyTo(std::vector<Vector2i>& edge_diff) { assert(0); };

				private:
				GRBEnv env = GRBEnv();
				GRBModel model = GRBModel(env);
				std::vector<GRBVar> vars;
				std::vector<std::pair<int, int>> edges;
				std::vector<std::vector<int>> nodes;
				};

				#endif

				class ECMaxFlowHelper : public MaxFlowHelper {
				public:
				struct FlowInfo {
				int id;
				int capacity, flow;
				int v, d;
				FlowInfo* rev;
				};
				struct SearchInfo {
				SearchInfo(int _id, int _prev_id, FlowInfo* _info)
				: id(_id), prev_id(_prev_id), info(_info) {}
				int id;
				int prev_id;
				FlowInfo* info;
				};
				ECMaxFlowHelper() { num = 0; }
				int num;
				std::vector<FlowInfo*> variable_to_edge;
				void resize(int n, int m) {
				graph.resize(n);
				variable_to_edge.resize(m, 0);
				num = n;
				}
				void addEdge(int x, int y, int c, int rc, int v, int cost = 0) {
				FlowInfo flow;
				flow.id = y;
				flow.capacity = c;
				flow.flow = 0;
				flow.v = v;
				flow.d = -1;
				graph[x].push_back(flow);
				auto& f1 = graph[x].back();
				flow.id = x;
				flow.capacity = rc;
				flow.flow = 0;
				flow.v = v;
				flow.d = 1;
				graph[y].push_back(flow);
				auto& f2 = graph[y].back();
				f2.rev = &f1;
				f1.rev = &f2;
				}
				int compute() {
				int total_flow = 0;
				int count = 0;
				while (true) {
				count += 1;
				std::vector<int> vhash(num, 0);
				std::vector<SearchInfo> q;
				q.push_back(SearchInfo(0, -1, 0));
				vhash[0] = 1;
				int q_front = 0;
				bool found = false;
				while (q_front < q.size()) {
				int vert = q[q_front].id;
				for (auto& l : graph[vert]) {
				if (vhash[l.id] \|\| l.capacity <= l.flow) continue;
				q.push_back(SearchInfo(l.id, q_front, &l));
				vhash[l.id] = 1;
				if (l.id == num - 1) {
				found = true;
				break;
				}
				}
				if (found) break;
				q_front += 1;
				}
				if (q_front == q.size()) break;
				int loc = q.size() - 1;
				while (q[loc].prev_id != -1) {
				q[loc].info->flow += 1;
				q[loc].info->rev->flow -= 1;
				loc = q[loc].prev_id;
				// int prev_v = q[loc].id;
				// applyFlow(prev_v, current_v, 1);
				// applyFlow(current_v, prev_v, -1);
				}
				total_flow += 1;
				}
				return total_flow;
				}
				void applyTo(std::vector<Vector2i>& edge_diff) {
				for (int i = 0; i < graph.size(); ++i) {
				for (auto& flow : graph[i]) {
				if (flow.flow > 0 && flow.v != -1) {
				if (flow.flow > 0) {
				edge_diff[flow.v / 2][flow.v % 2] += flow.d * flow.flow;
				if (abs(edge_diff[flow.v / 2][flow.v % 2]) > 2) {
				}
				}
				}
				}
				}
				}
				void applyFlow(int v1, int v2, int flow) {
				for (auto& it : graph[v1]) {
				if (it.id == v2) {
				it.flow += flow;
				break;
				}
				}
				}
				std::vector<std::list<FlowInfo>> graph;
				};

				} // namespace qflow

				#endif