# フロー増加法 (flow augumenting method)

# ラベリング法
labeling <- function(R) {
  V <- 1:ncol(R)
# (0)  
  L <- 1 # L の初期値は始点のみ要素
  S <- numeric(0) # Sの初期値は空
  p <- numeric(length(V)) # pの初期値はすべて0
  n <- 0
# (1)    
  while (is.element(V[length(V)], L) == FALSE) {
    cat(n, "L=", L, "S=", S, "p=", p, "\n")
  # t in L なら増加路が見つかったので終了
    if (setequal(L, S) == TRUE) {
      return(numeric(0))
      # L = S なら増加路は存在しないので修了
    }
    else {
      for (i in L) { 
        if (is.element(i, S) == FALSE) {
        # Lの中からSの要素でない点iを一つ選びSに加える
          S <- union(S, i)
          break
        }
      }
# (2)      
      for (j in 1:length(V)) {        
        # 残余ネットワークにおける点iを支援とする枝(i, j)のすべてに対して
        # Lの要素でないjをLに加え、p[j] <- i とする
#        cat("R[", i, ",", j, "]=", R[i, j], "\n")
        if ( is.element(j, L) == FALSE & R[i, j] > 0 ) {
          L <- union(L, j)
          p[j] <- i
        }
      }
    }
    n <- n + 1
  }
  cat(n, "L=", L, "S=", S, "p=", p, "\n")
  return(p)
}


flow_augumenting <- function(G) {
  V <- 1:ncol(G)
  t <- length(V)

# (0)  
  # 残余ネットワーク(residual network)
  R <- G
# 上三角部が残余容量
# 下三角部がフロー
  
  for (i in 1:nrow(R)) {
    cat(R[i,], "\n")
  }
# (1)
  p <- labeling(R)
  while(length(p) > 0) {
    minf <- 1e10 # minf ボトルネック
    d <- t
    j <- t
    while(p[j] > 0) {
      i <- p[j]
      d <- union(i, d)
      if (R[i, j] < minf) {
#        cat("R[", i, ",", j, "]=", R[i, j], "\n")
        minf <- R[i, j]
      }
      j <- i
    }
    cat("p=", p, "d=", d, ", flow=", minf, "\n")
    for (j in 2:length(d)) {
      R[d[j], d[j-1]] <- R[d[j], d[j-1]] + minf # フロー追加
      R[d[j-1], d[j]] <- R[d[j-1], d[j]] - minf # 残余ネットワーク更新
      
    }
    for (i in 1:nrow(R)) {
      cat(R[i,], "\n")
    }
    p <- labeling(R)
  }
  return(R)
}


# ネットワークおよび容量
G <- matrix(c(0, 3, 5, 0, 0, 0, 0,  0, 0, 0, 2, 2, 0, 0,  0, 0, 0, 3, 0, 3, 0,  0, 0, 0, 0, 0, 0, 2,  0, 0, 0, 0, 0, 0, 3,  0, 0, 0, 0, 0, 0, 4,  0, 0, 0, 0, 0, 0, 0), ncol=7, byrow=T)

#
R <- flow_augumenting(G)
cat("\n残余ネットワーク\n")
print(R)

cat("\nフロー\n")
F <- t(R)
for (i in 2:nrow(F)) {
  for (j in 1:(i-1)) {
    F[i, j] <- 0
  }
}
print(F)