再帰、スタックを使わないDFS(深さ優先探索)を書く

スタックを使用したDFSのコードがよくわからないのでスタックを使わないコードをまとめています。

準備

• グラフ G : 隣接リストで表現します。
• 変数 u : 現在の頂点を表現します。
• 配列 active : active[i]が頂点iの探索を終えた辺の数を管理します。
• 配列 parent : parent[i]が頂点iの親を管理します。訪問済みかどうかも管理します。根は自身を親にしてしまいます。

G[i][j]を頂点iでj個の辺の探索を終えたとき、次に探索する辺と読み替えます。

G,u,activeを組み合わせることにより、G[u][active[u]]は現在探索すべき辺として表現されます(以降、探索辺と呼びます)。これとparent,visitedを組み合わせることにより、DFSの主要な動作が表現できます。

1. すべての辺の探索が終わったか(次数と比較) : active[u] >= G[u].size
2. ある頂点が訪問済みか : parent[i]
3. 探索辺の他端へ移動する: u = G[u][active[u]]
4. 探索辺を探索済みとする: active[u] += 1
5. 現在の頂点の親へ移動する: u = parent[u]

注意/補足

辺にコストがある場合など、辺に情報を付与する場合を考えると本来的には辺として理解すべきものだとおもいます。フローのコードなどでは気をつけましょう。

# Edge = Struct.new(:to, :cap)
e = G[u][active[u]]
u = e.to

探索を終えた辺の数を管理するactiveの名前ですが、パスを復元するときには、アクティブな辺として理解するほうが易しいだろうと考えたからです。
始点Sから現在地uまでのパスは以下のようにSからアクティブな辺を辿るイメージで復元できます。

path = []
w = S
until w == u
  path << w
  w = G[w][active[w]]
end

全体のステップ

探索辺があるか?ある場合は探索辺にそって移動するか?
これをループの中心に据えることで、現在地を示すuの場所が移動しDFSを実現します。

1. (初期化処理)現在の頂点uを始点Sに設定し、訪問済とする。
2. 繰り返し(終点Tがある場合、u == T の場合にループを終了)
   1. 現在の頂点の辺がすべて探索済みの場合
      • 現在の頂点が始点Sの場合、すべての探索が終了
      • 現在の移動元(親)に移動し、探索辺を探索済とする。
   2. 探索辺がある場合、隣接する頂点が存在する。
      • 未訪問の場合、現在の頂点を移動先の頂点の親とし、訪問済にして移動する。
      • 訪問済の場合、探索辺を探索済とする。

このステップをrubyのコードにしたものが以下になります。
スタックをつかったわかりにくさがなくなり、DFSの探索をそのまま表現できていると思っています。

# Gは隣接リストとして表現されたグラフ
G = Array.new(4){ [] }
G[0] << 1 << 2; G[1] << 3; G[2] << 3

# Sは始点
S = 0

# DFSで使用するparent, active, u
active = Array.new(G.size, 0)
parent = Array.new(G.size)

u = S
parent[S] = S # 根の親の初期化は忘れない

# 終点Tがある場合はループの条件を設定する
# until u == T
# 注) loopで書いているが、while trueの方が速い。
loop do 
  # 探索済みの辺の数を確認する
  if (j = active[u]) >= G[u].size
    # 次数以上なら、現在の頂点の探索は終了
    break if u == S # 始点の場合は探索が終了
    # 親に戻り、辺の探索を終了する
    u = parent[u] 
    active[u] += 1
  else
    v = G[u][j]
    if !parent[v]
      # 探索辺があり、隣接する頂点が未訪問なので移動する
      # 移動の際、親を設定し、訪問済にする
      parent[v] = u
      u = v
    else
      # 探索辺はあるが、隣接する頂点は既に訪問済
      # この辺の探索を終了する
      active[u] += 1
    end
  end
end

例.1 オイラーツアー、訪問順、帰りがけ順など各種順序の構築

