最小費用流問題

最小費用流問題を最短路を使って解く方法。
LPで解けるが、LPソルバを使いたくないときに使える。

• 最短路を求めて、流せるだけ流す。
• 流したら、減らせるので逆向きのアークにもなるが、コストは負になる。
• コストが負だとDijikstraが使えないので、開始ノードからのコストで補正する (補正コスト＝コスト＋元ノードのP値ー先ノードのP値)。P値は開始ノードからのコストの累積値。
• ロジックがシンプルになるので、開始用のノード(st)と終了用のノード(en)を追加した。
• 計算速度よりわかりやすさを優先。

Graphクラスは、ORToolBoxより。

using McfGraph = Graph<McfNode, McfArc>;
using McfNTag = NodeTag<McfNode, McfArc>;
using McfATag = ArcTag<McfNode, McfArc>;
/// <summary>
/// MinCostFlowProblem用
/// </summary>
public sealed class McfNode
{
	public int Request; // 要求量(正なら供給、負なら需要)
	public int Result; // 結果量(正なら供給、負なら需要)
	/// <summary>コンストラクタ</summary>
	public McfNode(int r) { Request = r; }
}
/// <summary>
/// MinCostFlowProblem用
/// </summary>
public sealed class McfArc
{
	public int Capa; // 容量
	public int Flow; // 流量(結果)
	public double Cost; // コスト
	/// <summary>コンストラクタ</summary>
	public McfArc(int cap, double cst)
	{
		Capa = cap;
		Flow = 0;
		Cost = cst;
	}
	public int Residue(bool b)
	{
		return b ? Capa - Flow : Flow;
	}
}
public static class MinCostFlowProblem
{
	/// <summary>
	/// 最小費用流問題を解く
	/// </summary>
	/// <param name="g">グラフ</param>
	/// <param name="cost">コストの総和</param>
	/// <returns>全て流せたか</returns>
	public static bool Calc(McfGraph g, out double cost)
	{
		foreach (McfATag at in g.ArcTags)
			if (!at.IsDirected) throw new ApplicationException("Must be directed");
		// ダミーノード追加
		McfNTag st = g.AddNode(new McfNode(0)); // 供給ノード
		McfNTag en = g.AddNode(new McfNode(0)); // 需要ノード
		int n = g.NodesCount; // ノード数
		// ダミーアーク追加
		foreach (McfNTag nt in g.NodeTags)
		{
			if (nt == st) break;
			nt.Node.Result = 0;
			if (nt.Node.Request > 0)
			{
				st.Node.Request += nt.Node.Request;
				g.AddArc(new McfArc(nt.Node.Request, 0), st.Index, nt.Index, true);
			}
			else if (nt.Node.Request < 0)
			{
				en.Node.Request += nt.Node.Request;
				g.AddArc(new McfArc(-nt.Node.Request, 0), nt.Index, en.Index, true);
			}
		}
		double[] p = new double[n]; // 供給点からの最短値の累積
		while (st.Node.Request > 0)
		{
			List<KeyValuePair<McfATag, bool>> minrt
				= MinRoute(g, p, st.Index, en.Index);
			if (minrt == null) break;
			// 流量計算
			int flow = int.MaxValue;
			foreach (KeyValuePair<McfATag, bool> pr in minrt)
				flow = Math.Min(flow, pr.Key.Arc.Residue(pr.Value));
			foreach (KeyValuePair<McfATag, bool> pr in minrt)
				pr.Key.Arc.Flow += pr.Value ? flow : -flow;
			st.Node.Request -= flow;
			en.Node.Request += flow;
		}
		// コスト計算
		bool bCan = st.Node.Request == 0 && en.Node.Request == 0;
		foreach (McfATag at in st.OutArcTags)
			at.RightTag.Node.Result = at.Arc.Flow;
		foreach (McfATag at in en.InArcTags)
			at.LeftTag.Node.Result = -at.Arc.Flow;
		g.DelNode(en.Index);
		g.DelNode(st.Index);
		cost = 0;
		foreach (McfATag at in g.ArcTags)
			cost += at.Arc.Cost * at.Arc.Flow;
		return bCan;
	}
	private static List<KeyValuePair<McfATag, bool>> MinRoute
		(McfGraph g, double[] p, int n1, int n2)
	{
		if (n1 == n2) return null;
		McfATag[] prev = new McfATag[g.NodesCount];
		double[] dist = new double[g.NodesCount];
		for (int k = 0; k < dist.Length; ++k) dist[k] = double.PositiveInfinity;
		dist[n1] = 0;
		bool[] found = new bool[g.NodesCount];
		int nxt = n1, nfound = 1;
		while (nfound < found.Length)
		{
			found[nxt] = true;
			++nfound;
			foreach (McfATag at in g.NodeTag(nxt).OutArcTags)
			{
				int k = at.AnotherID(nxt);
				if (found[k] || at.Arc.Residue(true) == 0) continue;
				McfNTag nt = g.NodeTag(k);
				double d = at.Arc.Cost + p[at.LeftID] - p[at.RightID];
				d += dist[nxt];
				if (d < dist[k] || (d == dist[k] && prev[k] != null))
				{
					prev[k] = at;
					dist[k] = d;
				}
			}
			foreach (McfATag at in g.NodeTag(nxt).InArcTags)
			{
				int k = at.AnotherID(nxt);
				if (found[k] || at.Arc.Residue(false) == 0) continue;
				McfNTag nt = g.NodeTag(k);
				double d = -at.Arc.Cost - p[at.LeftID] + p[at.RightID];
				d += dist[nxt];
				if (d < dist[k] || (d == dist[k] && prev[k] != null))
				{
					prev[k] = at;
					dist[k] = d;
				}
			}
			double min = double.PositiveInfinity;
			for (int i = 0; i < dist.Length; ++i)
			{
				if (!found[i] && dist[i] < min)
				{
					min = dist[i];
					nxt = i;
				}
			}
			if (min == double.PositiveInfinity) break;
		}
		if (dist[n2] == double.PositiveInfinity) return null;
		for (int i = 0; i < p.Length; ++i) p[i] += dist[i];
		List<KeyValuePair<McfATag, bool>> res
			= new List<KeyValuePair<McfATag, bool>>();
		int m = n2;
		while (m != n1)
		{
			res.Add(new KeyValuePair<McfATag, bool>
				(prev[m], prev[m].RightID == m));
			m = prev[m].AnotherID(m);
		}
		res.Reverse();
		return res;
	}
	public static void Test()
	{
		string[] ss = {
			"1,-1,-1,1/1,5,0,1;1,3,0,2;1,2,3,1;1,1,3,2/5",
			"1,0,0,-1,0/1,2,0,1;1,2,1,2;1,2,2,3;1,2.5,0,4;1,2.5,4,3/5",
			"70,0,0,0,-70/30,3,0,1;50,7,1,3;60,7,3,4;40,5,1,2;20,8,2,3;60,9,0,2;50,6,2,4/1080",
		};
		foreach (string t in ss)
		{
			McfGraph g = new McfGraph();
			string[] tt = t.Split('/');
			string[] nn = tt[0].Split(',');
			string[] aa = tt[1].Split(';');
			double cost, exp = double.Parse(tt[2]);
			foreach (string s in nn)
				g.AddNode(new McfNode(int.Parse(s)));
			foreach (string s in aa)
			{
				string[] a = s.Split(',');
				g.AddArc(new McfArc(int.Parse(a[0]), double.Parse(a[1])),
					int.Parse(a[2]), int.Parse(a[3]), true);
			}
			bool b = MinCostFlowProblem.Calc(g, out cost);
			Console.WriteLine("{0}, {1}", b, cost == exp);
		}
	}
}