Skip to content

Commit c91556d

Browse files
kzrnmkzrnm
committed
1 parent 8947720 commit c91556d

File tree

3 files changed

+164
-132
lines changed

3 files changed

+164
-132
lines changed

‎Source/AtCoderLibrary/Graph/Internal/CSR.cs

+54
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
using System.ComponentModel;
2+
using System.Collections.Generic;
3+
4+
namespace AtCoder.Internal
5+
{
6+
/// <summary>
7+
/// 有向グラフの辺集合を表します。
8+
/// </summary>
9+
/// <example>
10+
/// <code>
11+
/// for (int i = graph.Starts[v]; i &gt; graph.Starts[v + 1]; i++)
12+
/// {
13+
/// int to = graph.Edges[i];
14+
/// }
15+
/// </code>
16+
/// </example>
17+
[EditorBrowsable(EditorBrowsableState.Never)]
18+
public class CSR<TEdge>
19+
{
20+
/// <summary>
21+
/// 各頂点から伸びる有向辺数の累積和を取得します。
22+
/// </summary>
23+
public readonly int[] Start;
24+
25+
/// <summary>
26+
/// 有向辺の終点の配列を取得します。
27+
/// </summary>
28+
public readonly TEdge[] EList;
29+
30+
public CSR(int n, List<(int from, TEdge e)> edges)
31+
{
32+
// 本家 C++ 版 ACL を参考に実装。通常の隣接リストと比較して高速か否かは未検証。
33+
Start = new int[n + 1];
34+
EList = new TEdge[edges.Count];
35+
36+
foreach (var (from, _) in edges)
37+
{
38+
Start[from + 1]++;
39+
}
40+
41+
for (int i = 1; i <= n; i++)
42+
{
43+
Start[i] += Start[i - 1];
44+
}
45+
46+
var counter = (int[])Start.Clone();
47+
foreach (var (from, e) in edges)
48+
{
49+
EList[counter[from]++] = e;
50+
}
51+
}
52+
}
53+
54+
}

‎Source/AtCoderLibrary/Graph/MinCostFlow.cs

+103-75
Original file line numberDiff line numberDiff line change
@@ -59,12 +59,6 @@ public class McfGraph<TCap, TCapOp, TCost, TCostOp, TCast>
5959
public McfGraph(int n)
6060
{
6161
_n = n;
62-
_g = new List<EdgeInternal>[n];
63-
for (int i = 0; i < n; i++)
64-
{
65-
_g[i] = new List<EdgeInternal>();
66-
}
67-
_pos = new List<(int first, int second)>();
6862
}
6963

7064
/// <summary>
@@ -90,16 +84,8 @@ public int AddEdge(int from, int to, TCap cap, TCost cost)
9084
DebugUtil.Assert(0 <= to && to < _n);
9185
DebugUtil.Assert(capOp.LessThanOrEqual(default, cap));
9286
DebugUtil.Assert(costOp.LessThanOrEqual(default, cost));
93-
int m = _pos.Count;
94-
_pos.Add((from, _g[from].Count));
95-
96-
int fromId = _g[from].Count;
97-
int toId = _g[to].Count;
98-
99-
if (from == to) toId++;
100-
101-
_g[from].Add(new EdgeInternal(to, toId, cap, cost));
102-
_g[to].Add(new EdgeInternal(from, fromId, default, costOp.Minus(cost)));
87+
int m = _edges.Count;
88+
_edges.Add(new Edge(from, to, cap, default, cost));
10389
return m;
10490
}
10591

@@ -113,11 +99,9 @@ public int AddEdge(int from, int to, TCap cap, TCost cost)
11399
/// </remarks>
114100
public Edge GetEdge(int i)
115101
{
116-
int m = _pos.Count;
102+
int m = _edges.Count;
117103
DebugUtil.Assert(0 <= i && i < m);
118-
var _e = _g[_pos[i].first][_pos[i].second];
119-
var _re = _g[_e.To][_e.Rev];
120-
return new Edge(_pos[i].first, _e.To, capOp.Add(_e.Cap, _re.Cap), _re.Cap, _e.Cost);
104+
return _edges[i];
121105
}
122106