次のコードではオイラーツアー、訪問順、帰りがけ順、DFS木での深さを求めます。
訪問順は初めて訪問するとき、帰りがけ順は親に戻るとき、オイラーツアーは移動のたびにそれぞれ頂点をリストに追加することで実現しています。

LCA(最深共通先祖)の下準備などで使用します。
スタックを使った実装だとスタック上に頂点の重複が発生することもあり、帰りがけ順やオイラーパスをつくるときに不安な部分があったのですが、その点が払拭されました。

G = Array.new(4){ [] }
G[0] << 1 << 2
G[1] << 3
G[2] << 3

parent = Array.new(G.size)
active = Array.new(G.size, 0)

# オプション
# depth: 各頂点の深さ
# euler_tour, pre_order, post_order: 各種順序
depth = Array.new(G.size) 
euler_tour = [] 
pre_order = []
post_order = []

s = G.size

# 今回はGは連結とは限らない前提とし、各頂点からDFSを行う
G.each_index do |s|
  next if parent[s] # 既に訪問済みなら探索は始めない

  parent[s] = s 
  depth[s] = 0 # ルートの深さは 0

  u = s # 現在地を s にする

  # 訪問順、オイラーツアーの最初にsを追加
  pre_order << u 
  euler_tour << u

  loop do
    if (j = active[u]) >= G[u].size
      # 帰りがけ順に現在の頂点を追加
      post_order << u

      break if u == s
      u = parent[u]
      active[u] += 1

      # 移動先の頂点をオイラーツアーに追加
      euler_tour << u
    else
      v = G[u][j]
      if !parent[v]
        # 深さを設定
        depth[v] = depth[u] + 1
        parent[v] = u
        u = v
  
        # 移動先の頂点を訪問順とオイラーツアーに追加
        pre_order << u
        euler_tour << u
      else
        active[u] += 1
      end
    end
  end
end

puts "pre_order : #{pre_order}" 
puts "post_order : #{post_order}" 
puts "euler_tour : #{euler_tour}"
puts "parent: #{parent}"
puts "depth: #{depth}"

例2. 最大流問題

Dinic法のコードです。

何度も増加道を探す都合上、再帰で書いたDFSでも継続して探索できるような工夫が必要ですが、このDFSだとactiveをそのまま利用すればよいだけです。

• 最大流問題について

蟻本(p. 194) に書いてある
無駄な辺を何度も調べない
というのは, dfs内のループが iter[v] をインデックスとして回っている事に対応しています.

class MaxFlow
  Edge = Struct.new(:to, :cap, :rev)
  def max_flow(s,t)
    flow = 0
    active = Array.new(@g.size)
    level = Array.new(@g.size)
    parent = Array.new(@g.size)
    while levelize(s,t,level)
      active.fill(0)
      while f = find_flow(s, t, level, active, parent)
        flow += f
      end
    end
    flow
  end
  
  def find_flow(s, t, level, active, parent)
    parent.fill(nil)
    u = parent[s] = s
    until u == t
      if (j = active[u]) >= @g[u].size
        return false if u == s
        u = parent[u]
        active[u] += 1
      else
        e = @g[u][j]
        if e.cap > 0 && level[e.to] < level[u]
          parent[e.to] = u
          u = e.to
        else
          active[u] += 1
        end
      end
    end
    u, f = s, Float::INFINITY
    until u == t
      fwd = @g[u][active[u]]
      f = _min(fwd.cap, f)
      u = fwd.to
    end
    u = s
    until u == t
      fwd = @g[u][active[u]]
      rev = @g[fwd.to][fwd.rev]
      fwd.cap -= f
      rev.cap += f
      u = fwd.to
    end
    return f
  end

  def levelize(s, t, level)
    level.fill(@g.size * 2)
    level[t] = 0
    q = [t]
    until q.empty?
      u = q.shift
      d = level[u]
      next if u == s
      @g[u].each do |e|
        v = e.to
        rev = @g[e.to][e.rev] 
        # tからbfsしているので逆辺の容量を確認
        if rev.cap > 0 && level[v] > d + 1
          level[v] = d + 1
          q << v
        end
      end
    end
    return false if level[s] > @g.size
    level[s] = @g.size
  end
  # 残余ネットワークで到達可能かどうかを返す
  # G.max_flow(s,t)
  # reachables = G.min_cust(s)
  # blocked = E.select{|from,to| reachables[from] & !reacheables[to] }
  # カットの復元につかう(see. ABC239 G)
  def min_cut(s)
    used = Array.new(@g.size, false)
    q = [s]
    used[s] = true
    until q.empty?
      u = q.shift
      @g[u].each do |e|
        next if e.cap.zero? || used[e.to]
        used[e.to] = true
        q << e.to
      end
    end
    used
  end

  def initialize(n); @g = Array.new(n){ [] }; end
  def size; @g.size; end
    
  def add_edge(from, to, cap)
    fwd, rev = @g[from].size, @g[to].size
    @g[from] << Edge.new(to, cap, rev)
    @g[to] << Edge.new(from, 0, fwd)
    self
  end
  
  def link(u, v, cap)
    fwd, rev = @g[u].size, @g[v].size
    @g[u] << Edge.new(v, cap, rev)
    @g[v] << Edge.new(u, cap, fwd)
    self
  end
  def _min(a,b); a < b ? a : b; end
