# install.packages(c("dodgr", "igraph"))  # run once if needed
library(dodgr)
library(igraph)

## 1. Define the tree as a directed graph (X upstream of A)
#              X
#              |
#              A
#            /   \
#           B     C
#          / \   / \
#         D   E F   G

graph <- data.frame(
  edge_id  = 1:7,
  from     = c("X", "A","A",  "B","B",  "C","C"),
  to       = c("A", "B","C",  "D","E",  "F","G"),
  weight   = rep(1, 7),   # used by dodgr for routing
  distance = rep(1, 7)    # distance-like metric
)

## 2. Define OD flows at the leaf nodes (demands at D, E, F, G)
from  <- "X"                            # origin (root)
to    <- c("D", "E", "F", "G")          # leaves
flows <- matrix(c(5, 5, 10, 10),        # X->D, X->E, X->F, X->G
                nrow = length(from))

## 3. Use dodgr to aggregate flows on each edge
agg <- dodgr_flows_aggregate(
  graph     = graph,
  from      = from,
  to        = to,
  flows     = flows,
  contract  = FALSE,    # avoid contraction issues on tiny graphs
  norm_sums = FALSE     # keep raw flow sums (not normalised)
)

## 4. Compute pressure drop per edge from the flows
# Simple model: ΔP_e = R_e * Q_e^2
# Here we just set R_e = 1 Pa / (flow_unit^2) for all edges
# (You’d normally use real R values from duct calcs)

agg$R_coeff    <- 1
agg$dp_edge_Pa <- agg$R_coeff * agg$flow^2

# Inspect flows and pressure drops per edge
print(agg[, c("edge_id", "from", "to", "flow", "dp_edge_Pa")])

## 5. Build an igraph object for path analysis
edge_dp <- agg[, c("from", "to", "dp_edge_Pa")]
g <- graph_from_data_frame(edge_dp, directed = TRUE)

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

## 6. Enumerate all root → leaf paths and sum ΔP along each
all_paths <- unlist(lapply(leaves, function(leaf) {
  all_simple_paths(g, from = root, to = leaf)
}), recursive = FALSE)

path_info <- lapply(all_paths, function(p) {
  nodes <- names(p)
  # consecutive node pairs = edges along the path
  path_edges <- cbind(nodes[-length(nodes)], nodes[-1])

  # look up dp_edge_Pa for each edge in this path
  dp_vals <- apply(path_edges, 1, function(e) {
    edge_dp$dp_edge_Pa[edge_dp$from == e[1] & edge_dp$to == e[2]]
  })

  list(
    path       = nodes,
    dp_per_edge = dp_vals,
    total_dp   = sum(dp_vals)
  )
})

## 7. Identify the index run (path with max total ΔP)
total_dps <- sapply(path_info, function(x) x$total_dp)
index_run <- path_info[[which.max(total_dps)]]

cat("\nIndex run (max total ΔP): ",
    paste(index_run$path, collapse = " -> "), "\n", sep = "")
cat("Total ΔP along index run: ",
    index_run$total_dp, " Pa\n", sep = "")