///|
/// Immutable snapshot of an original edge.
pub struct McfEdge[Cap, Cost] {
  from : Int
  to : Int
  cap : Cap
  flow : Cap
  cost : Cost
} derive(Debug, Eq)

///|
priv struct Residual {
  to : Int
  mut cap : Int64
  cost : Int64
} derive(Debug)

///|
/// Minimum-cost flow with nonnegative original costs. Capacity and cost can independently be Int or Int64.
pub struct McfGraph[_, _] {
  priv adjacency : Array[Array[Int]]
  priv residual : Array[Residual]
  priv sources : Array[Int]
  priv mut used : Bool
} derive(Debug)

///|
pub fn[C, W] McfGraph::new(n : Int) -> McfGraph[C, W] {
  @internal.require(0 <= n && n <= 100000000)
  {
    adjacency: Array::makei(n, _ => []),
    residual: [],
    sources: [],
    used: false,
  }
}

///|
pub fn[C : @algebra.FlowInt, W : @algebra.FlowInt] McfGraph::add_edge(
  self : McfGraph[C, W],
  from : Int,
  to : Int,
  cap : C,
  cost : W,
) -> Int {
  let n = self.adjacency.length()
  let cap = @algebra.FlowInt::to_i64(cap)
  let cost = @algebra.FlowInt::to_i64(cost)
  @internal.require(
    0 <= from && from < n && 0 <= to && to < n && cap >= 0L && cost >= 0L,
  )
  let index = self.residual.length()
  self.adjacency[from].push(index)
  self.adjacency[to].push(index + 1)
  self.residual.push({ to, cap, cost, })
  self.residual.push({ to: from, cap: 0L, cost: -cost, })
  self.sources.push(from)
  index / 2
}

///|
pub fn[C : @algebra.FlowInt, W : @algebra.FlowInt] McfGraph::get_edge(
  self : McfGraph[C, W],
  i : Int,
) -> McfEdge[C, W] {
  @internal.require(0 <= i && i < self.sources.length())
  let forward = self.residual[2 * i]
  let reverse = self.residual[2 * i + 1]
  {
    from: self.sources[i],
    to: forward.to,
    cap: @algebra.FlowInt::from_i64(forward.cap + reverse.cap),
    flow: @algebra.FlowInt::from_i64(reverse.cap),
    cost: @algebra.FlowInt::from_i64(forward.cost),
  }
}

///|
pub fn[C : @algebra.FlowInt, W : @algebra.FlowInt] McfGraph::edges(
  self : McfGraph[C, W],
) -> Array[McfEdge[C, W]] {
  Array::makei(self.sources.length(), i => self.get_edge(i))
}

///|
/// Piecewise-linear minimum cost curve, starting at (0,0), with collinear points merged.
/// Call flow/slope at most once per graph. Total flow/cost must fit the chosen types.
pub fn[C : @algebra.FlowInt, W : @algebra.FlowInt] McfGraph::slope(
  self : McfGraph[C, W],
  s : Int,
  t : Int,
  flow_limit? : C,
) -> Array[(C, W)] {
  let n = self.adjacency.length()
  let limit = match flow_limit {
    Some(x) => @algebra.FlowInt::to_i64(x)
    None => C::maximum()
  }
  @internal.require(
    0 <= s && s < n && 0 <= t && t < n && s != t && limit >= 0L && !self.used,
  )
  self.used = true
  let dual = Array::make(n, 0L)
  let previous = Array::make(n, -1)
  let mut flow = 0L
  let mut cost = 0L
  let mut previous_unit_cost = -1L
  let result : Array[(C, W)] = [
    (@algebra.FlowInt::from_i64(0L), @algebra.FlowInt::from_i64(0L)),
  ]
  while flow < limit {
    let distance = Array::make(n, 9223372036854775807L)
    let visited = Array::make(n, false)
    let heap = @internal.MinHeap::new()
    distance[s] = 0L
    heap.push((0L, s))
    while heap.pop() is Some((_, v)) {
      if visited[v] {
        continue
      }
      visited[v] = true
      if v == t {
        break
      }
      for index in self.adjacency[v] {
        let edge = self.residual[index]
        if edge.cap == 0L || visited[edge.to] {
          continue
        }
        let reduced = edge.cost - dual[edge.to] + dual[v]
        // Subtract before comparing so unreachable sentinel + reduced cannot overflow.
        if distance[edge.to] - distance[v] > reduced {
          distance[edge.to] = distance[v] + reduced
          previous[edge.to] = index
          heap.push((distance[edge.to], edge.to))
        }
      }
    }
    if !visited[t] {
      break
    }
    for v = 0; v < n; v = v + 1 {
      if visited[v] {
        dual[v] -= distance[t] - distance[v]
      }
    }
    let mut amount = limit - flow
    let mut v = t
    while v != s {
      let index = previous[v]
      amount = Int64::min(amount, self.residual[index].cap)
      v = self.residual[index ^ 1].to
    }
    v = t
    while v != s {
      let index = previous[v]
      self.residual[index].cap -= amount
      self.residual[index ^ 1].cap += amount
      v = self.residual[index ^ 1].to
    }
    let unit_cost = -dual[s]
    flow += amount
    cost += amount * unit_cost
    if previous_unit_cost == unit_cost {
      ignore(result.pop())
    }
    result.push(
      (@algebra.FlowInt::from_i64(flow), @algebra.FlowInt::from_i64(cost)),
    )
    previous_unit_cost = unit_cost
  }
  result
}

///|
pub fn[C : @algebra.FlowInt, W : @algebra.FlowInt] McfGraph::flow(
  self : McfGraph[C, W],
  s : Int,
  t : Int,
  flow_limit? : C,
) -> (C, W) {
  let points = self.slope(s, t, flow_limit?)
  points[points.length() - 1]
}