end

例3. 2部グラフマッチング

2部グラフマッチングは交差路のため、頂点２つ分戻る点が異なります
探索部分は定義をそのままコードに落とした感じです。

マッチング M について以下を満たすとき道 P は増加道であるという。

P の出発点と到着点がマッチングの端点でない
交互路である（M に含まれていない枝と含まれている枝を交互に使う）
参考 : 二部グラフの最大マッチングと増加道

class BiparateMatching
  def initialize(n)
    @g = Array.new(n){[]}
  end
  def add_edge(u,v)
    raise if u == v
    u = (u.is_a?(Enumerable) ? u.to_a : [u])
    v = (v.is_a?(Enumerable) ? v.to_a : [v])
    u.product(v){|a,b| @g[a] << b; @g[b] << a }
    self
  end
  
  def biparate_matching
    matching = Array.new(@g.size) # マッチングを記録
    active = Array.new(@g.size, 0)
    parent = Array.new(@g.size)

    @g.each_index do |s|
      next if matching[s] # 頂点sにペアがいれば探さない 
      find_matching(s, matching, active.fill(0), parent.fill(nil))
    end
    matching
  end
  alias_method :bi_matching, :biparate_matching
  
  def find_matching(s, matching, active, parent)
    u = parent[s] =  s
  
    # 現在の頂点uがs以外のマッチングのない端点ならば、増加道が見つかった。
    until u != s && !matching[u]
      if (j = active[u]) >= @g[u].size
        return false if u == s
        # 増加道の探索中は必ずmatchingがあるため一緒にもどる。
        u = parent[matching[u]]
        active[u] += 1
      else
        v = @g[u][j]
        if !parent[v]
          # 移動先にマッチングがある場合は、
          # 交差路となるように移動する必要があるため、マッチング辺もあわせて移動する。
          if w = matching[u]
            parent[v] = u
            parent[w] = v
            u = w
          else
            parent[v] = u
            u = v
          end
        else
          active[u] += 1
        end
      end
    end
    # 増加道の終端をtとし、s - tパスに沿ってマッチングを更新
    # 交差路なので、以前のmatchingと交互に移動することに注意
    t,u,v = u,s,nil
    until v == t
      v = @g[u][active[u]]
      w = matching[v]
      matching[u] = v
      matching[v] = u
      u = w
    end
    return true
  end
end

G = BiparateMatching.new(9)
G.add_edge(1,[5,6,7,8])
   .add_edge([0,2,3],6)
   .add_edge(4,5)
   
p G.biparate_matching

例4. 強連結成分分解

元のグラフと逆辺のグラフで２回のDFSを行います。
始点毎に帰りがけ順の配列で返すことによってスッキリかけます。
帰りがけ順の作り方にも不安はありません。