123107
/// <summary>
@@ -127,16 +111,7 @@ public Edge GetEdge(int i)
127111
/// <para>辺の順番はadd_edgeで追加された順番と同一。</para>
128112
/// <para>計算量: m を追加された辺数として O(m)</para>
129113
/// </remarks>
130-
public List<Edge> Edges()
131-
{
132-
int m = _pos.Count;
133-
var result = new List<Edge>();
134-
for (int i = 0; i < m; i++)
135-
{
136-
result.Add(GetEdge(i));
137-
}
138-
return result;
139-
}
114+
public IReadOnlyList<Edge> Edges() => _edges;
140115

141116
/// <summary>
142117
/// 頂点 <paramref name="s"/> から <paramref name="t"/> へ流せる限り流し、
@@ -271,51 +246,105 @@ public List<Edge> Edges()
271246
DebugUtil.Assert(0 <= s && s < _n);
272247
DebugUtil.Assert(0 <= t && t < _n);
273248
DebugUtil.Assert(s != t);
249+
250+
251+
int m = _edges.Count;
252+
var edgeIdx = new int[m];
253+
254+
CSR<EdgeInternal> g;
255+
{
256+
var degree = new int[_n];
257+
var redgeIdx = new int[m];
258+
var elist = new List<(int from, EdgeInternal edge)>(2 * m);
259+
260+
for (int i = 0; i < m; i++)
261+
{
262+
var e = _edges[i];
263+
edgeIdx[i] = degree[e.From]++;
264+
redgeIdx[i] = degree[e.To]++;
265+
elist.Add((e.From, new EdgeInternal(e.To, -1, capOp.Subtract(e.Cap, e.Flow), e.Cost)));
266+
elist.Add((e.To, new EdgeInternal(e.From, -1, e.Flow, costOp.Minus(e.Cost))));
267+
}
268+
g = new CSR<EdgeInternal>(_n, elist);
269+
for (int i = 0; i < m; i++)
270+
{
271+
var e = _edges[i];
272+
edgeIdx[i] += g.Start[e.From];
273+
redgeIdx[i] += g.Start[e.To];
274+
g.EList[edgeIdx[i]].Rev = redgeIdx[i];
275+
g.EList[redgeIdx[i]].Rev = edgeIdx[i];
276+
}
277+
}
278+
279+
var result = Slope(g, s, t, flowLimit);
280+
281+
for (int i = 0; i < m; i++)
282+
{
283+
var e = g.EList[edgeIdx[i]];
284+
_edges[i].Flow = capOp.Subtract(_edges[i].Cap, e.Cap);
285+
}
286+
287+
return result;
288+
}
289+
290+
291+
private List<(TCap cap, TCost cost)> Slope(CSR<EdgeInternal> g, int s, int t, TCap flowLimit)
292+
{
274293
// variants (C = maxcost):
275294
// -(n-1)C <= dual[s] <= dual[i] <= dual[t] = 0
276-
// reduced cost (= e.cost + dual[e.from] - dual[e.to]) >= 0 for all edge
295+
// reduced cost (= e.cost + dual[e.from] - dual[e.To]) >= 0 for all edge
277296
var dual = new TCost[_n];
278297
var dist = new TCost[_n];
279-
var pv = new int[_n];
280-
var pe = new int[_n];
281-
var vis = new bool[_n];
298+
var prevE = new int[_n];
299+
282300

283301
bool DualRef()
284302
{
285303
dist.AsSpan().Fill(costOp.MaxValue);
286-
pv.AsSpan().Fill(-1);
287-
pe.AsSpan().Fill(-1);
288-
vis.AsSpan().Fill(false);
304+
var vis = new bool[_n];
289305

306+
var queMin = new Stack<int>();
290307
var que = new PriorityQueueForMcf();
308+
291309
dist[s] = default;
292-
que.Enqueue(default, s);
293-
while (que.Count > 0)
310+
queMin.Push(s);
311+
while (queMin.Count > 0 || que.Count > 0)
294312
{
295-
int v = que.Dequeue().to;
313+
int v;
314+
if (queMin.Count > 0)
315+
v = queMin.Pop();
316+
else
317+
v = que.Dequeue().to;
296318
if (vis[v]) continue;
297319
vis[v] = true;
298320
if (v == t) break;
299321
// dist[v] = shortest(s, v) + dual[s] - dual[v]
300322
// dist[v] >= 0 (all reduced cost are positive)
301323
// dist[v] <= (n-1)C
302-
for (int i = 0; i < _g[v].Count; i++)
324+
var dualV = dual[v];
325+
var distV = dist[v];
326+
327+
var gStartCur = g.Start[v];
328+
var gStartNext = g.Start[v + 1];
329+
for (int i = gStartCur; i < gStartNext; i++)
303330
{
304-
var e = _g[v][i];
305-
if (vis[e.To] || capOp.Equals(e.Cap, default)) continue;
306-
// |-dual[e.to] + dual[v]| <= (n-1)C
331+
var e = g.EList[i];
332+
if (EqualityComparer<TCap>.Default.Equals(e.Cap, default)) continue;
333+
// |-dual[e.To] + dual[v]| <= (n-1)C
307334
// cost <= C - -(n-1)C + 0 = nC
308-
TCost cost = costOp.Add(costOp.Subtract(e.Cost, dual[e.To]), dual[v]);
309-
if (costOp.GreaterThan(costOp.Subtract(dist[e.To], dist[v]), cost))
335+
var cost = costOp.Add(costOp.Subtract(e.Cost, dual[e.To]), dualV);
336+
if (costOp.GreaterThan(costOp.Subtract(dist[e.To], distV), cost))
310337
{
311-
dist[e.To] = costOp.Add(dist[v], cost);
312-
pv[e.To] = v;
313-
pe[e.To] = i;
314-
que.Enqueue(dist[e.To], e.To);
338+
var distTo = costOp.Add(distV, cost);
339+
dist[e.To] = distTo;
340+
prevE[e.To] = e.Rev;
341+
if (EqualityComparer<TCost>.Default.Equals(distTo, distV))
342+
queMin.Push(e.To);
343+
else
344+
que.Enqueue(distTo, e.To);
315345
}
316346
}
317347
}
318-
319348
if (!vis[t])
320349
{
321350
return false;
@@ -325,42 +354,42 @@ bool DualRef()
325354
{
326355
if (!vis[v]) continue;
327356
// dual[v] = dual[v] - dist[t] + dist[v]
328-
// = dual[v] - (shortest(s, t) + dual[s] - dual[t]) + (shortest(s, v) + dual[s] - dual[v])
329-
// = - shortest(s, t) + dual[t] + shortest(s, v)
330-
// = shortest(s, v) - shortest(s, t) >= 0 - (n-1)C
357+
// = dual[v] - (shortest(s, t) + dual[s] - dual[t]) +
358+
// (shortest(s, v) + dual[s] - dual[v]) = - shortest(s,
359+
// t) + dual[t] + shortest(s, v) = shortest(s, v) -
360+
// shortest(s, t) >= 0 - (n-1)C
331361
dual[v] = costOp.Subtract(dual[v], costOp.Subtract(dist[t], dist[v]));
332362
}
333-
334363
return true;
335364
}
336365

366+
337367
TCap flow = default;
338368
TCost cost = default;
339-
TCost prevCostPerFlow = costOp.Decrement(default); //-1
340-
var result = new List<(TCap cap, TCost cost)>
341-
{
342-
(flow, cost)
343-
};
369+
TCost prevCostPerFlow = costOp.Decrement(default);
370+
var result = new List<(TCap cap, TCost cost)> { (flow, cost) };
344371
while (capOp.LessThan(flow, flowLimit))
345372
{
346373
if (!DualRef()) break;
347-
TCap c = capOp.Subtract(flowLimit, flow);
348-
for (int v = t; v != s; v = pv[v])
374+
var c = capOp.Subtract(flowLimit, flow);
375+
for (int v = t; v != s; v = g.EList[prevE[v]].To)
349376
{
350-
if (capOp.LessThan(_g[pv[v]][pe[v]].Cap, c))
377+
var c2 = g.EList[g.EList[prevE[v]].Rev].Cap;
378+
if (capOp.LessThan(c2, c))
351379
{
352-
c = _g[pv[v]][pe[v]].Cap;
380+
c = c2;
353381
}
354382
}
355-
for (int v = t; v != s; v = pv[v])
383+
for (int v = t; v != s; v = g.EList[prevE[v]].To)
356384
{
357-
_g[pv[v]][pe[v]].Cap = capOp.Subtract(_g[pv[v]][pe[v]].Cap, c);
358-
_g[v][_g[pv[v]][pe[v]].Rev].Cap = capOp.Add(_g[v][_g[pv[v]][pe[v]].Rev].Cap, c);
385+
var e = g.EList[prevE[v]];
386+
e.Cap = capOp.Add(e.Cap, c);
387+
g.EList[e.Rev].Cap = capOp.Subtract(g.EList[e.Rev].Cap, c);
359388
}
360-
TCost d = costOp.Minus(dual[s]);
389+
var d = costOp.Minus(dual[s]);
361390
flow = capOp.Add(flow, c);
362391
cost = costOp.Add(cost, costOp.Multiply(cast.Cast(c), d));
363-
if (costOp.Equals(prevCostPerFlow, d))
392+
if (EqualityComparer<TCost>.Default.Equals(prevCostPerFlow, d))
364393
{
365394
result.RemoveAt(result.Count - 1);
366395
}
@@ -373,7 +402,7 @@ bool DualRef()
373402
/// <summary>
374403
/// フローを流すグラフの各辺に対応した情報を持ちます。
375404
/// </summary>
376-
public struct Edge : IEquatable<Edge>
405+
public class Edge : IEquatable<Edge>
377406
{
378407
/// <summary>フローが流出する頂点。</summary>
379408
public int From { get; set; }
@@ -422,8 +451,7 @@ public EdgeInternal(int to, int rev, TCap cap, TCost cost)
422451
};
423452

424453
private readonly int _n;
425-
private readonly List<(int first, int second)> _pos;
426-
private readonly List<EdgeInternal>[] _g;
454+
private readonly List<Edge> _edges = new List<Edge>();
427455

428456
private class PriorityQueueForMcf
429457
{
@@ -452,7 +480,7 @@ public void Enqueue(TCost cost, int to)
452480
while (c > 0)
453481
{
454482
int p = (c - 1) >> 1;
455-
if (costOp.Compare(cost, _heap[p].cost) < 0)
483+
if (costOp.LessThan(cost, _heap[p].cost))
456484
{
457485
_heap[c] = _heap[p];
458486
c = p;
@@ -476,12 +504,12 @@ public void Enqueue(TCost cost, int to)
476504
int c = (p << 1) + 1;
477505
while (c < n)
478506
{
479-
if (c != n - 1 && costOp.Compare(_heap[c].cost, _heap[c + 1].cost) > 0)
507+
if (c != n - 1 && costOp.GreaterThan(_heap[c].cost, _heap[c + 1].cost))
480508
{
481509
++c;
482510
}
483511

484-
if (costOp.Compare(_heap[c].cost, item.cost) < 0)
512+
if (costOp.LessThan(_heap[c].cost, item.cost))
485513
{
486514
_heap[p] = _heap[c];
487515
p = c;

0 commit comments

Comments
 (0)