library(igraph)

# ------------------------------------------------------------
# 1. Define your tree edges with pressure drops already calculated
# ------------------------------------------------------------

edges <- data.frame(
  from = c("X","A","A","B","B","C","C"),
  to   = c("A","B","C","D","E","F","G"),
  dp   = c(2, 1, 3, 0.5, 0.5, 1.2, 1.2)  # example ΔP values per edge
)

# Build directed graph
g <- graph_from_data_frame(edges, directed = TRUE)

# ------------------------------------------------------------
# 2. Identify root and leaves
# ------------------------------------------------------------

root <- "X"
# leaves = nodes with no outgoing edges
leaves <- V(g)[degree(g, mode = "out") == 0]$name

# ------------------------------------------------------------
# 3. Extract all root → leaf paths
# ------------------------------------------------------------

all_paths <- lapply(leaves, function(leaf) {
  # list of igraph vertex sequences (usually just 1 path in a tree)
  paths <- all_simple_paths(g, from = root, to = leaf)
  return(paths)
})

# Flatten list structure
all_paths <- unlist(all_paths, recursive = FALSE)

# ------------------------------------------------------------
# 4. Convert each path to edge list and sum pressure drops
# ------------------------------------------------------------

path_info <- lapply(all_paths, function(p) {
  # convert vertex sequence to character vector
  nodes <- names(p)
  
  # get matrix of consecutive edges along the path
  path_edges <- cbind(nodes[-length(nodes)], nodes[-1])
  
  # find pressure drop for each path edge
  dp_vals <- apply(path_edges, 1, function(e) {
    edges$dp[edges$from == e[1] & edges$to == e[2]]
  })
  
  # total pressure drop for this path
  total_dp <- sum(dp_vals)
  
  list(
    path = nodes,
    dp_per_edge = dp_vals,
    total_dp = total_dp
  )
})

# ------------------------------------------------------------
# 5. Identify **index run** = path with maximum ∑ΔP
# ------------------------------------------------------------

total_dps <- sapply(path_info, function(x) x$total_dp)
index_run <- path_info[[ which.max(total_dps) ]]

# ------------------------------------------------------------
# 6. Print index run
# ------------------------------------------------------------

cat("Index run: ", paste(index_run$path, collapse=" → "), "\n")
cat("Total ΔP: ", index_run$total_dp, "Pa\n")