class SCCGraph
  def initialize(n)
    @fwd = Array.new(n){ [] }
    @rev = Array.new(n){ [] }
  end

  def add_edge(from, to)
    @fwd[from] << to
    @rev[to] << from
    self
  end

  # 強連結成分分解を行います
  # 実行後、components, idsから結果を取得できる
  def scc!
    # fwd_components.flatten!.reverse!が
    # 元のグラフ(@fwd)のreverse_post_order
    fwd_components = decomposite(@fwd, @fwd.each_index)

    # reverse_post_orderにしたがって逆辺のグラフ(@rev)を分解
    @components = decomposite(@rev, fwd_components.flatten!.reverse!)
    
    # 頂点ごとの強連結成分のインデクスをまとめる
    @ids = Array.new(@rev.size)
    @components.each_with_index do |connected, i|
      connected.each do |u|
        @ids[u] = i
      end
    end
    self
  end
  # 強連結成分毎にまとまった頂点の配列 : [[1,2,0],[3],[4]]
  attr_reader :components

  # 頂点を含む強連結成分のインデクス : [0,0,0,1,2]
  attr_reader :ids

  # 与えられたグラフ(g)をorderに従って、DFSした結果を返します。
  def decomposite(g, order)
    active = Array.new(g.size, 0)
    parent = Array.new(g.size)
    components = []
    order.each do |s| 
      next if parent[s]
      components << post_ordered_component(g, s, active, parent)
    end
    components
  end

  # s から到達できる頂点を帰りがけ順の配列にして返します
  def post_ordered_component(g, s, active, parent)
    connected = []
    u = parent[s] = s
    while true
      if (j = active[u]) >= g[u].size
        connected << u
        break if u == s
        u = parent[u]
        active[u] += 1
      else
        v = g[u][j]
        if !parent[v]
          parent[v] = u
          u = v
        else
          active[u] += 1
        end
      end
    end
    return connected
  end

end

N, M = gets.split.map(&:to_i)
G = SCCGraph.new(N)
M.times do
  a,b = gets.split.map(&:to_i)
  G.add_edge(a, b)
end

scc = G.scc!.components
puts scc.size
puts scc.map{|connected| "#{connected.size} #{connected.join(' ')}" }

例5. 2-SAT

SCCができると2-SATも解けます。

class TwoSAT
  def initialize(n)
    @size = n
    @fwd = Array.new(2 * n){ [] }
    @rev = Array.new(2 * n){ [] }
  end

  def add_clause(i, f, j, g)
    add_edge(2 * i + (f ? 0 : 1), 2 * j + (g ? 1 : 0))
    add_edge(2 * j + (g ? 0 : 1), 2 * i + (f ? 1 : 0))
  end
  
  def satisfiable
    scc!
    # 変数の真偽は連続した2頂点となっているので、t,fは真,偽それぞれが含まれる強連結成分
    ids.each_slice(2).map do |t, f|
      return false if t == f # ある変数の真偽が同じ強連結成分に含まれている場合は矛盾
      t < f
    end
  end

  private

  def add_edge(from, to)
    @fwd[from] << to
    @rev[to] << from
    self
  end

  def scc!
    fwd_components = decomposite(@fwd, @fwd.each_index)
    @components = decomposite(@rev, fwd_components.flatten!.reverse!)
    
    # 各頂点がどの強連結成分に含まれているかを示す配列
    @ids = Array.new(@rev.size)
    @components.each_with_index do |connected, i|
      connected.each do |u|
        @ids[u] = i
      end
    end
    self
  end

  attr_reader :components
  attr_reader :ids

  def decomposite(g, order)
    active = Array.new(g.size, 0)
    parent = Array.new(g.size)
    components = []
    for s in order 
      next if parent[s]
      components << post_ordered_component(g, s, active, parent)
    end
    components
  end

  def post_ordered_component(g, s, active, parent)
    connected = []
    u = parent[s] = s
    while true
      if (j = active[u]) >= g[u].size
        connected << u
        break if u == s
        u = parent[u]
        active[u] += 1
      else
        v = g[u][j]
        if !parent[v]
          parent[v] = u
          u = v
        else
          active[u] += 1
        end
      end
    end
    return connected
  end

end

n, d = gets.split.map(&:to_i)
x, y = n.times.map { gets.split.map(&:to_i) }.transpose
ts = TwoSAT.new(n)

n.times do |i|
  (i + 1...n).each do |j|
    ts.add_clause(i, false, j, false) if (x[i] - x[j]).abs < d
    ts.add_clause(i, false, j, true ) if (x[i] - y[j]).abs < d
    ts.add_clause(i, true,  j, false) if (y[i] - x[j]).abs < d
    ts.add_clause(i, true,  j, true ) if (y[i] - y[j]).abs < d
  end
end

if answer = ts.satisfiable
  puts 'Yes'
  puts answer.each_with_index.map{|ans, i| ans ? x[i] : y[i] }
else
  puts 'No'
end

例6. 2辺連結成分分解

2辺連結成分分解も強連結成分分解と同じように書いてみた。
一言でいうと、無向グラフを(必ず逆辺とセットの)有向グラフに変換して、DFS木を取り除いたグラフで強連結成分分解すると2辺連結成分分解

通常2辺連結成分分解は low-link法で実装するとおもいますし、そちらの方が実装も簡単です。

class TeCCGraph
  
  # 辺を使用して移動した場合 usedをtrueにする
  Edge = Struct.new(:to, :used)
  
  def initialize(n); @g = Array.new(n){ [] }; end
  def add_edge(u, v)
    @g[u] << Edge.new(v, false)
    @g[v] << Edge.new(u, false)
    self
  end

  
  # ２辺連結成分分解
  # 実行後、components, idsから結果を取得できる
  def tecc!
    fwd_components = decomposite(@g, @g.each_index)
    @components = decomposite(@g, fwd_components.flatten!.reverse!)
    @ids = Array.new(@g.size)
    @components.each_with_index do |connected, i|
      connected.each do |u|
        @ids[u] = i
      end
    end
    self
  end
  # 2辺連結成分毎にまとまった頂点の配列 : [[1,2,0],[3],[4]]
  attr_reader :components

  # 頂点を含む2辺連結成分のインデクス : [0,0,0,1,2]
  attr_reader :ids

  def decomposite(g, order)
    parent = Array.new(g.size)
    active = Array.new(g.size, 0)
    components = []
    for s in order
      next if parent[s]
      components << post_ordered_component(g, s, parent, active)
    end
    components
  end

  def post_ordered_component(g, s, parent, active)
    u = parent[s] = s
    connected = []
    while true
      if (j = active[u]) >= g[u].size
        connected << u
        break if u == s
        u = parent[u]
        active[u] += 1
      else
        e = g[u][j]
        v = e.to
        # 未訪問かつ、この辺が未使用の場合に移動し、usedをtrueに設定する。
        # このため２回目のDFSでは親か1回目のDFSでの後退辺でしか移動できない。
        # ２回目のDFSでは後退辺を辿ることで２辺連結成分となる
        if !parent[v] && !e.used
          parent[v] = e.used = u
          u = v
        else
          active[u] += 1
        end
      end
    end
    connected
  end
end

G = TeCCGraph.new(5)
G.add_edge(0,1)
  .add_edge(1,2)
  .add_edge(2,0)
  .add_edge(2,3)
  .add_edge(3,4)

G.tecc!
p G.components
p G.ids

例7. オイラーツアーとセグメント木のLCA(最深共通祖先)

セグメント木も非再帰で書いたら、2べきにあわせなくてよいし、とても楽ちん

class EulerTourLCA
  class << self
    # オイラーツアー(tour)と各頂点が初めて現れたときのインデクス(ord)
    # 深さ(depth)を求め、セグメント木を作る
    def create(g, s = 0)
      ord = Array.new(g.size)
      parent = Array.new(g.size)
      active = Array.new(g.size, 0)
      depth = Array.new(g.size)
      tour = []
      
      u = parent[s] = s
      depth[s] = 0
      ord[s] = tour.size
      tour << s
      
      while true
        if (j = active[u]) >= g[u].size
          break if u == s
          u = parent[u]
          tour << u
          active[u] += 1
        else
          v = g[u][j]
          if !parent[v]
            parent[v] = u
            depth[v] = depth[u] + 1
            u = v
            ord[u] = tour.size
            tour << u
          else
            active[u] += 1
          end
        end
      end
      new(tour, depth, ord)
    end
  end

  def initialize(tour, depth, ord)
    @depth, @ord, @offset, k = depth, ord, tour.size, tour.size * 2
    # オイラーツアーを元にセグメント木のデータを作る2冪でなくてよい
    @data = (tour * 2).fill(-1, 0, @offset)
    @data[k >> 1] = min_depth(@data[k >> 1], @data[k]) while (k -= 1) > 0
  end
  # a,b間の距離 
  def dist(a,b); @depth[a] + @depth[b] - 2 * @depth[lca(a,b)]; end

  # オイラーツアー上の深さ最小の頂点をセグメント木で求める
  def lca(a,b)
    a,b = b,a if @ord[a] > @ord[b] 
    find_min(@ord[a], @ord[b] + 1)
  end

  private
  # a,bのうち浅い方を返す
  def min_depth(a,b)
    raise if a >= @depth.size || b >= @depth.size
    return b if a < 0
    return a if b < 0
    @depth[a] < @depth[b] ? a : b
  end
 
  # 非再帰のセグメント木探索
  def find_min(a, b)
    l, r, x = @offset + a, @offset + b, -1
    while l < r
      if l.odd?
        x = min_depth(x, @data[l]) 
        l += 1
      end
      if r.odd?
        r -= 1
        x = min_depth(x, @data[r])
      end
      l >>= 1
      r >>= 1
    end
    return x
  end

end

N = gets.to_i
G = Array.new(N + 1){ [] }
(N - 1).times do
  u, v, d = gets.split.map(&:to_i)
  G[u] << [v, d]
  G[v] << [u, d]
end

tour = EulerTourLCA.create(G, 1)

Q = gets.to_i
puts Q.times.map{ tour.dist(*gets.split.map(&:to_i)) }

感想など

rubyで競プロやると再帰で結構こけることが多いのでなるべくループで書けるように練習しています。
木を探索しようとして入力が1本道になっているような入力で再帰のDFSは簡単に落とされていました。
(昔はスタックサイズの指定をしていなかったのかな? いまはたぶん落ちないです。)

DFSの非再帰版はスタック使うと書いてあるのですが、スタックを使った場合は訪問順は分かるのですが、(頭悪いだけですが)どのタイミングで親に戻っているのか、よくわかっていませんでした。

以前は以下のようなコードで理解しようとしていたんだとおもいます。BFSと似たような書き方を意識しすぎていたのが今となっては分かります。
(再帰とは逆ですが)訪問順で処理できるのですが、親や深さなど全然わからない状態でした。

  def dfs(s, &block)
    stack = [s]
    used = Array.new(@g.size)
    until u = stack.pop
      next if used[u]
      used[u] = true
      // ここにそのノードの処理
      // block.call(u)
      @g[u].each do |e|
        next if !used[e.to]
        stack << e.to 
      end
    end
  end

いまWikipediaの説明をみるとスタックに積むのはイテレータであってノードではないというのが重要だとわかります。
疑似コードのほうはそれでよいのか?という気もしますが、Pythonのコードは理解できる気がします。

最近やっとスタックを使った書き方がしっくりしてきました。
現在は以下のような形の理解で落ち着いています。

class Graph
  Edge = Struct.new(:to)
  def initialize(n)
    @g = Array.new(n){ [] }
  end
  def add_edge(a,b)
    @g[a] << Edge.new(b)
    @g[b] << Edge.new(a)
    self
  end
  def dfs(s)
    stack = [s]
    active = Array.new(@g.size, -1)
    while u = stack[-1]
      active[u] += 1
      e = @g[u][active[u]]
      if e.nil?
        stack.pop
      elsif active[e.to] < 0
        stack << e.to
      end
    end
  end
end