tests/testthat/test-operations.R

In caugi: Causal Graph Interface

# ──────────────────────────────────────────────────────────────────────────────
# ───────────────────────────────── Skeleton ───────────────────────────────────
# ──────────────────────────────────────────────────────────────────────────────

test_that("DAG --> UG with skeleton keeps the skeleton", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    C %-->% D,
    D %-->% E,
    F,
    G,
    class = "DAG"
  )

  skel_cg <- skeleton(cg)

  expect_equal(edge_types(skel_cg), "---")
  expect_equal(length(skel_cg), 7)
  expect_equal(nodes(cg), nodes(skel_cg))
  expect_equal(edges(cg)$from, edges(skel_cg)$from)
  expect_equal(edges(cg)$to, edges(skel_cg)$to)
})

test_that("PDAG --> UG with skeleton keeps the skeleton", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    C %---% D,
    D %-->% E,
    F,
    G,
    class = "PDAG"
  )

  skel_cg <- skeleton(cg)

  expect_equal(edge_types(skel_cg), "---")
  expect_equal(length(skel_cg), 7)
  expect_equal(nodes(cg), nodes(skel_cg))
  expect_equal(edges(cg)$from, edges(skel_cg)$from)
  expect_equal(edges(cg)$to, edges(skel_cg)$to)
})

# ──────────────────────────────────────────────────────────────────────────────
# ──────────────────────────────── Moralizing ──────────────────────────────────
# ──────────────────────────────────────────────────────────────────────────────

test_that("moralize works on DAGs", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    D %-->% C,
    E %-->% C,
    D %-->% B,
    F,
    G,
    nodes = c("A", "B", "C", "D", "E", "F", "G"),
    class = "DAG"
  )

  moralized_cg <- moralize(cg)

  expect_equal(edge_types(moralized_cg), "---")
  expect_equal(length(moralized_cg), 7)
  expect_equal(nodes(cg), nodes(moralized_cg))
  expect_equal(moralized_cg@graph_class, "UG")

  expected_moral_cg <- caugi(
    A %---% B,
    A %---% D,
    B %---% C,
    B %---% D,
    B %---% E,
    C %---% D,
    C %---% E,
    D %---% E,
    F,
    G,
    nodes = c("A", "B", "C", "D", "E", "F", "G"),
    class = "UG"
  )
  expect_equal(edges(moralized_cg), edges(expected_moral_cg))
})

test_that("moralize fails on non-DAGs", {
  cg_pdag <- caugi(
    A %-->% B,
    B %-->% C,
    C %---% D,
    class = "PDAG"
  )

  expect_error(
    moralize(cg_pdag),
    "moralize\\(\\) can only be applied to DAGs\\."
  )

  cg_ug <- caugi(
    A %---% B,
    B %---% C,
    class = "UG"
  )

  expect_error(
    moralize(cg_ug),
    "moralize\\(\\) can only be applied to DAGs\\."
  )
})

############# NetworkX tests for moralization #############
# https://github.com/networkx/networkx/blob/main/networkx/algorithms/tests/test_moral.py

test_that("NetworkX moralize test 1", {
  cg <- caugi(
    A %-->% B,
    C %-->% B,
    D %-->% A,
    D %-->% E,
    F %-->% E,
    G %-->% E,
    class = "DAG"
  )
  moral_cg <- moralize(cg)
  expect_equal(moral_cg@graph_class, "UG")

  has_edge <- function(edges, u, v) {
    any(
      (edges$from == u & edges$to == v) |
        (edges$from == v & edges$to == u)
    )
  }

  expect_true(has_edge(moral_cg@edges, "A", "C"))
  expect_true(has_edge(moral_cg@edges, "D", "F"))
  expect_true(has_edge(moral_cg@edges, "F", "G"))
  expect_true(has_edge(moral_cg@edges, "D", "G"))
  expect_false(has_edge(moral_cg@edges, "A", "E"))
})

# ──────────────────────────────────────────────────────────────────────────────
# ───────────────────────────────── Mutation ───────────────────────────────────
# ──────────────────────────────────────────────────────────────────────────────

test_that("mutate_caugi throws error, when mutation isn't possible", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    C %---% D,
    D %-->% E,
    F,
    G,
    class = "PDAG"
  )

  is_it <- is_dag(cg)
  expect_false(is_it)
  expect_error(mutate_caugi(cg, class = "DAG"), "Cannot convert caugi of class")
})

test_that("mutate_caugi keeps simple = FALSE where the target class allows it", {
  # A bow is a valid, non-simple ADMG.
  cg <- caugi(A %-->% B, A %<->% B, class = "UNKNOWN", simple = FALSE)
  cg_admg <- mutate_caugi(cg, class = "ADMG")
  expect_false(cg_admg@simple)
  expect_silent(build(cg_admg))
  expect_equal(children(cg_admg, "A"), "B")
  expect_equal(spouses(cg_admg, "A"), "B")

  cg_back <- mutate_caugi(cg_admg, class = "UNKNOWN")
  expect_false(cg_back@simple)
  expect_silent(build(cg_back))

  # Simple-only targets become simple.
  cg_dag_like <- caugi(A %-->% B, B %-->% C, class = "UNKNOWN", simple = FALSE)
  cg_dag <- mutate_caugi(cg_dag_like, class = "DAG")
  expect_true(cg_dag@simple)
  expect_silent(build(cg_dag))
  expect_equal(children(cg_dag, "A"), "B")

  # Simple sources stay simple.
  expect_true(mutate_caugi(caugi(A %-->% B), class = "ADMG")@simple)
})

test_that("mutate_caugi works for all classes to UNKNOWN", {
  cg_pdag <- caugi(
    A %-->% B,
    B %-->% C,
    C %---% D,
    D %-->% E,
    F,
    G,
    class = "PDAG"
  )
  cg_dag <- caugi(
    A %-->% B,
    B %-->% C,
    D %-->% E,
    F,
    G,
    class = "DAG"
  )
  cg_ug <- caugi(
    A %---% B,
    B %---% C,
    C %---% D,
    D %---% E,
    F,
    G,
    class = "UG"
  )
  cg_unknown_pdag <- mutate_caugi(cg_pdag, class = "UNKNOWN")
  cg_unknown_dag <- mutate_caugi(cg_dag, class = "UNKNOWN")
  cg_unknown_ug <- mutate_caugi(cg_ug, class = "UNKNOWN")
  expect_equal(cg_unknown_pdag@graph_class, "UNKNOWN")
  expect_equal(cg_unknown_dag@graph_class, "UNKNOWN")
  expect_equal(cg_unknown_ug@graph_class, "UNKNOWN")

  expect_equal(edges(cg_unknown_pdag), edges(cg_pdag))
  expect_equal(edges(cg_unknown_dag), edges(cg_dag))
  expect_equal(edges(cg_unknown_ug), edges(cg_ug))

  expect_equal(nodes(cg_unknown_pdag), nodes(cg_pdag))
  expect_equal(nodes(cg_unknown_dag), nodes(cg_dag))
  expect_equal(nodes(cg_unknown_ug), nodes(cg_ug))
})

test_that("mutate_caugi works from UNKNOWN to DAG", {
  cg_unknown <- caugi(
    A %-->% B,
    B %-->% C,
    D %-->% E,
    F,
    G,
    class = "UNKNOWN"
  )
  is_it <- is_dag(cg_unknown)
  expect_true(is_it)

  cg_dag <- mutate_caugi(cg_unknown, class = "DAG")
  expect_equal(cg_dag@graph_class, "DAG")
  expect_equal(edges(cg_dag), edges(cg_unknown))
  expect_equal(nodes(cg_dag), nodes(cg_unknown))
})

test_that("mutate_caugi works from DAG to PDAG", {
  cg_dag <- caugi(
    A %-->% B,
    B %-->% C,
    D %-->% E,
    F,
    G,
    class = "DAG"
  )

  cg_pdag <- mutate_caugi(cg_dag, class = "PDAG")
  expect_equal(cg_pdag@graph_class, "PDAG")
  expect_equal(edges(cg_pdag), edges(cg_dag))
  expect_equal(nodes(cg_pdag), nodes(cg_dag))
})

test_that("mutate_caugi supports MPDAG target when valid", {
  cg_pdag <- caugi(
    A %---% B,
    A %-->% C,
    B %-->% C,
    class = "PDAG"
  )
  expect_true(is_mpdag(cg_pdag))
  cg_mpdag <- mutate_caugi(cg_pdag, class = "MPDAG")
  expect_equal(cg_mpdag@graph_class, "MPDAG")
  expect_equal(edges(cg_mpdag), edges(cg_pdag))

  cg_not_mpdag <- caugi(
    A %-->% B,
    B %---% C,
    class = "PDAG"
  )
  expect_false(is_mpdag(cg_not_mpdag))
  expect_error(
    mutate_caugi(cg_not_mpdag, class = "MPDAG"),
    "Cannot convert caugi"
  )
})

test_that("mutate_caugi works from PDAG to DAG if PDAG is a DAG", {
  cg_pdag <- caugi(
    A %-->% B,
    B %-->% C,
    D %-->% E,
    F,
    G,
    class = "PDAG"
  )
  is_it <- is_dag(cg_pdag)
  expect_true(is_it)

  cg_dag <- mutate_caugi(cg_pdag, class = "DAG")
  expect_equal(cg_dag@graph_class, "DAG")
  expect_equal(edges(cg_dag), edges(cg_pdag))
  expect_equal(nodes(cg_dag), nodes(cg_pdag))
})

test_that("mutate_caugi works on empty caugi", {
  cg_empty <- caugi(class = "DAG")
  cg_empty_ug <- mutate_caugi(cg_empty, class = "UG")
  expect_equal(length(cg_empty_ug), 0)
  expect_equal(cg_empty_ug@graph_class, "UG")

  cg_empty_mpdag <- mutate_caugi(cg_empty, class = "MPDAG")
  expect_equal(length(cg_empty_mpdag), 0)
  expect_equal(cg_empty_mpdag@graph_class, "MPDAG")
})

test_that("mutate_caugi doesn't change class if old class is equal to new class", {
  cg_dag <- caugi(
    A %-->% B,
    B %-->% C,
    D %-->% E,
    F,
    G,
    class = "DAG"
  )

  cg_dag_mutated <- mutate_caugi(cg_dag, class = "DAG")
  expect_equal(cg_dag_mutated@graph_class, "DAG")
  expect_equal(edges(cg_dag_mutated), edges(cg_dag))
  expect_equal(nodes(cg_dag_mutated), nodes(cg_dag))
})

# ──────────────────────────────────────────────────────────────────────────────
# ──────────────────────────── Latent Projection ───────────────────────────────
# ──────────────────────────────────────────────────────────────────────────────

test_that("latent_project basic confounding: U -> X, U -> Y, X -> Y", {
  # DAG: U -> X, U -> Y, X -> Y
  # Project out U:
  # - U has Ch(U) = {X, Y}
  # - Step 3: Add X <-> Y (children of U pair up)
  # Result: X -> Y AND X <-> Y
  dag <- caugi(
    U %-->% X,
    U %-->% Y,
    X %-->% Y,
    class = "DAG"
  )

  admg <- latent_project(dag, latents = "U")

  expect_equal(admg@graph_class, "ADMG")
  expect_equal(length(admg), 2)
  expect_equal(sort(nodes(admg)$name), c("X", "Y"))

  # Check directed edge X -> Y preserved
  expect_equal(parents(admg, "Y"), "X")
  expect_equal(children(admg, "X"), "Y")

  # X <-> Y added
  expect_equal(spouses(admg, "X"), "Y")
  expect_equal(spouses(admg, "Y"), "X")
})

test_that("latent_project with no latents returns the same graph as ADMG", {
  dag <- caugi(
    X %-->% Y,
    Y %-->% Z,
    class = "DAG"
  )

  admg <- latent_project(dag, latents = character(0))

  expect_equal(admg@graph_class, "ADMG")
  expect_equal(length(admg), 3)
  expect_equal(nodes(admg)$name, c("X", "Y", "Z"))

  # Same directed edges
  expect_equal(parents(admg, "Y"), "X")
  expect_equal(parents(admg, "Z"), "Y")

  # No bidirected edges (returns NULL when no spouses)
  expect_null(spouses(admg, "X"))
  expect_null(spouses(admg, "Y"))
  expect_null(spouses(admg, "Z"))
})

test_that("latent_project with multiple latents", {
  # DAG: L1 -> X, L1 -> Y, L2 -> Y, L2 -> Z, X -> Y, Y -> Z
  # Project out L1, L2
  #
  # Eliminate L1:
  # - Ch(L1) = {X, Y}
  # - Step 3: Add X <-> Y
  #
  # Eliminate L2:
  # - Ch(L2) = {Y, Z}
  # - Step 3: Add Y <-> Z
  #
  # Result: X -> Y -> Z (directed) + X <-> Y, Y <-> Z (bidirected)
  dag <- caugi(
    L1 %-->% X,
    L1 %-->% Y,
    L2 %-->% Y,
    L2 %-->% Z,
    X %-->% Y,
    Y %-->% Z,
    class = "DAG"
  )

  admg <- latent_project(dag, latents = c("L1", "L2"))

  expect_equal(admg@graph_class, "ADMG")
  expect_equal(length(admg), 3)
  expect_equal(nodes(admg)$name, c("X", "Y", "Z"))

  # Directed edges preserved
  expect_equal(parents(admg, "Y"), "X")
  expect_equal(parents(admg, "Z"), "Y")

  # - X <-> Y: added (children of L1 pair up)
  # - Y <-> Z: added (children of L2 pair up)
  expect_equal(spouses(admg, "X"), "Y")
  expect_equal(sort(spouses(admg, "Y")), c("X", "Z"))
  expect_equal(spouses(admg, "Z"), "Y")
})

test_that("latent_project with latent chain", {
  # DAG: L1 -> L2 -> X, L2 -> Y
  # Project out L1, L2
  # X and Y share latent ancestor L2, so X <-> Y
  dag <- caugi(
    L1 %-->% L2,
    L2 %-->% X,
    L2 %-->% Y,
    class = "DAG"
  )

  admg <- latent_project(dag, latents = c("L1", "L2"))

  expect_equal(admg@graph_class, "ADMG")
  expect_equal(length(admg), 2)
  expect_equal(nodes(admg)$name, c("X", "Y"))

  # No directed edges between X and Y (returns NULL when no children)
  expect_null(children(admg, "X"))
  expect_null(children(admg, "Y"))

  # X <-> Y from shared ancestor L2
  expect_equal(spouses(admg, "X"), "Y")
  expect_equal(spouses(admg, "Y"), "X")
})

test_that("latent_project all nodes latent returns empty ADMG", {
  dag <- caugi(
    L1 %-->% L2,
    class = "DAG"
  )

  admg <- latent_project(dag, latents = c("L1", "L2"))

  expect_equal(admg@graph_class, "ADMG")
  expect_equal(length(admg), 0)
})

test_that("latent_project no shared latent ancestors (no bidirected edges)", {
  # DAG: U -> X, V -> Y (no shared ancestor)
  # Project out U, V
  # Result: X, Y (isolated nodes, no edges)
  dag <- caugi(
    U %-->% X,
    V %-->% Y,
    class = "DAG"
  )

  admg <- latent_project(dag, latents = c("U", "V"))

  expect_equal(admg@graph_class, "ADMG")
  expect_equal(length(admg), 2)

  # No edges at all (returns NULL when no children/spouses)
  expect_null(children(admg, "X"))
  expect_null(children(admg, "Y"))
  expect_null(spouses(admg, "X"))
  expect_null(spouses(admg, "Y"))
})

test_that("latent_project fails on non-DAGs", {
  cg_pdag <- caugi(
    A %-->% B,
    B %---% C,
    class = "PDAG"
  )

  expect_error(
    latent_project(cg_pdag, latents = "A"),
    "latent_project\\(\\) can only be applied to DAGs\\."
  )

  cg_ug <- caugi(
    A %---% B,
    class = "UG"
  )

  expect_error(
    latent_project(cg_ug, latents = "A"),
    "latent_project\\(\\) can only be applied to DAGs\\."
  )
})

test_that("latent_project fails with unknown latent names", {
  dag <- caugi(
    X %-->% Y,
    class = "DAG"
  )

  expect_error(
    latent_project(dag, latents = "Z"),
    "Unknown latent node\\(s\\): Z"
  )

  expect_error(
    latent_project(dag, latents = c("X", "Z", "W")),
    "Unknown latent node\\(s\\): Z, W"
  )
})

test_that("latent_project fails with non-character latents", {
  dag <- caugi(
    X %-->% Y,
    class = "DAG"
  )

  expect_error(
    latent_project(dag, latents = 1),
    "`latents` must be a character vector"
  )
})

# ──────────────────────────────────────────────────────────────────────────────
# ───────────────────────────────── Exogenize ──────────────────────────────────
# ──────────────────────────────────────────────────────────────────────────────

test_that("exogenize works", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    class = "DAG"
  )
  exogenized_cg <- exogenize(cg, nodes = "B")
  cg_expected <- caugi(
    A %-->% C,
    B %-->% C,
    class = "DAG"
  )
  # Compare edges irrespective of order
  actual <- edges(exogenized_cg)[order(from, to, edge)]
  expected <- edges(cg_expected)[order(from, to, edge)]
  expect_equal(actual, expected)
  expect_setequal(nodes(exogenized_cg)$name, nodes(cg_expected)$name)
})

test_that("exogenize fails with nodes not in graph", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    class = "DAG"
  )
  expect_error(exogenize(cg, nodes = "D"), "Node D not in graph.")
})

test_that("exogenize fails with non-character nodes", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    class = "DAG"
  )
  expect_error(
    exogenize(cg, nodes = 1),
    "`nodes` must be a non-empty character vector of node names."
  )
  expect_error(
    exogenize(cg, nodes = character(0)),
    "`nodes` must be a non-empty character vector of node names."
  )
})

test_that("exogenize agrees with exogenous query", {
  cg <- generate_graph(n = 100, m = 10, class = "DAG")
  exogenous_nodes <- exogenous(cg)
  new_exogenous_nodes <- sample(nodes(cg)$name, size = 10)
  exogenized_cg <- exogenize(cg, nodes = new_exogenous_nodes)
  expect_setequal(
    exogenous(exogenized_cg),
    union(exogenous_nodes, new_exogenous_nodes)
  )
})

test_that("exogenize supports multiple nodes", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    C %-->% D,
    class = "DAG"
  )
  exogenized_cg <- exogenize(cg, nodes = c("B", "C"))

  expect_null(parents(exogenized_cg, "B"))
  expect_null(parents(exogenized_cg, "C"))
  expect_setequal(parents(exogenized_cg, "D"), c("B", "C", "A"))
})

test_that("exogenize handles duplicate node inputs", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    class = "DAG"
  )

  exo_once <- exogenize(cg, nodes = "B")
  exo_twice <- exogenize(cg, nodes = c("B", "B"))

  actual_once <- edges(exo_once)[order(from, to, edge)]
  actual_twice <- edges(exo_twice)[order(from, to, edge)]
  expect_equal(actual_once, actual_twice)
})

test_that("exogenize fails with non-simple graphs", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    class = "UNKNOWN"
  )
  expect_error(
    exogenize(cg, nodes = "B"),
    "`cg` must be a DAG"
  )
})


test_that("Marginalization and conditioning work", {
  ## Examples from Figure 10 of Richardson and Spirtes 2002

  mg <- caugi(U %-->% X + Y, A %-->% X, B %-->% Y, class = "DAG")

  mgmu <- condition_marginalize(mg, marg_vars = "U")
  resmu <- caugi(A %-->% X, B %-->% Y, X %<->% Y, class = "AG")

  same_graphs <- function(cg1, cg2) {
    if (!cg1@graph_class == cg2@graph_class) {
      FALSE
    }
    if (!setequal(nodes(cg1)$name, nodes(cg2)$name)) {
      FALSE
    }

    shd(cg1, cg2) == 0
  }

  expect_true(same_graphs(mgmu, resmu))

  rescu <- caugi(A %-->% X, B %-->% Y)
  mgcu <- condition_marginalize(mg, cond_vars = "U")

  expect_true(same_graphs(rescu, mgcu))

  expect_error(condition_marginalize(mg, cond_vars = NULL, marg_vars = NULL))

  # Figure 11

  f11 <- caugi(
    A %-->% L1,
    L1 %-->% B,
    L2 %-->% B + C,
    B %-->% S,
    S %-->% D,
    D %-->% C
  )

  f11.cS <- condition_marginalize(f11, cond_vars = "S")
  f11.mL1L2 <- condition_marginalize(f11, marg_vars = c("L1", "L2"))
  f11.cSmL1L2 <- condition_marginalize(
    f11,
    marg_vars = c("L1", "L2"),
    cond_vars = "S"
  )

  f11.ii <- caugi(
    A %---% L1,
    L1 %---% L2 + B,
    L2 %---% B,
    L2 %-->% C,
    D %-->% C,
    class = "AG"
  )
  f11.iii <- caugi(
    A %-->% B + C,
    B %-->% S + C,
    S %-->% D,
    D %-->% C,
    class = "AG"
  )
  f11.iv <- caugi(A %---% B, A %-->% C, B %-->% C, D %-->% C, class = "AG")

  expect_true(same_graphs(f11.cS, f11.ii))
  expect_true(same_graphs(f11.mL1L2, f11.iii))
  expect_true(same_graphs(f11.cSmL1L2, f11.iv))

  ## order doesn't matter

  expect_true(same_graphs(
    condition_marginalize(f11.cS, marg_vars = c("L1", "L2")),
    f11.cSmL1L2
  ))
  expect_true(same_graphs(
    condition_marginalize(f11.mL1L2, cond_vars = "S"),
    f11.cSmL1L2
  ))

  ## error checking

  emg <- caugi(U %-->% X + Y, A %-->% X, B %-->% Y, K %<->% Y, class = "AG")

  expect_no_error(condition_marginalize(emg, cond_vars = "U"))
  expect_error(condition_marginalize(
    f11,
    cond_vars = c("S", "L1"),
    marg_vars = c("L1")
  ))
})

# ──────────────────────────────────────────────────────────────────────────────
# ───────────────────────────────── DAG from PDAG ──────────────────────────────
# ──────────────────────────────────────────────────────────────────────────────

test_that("dag_from_pdag converts a valid PDAG to a DAG", {
  PDAG <- caugi(
    A %---% B,
    B %---% C,
    class = "PDAG"
  )
  DAG <- dag_from_pdag(PDAG)
  edges_DAG <- edges(DAG)
  expect_true(all(edges_DAG$edge == "-->"))
})

test_that("dag_from_pdag handles a sink with multiple undirected neighbors (#298)", {
  pdag <- caugi(
    V10 %---% V2 + V7,
    V2 %-->% V1,
    V2 %---% V3 + V7,
    V3 %-->% V5,
    V3 %---% V7,
    V4 %-->% V8,
    V5 %-->% V1,
    V6 %-->% V8,
    V6 %---% V9,
    V7 %-->% V4,
    V9 %-->% V1 + V4 + V5,
    class = "PDAG"
  )

  dag <- dag_from_pdag(pdag)

  expected <- caugi(
    V2 %-->% V1 + V10,
    V3 %-->% V2 + V5,
    V4 %-->% V8,
    V5 %-->% V1,
    V6 %-->% V8,
    V7 %-->% V10 + V2 + V3 + V4,
    V9 %-->% V1 + V4 + V5 + V6,
    class = "DAG"
  )

  expect_equal(dag, expected)
})


test_that("dag_from_pdag errors on non-PDAG input", {
  cg <- caugi(
    A %-->% B,
    B %-->% C,
    class = "DAG"
  )
  expect_error(
    dag_from_pdag(cg),
    "Input must be a caugi PDAG/MPDAG/CPDAG graph"
  )
})

test_that("dag_from_pdag errors if PDAG cannot be extended to a DAG", {
  PDAG <- caugi(
    A %---% B,
    A %---% D,
    B %---% C,
    C %---% D,
    class = "PDAG"
  )
  expect_error(dag_from_pdag(PDAG), "PDAG cannot be extended to a DAG")
})

test_that("dag_from_pdag accepts the deprecated `PDAG` argument with a warning", {
  pdag <- caugi(
    A %---% B,
    B %---% C,
    class = "PDAG"
  )
  sig <- function(cg) {
    e <- edges(cg)
    sort(paste0(as.character(e$from), "->", as.character(e$to)))
  }
  expect_warning(
    dag <- dag_from_pdag(PDAG = pdag),
    "`PDAG` argument.*deprecated"
  )
  expect_equal(dag@graph_class, "DAG")
  # The deprecated alias yields the same result as the current argument.
  expect_setequal(sig(dag), sig(dag_from_pdag(pdag)))
})

# ──────────────────────────────────────────────────────────────────────────────
# ──────────────────────────────── Meek closure ───────────────────────────────
# ──────────────────────────────────────────────────────────────────────────────

test_that("meek_closure orients compelled edges and returns an MPDAG", {
  check_oriented <- function(g, from, to) {
    expect_true(to %in% children(g, from))
    if (is_dag(g)) {
      return()
    } else {
      und <- neighbors(g, from, mode = "undirected")
      if (is.null(und)) {
        und <- character(0)
      }
      expect_false(to %in% und)
    }
  }

  # R1: A -> B, C -> B, B -- D, and C !~ D  =>  B -> D
  g_r1 <- caugi(
    A %-->% B,
    C %-->% B,
    B %---% D,
    A %---% D,
    class = "PDAG"
  )
  c_r1 <- meek_closure(g_r1)
  check_oriented(c_r1, "B", "D")
  expect_true(is_mpdag(c_r1))

  # R2: A -- B and A -> C -> B  =>  A -> B
  g_r2 <- caugi(
    A %---% B,
    A %-->% C,
    C %-->% B,
    class = "PDAG"
  )
  c_r2 <- meek_closure(g_r2)
  check_oriented(c_r2, "A", "B")
  expect_true(is_mpdag(c_r2))

  # R3: A -- B, C -> B, D -> B, C !~ D, A -- C, A -- D  =>  A -> B
  g_r3 <- caugi(
    A %---% B,
    C %-->% B,
    D %-->% B,
    A %---% C,
    A %---% D,
    class = "PDAG"
  )
  c_r3 <- meek_closure(g_r3)
  check_oriented(c_r3, "A", "B")
  expect_true(is_mpdag(c_r3))

  # R4: A -- B and A => B through A -> C -> D -> B  =>  A -> B
  g_r4 <- caugi(
    A %---% B,
    A %-->% C,
    C %-->% D,
    A %---% D,
    D %-->% B,
    class = "PDAG"
  )
  c_r4 <- meek_closure(g_r4)
  check_oriented(c_r4, "A", "B")
  expect_true(is_mpdag(c_r4))
})

test_that("meek_closure is idempotent and preserves node names", {
  g <- caugi(
    B %---% A,
    B %---% D,
    C %-->% E,
    B %-->% E,
    D %-->% F,
    E %-->% F,
    class = "PDAG"
  )

  c1 <- meek_closure(g)
  c2 <- meek_closure(c1)

  expect_equal(c1@graph_class, "MPDAG")
  expect_equal(nodes(c1)$name, nodes(g)$name)
  expect_equal(edges(c1), edges(c2))
})

test_that("meek_closure errors on non-PDAG-compatible graphs", {
  g <- caugi(A %o->% B, class = "UNKNOWN")
  expect_error(
    meek_closure(g),
    "meek_closure\\(\\) can only be applied to PDAGs\\."
  )
})

test_that("meek_closure matches causal-learn style multi-rule regression", {
  # Skeleton: A-B-C, A->D<-C, B-D, D-E, C-E
  # Meek progression:
  #   R1: D -> E
  #   R2: C -> E
  #   R3: B -> D
  g <- caugi(
    A %---% B,
    B %---% C,
    A %-->% D,
    C %-->% D,
    B %---% D,
    D %---% E,
    C %---% E,
    class = "PDAG"
  )

  g_closed <- meek_closure(g)

  expected <- caugi(
    A %---% B,
    B %---% C,
    A %-->% D,
    B %-->% D,
    C %-->% D,
    D %-->% E,
    C %-->% E,
    class = "PDAG"
  )

  norm_edges <- function(x) edges(x)[order(from, to, edge)]
  expect_equal(norm_edges(g_closed), norm_edges(expected))

  # Invariants inspired by CPDAG/Meek references.
  expect_equal(
    norm_edges(skeleton(g_closed)),
    norm_edges(skeleton(g))
  )
  expect_true(is_acyclic(g_closed, force_check = TRUE))
  expect_true(is_mpdag(g_closed))
})

test_that("dag_from_pdag preserves directed edges in mixed extension cases", {
  pdag <- caugi(
    A %-->% B,
    C %-->% B,
    C %---% D,
    D %---% A,
    class = "PDAG"
  )

  dag <- dag_from_pdag(pdag)
  ed <- edges(dag)
  has_dir <- function(from, to) {
    any(ed$from == from & ed$to == to & ed$edge == "-->")
  }

  expect_true(is_dag(dag))
  expect_true(has_dir("A", "B"))
  expect_true(has_dir("C", "B"))
  expect_true(xor(has_dir("A", "D"), has_dir("D", "A")))
  expect_true(xor(has_dir("C", "D"), has_dir("D", "C")))
  expect_false(has_dir("A", "D") && has_dir("C", "D"))
  expect_equal(nrow(ed), 4L)
})

test_that("dag_from_pdag adds no v-structure the PDAG did not have (#298)", {
  # `A --> B <-- C` with `D --- B`. Orienting `D --> B` would add two
  # v-structures, since D is adjacent to neither A nor C, so the only
  # consistent extension is `B --> D`. Dor-Tarsi gets this from condition
  # (b) quantifying over *all* neighbours of the candidate sink, not just
  # its undirected ones -- B has the parents A and C, which D is not
  # adjacent to, so B is not a sink here and D is.
  pdag <- caugi(
    A %-->% B,
    C %-->% B,
    B %---% D,
    class = "PDAG"
  )

  dag <- dag_from_pdag(pdag)
  ed <- as.data.frame(edges(dag))
  expect_true(is_dag(dag))
  expect_equal(nrow(ed), 3L)
  expect_true(any(ed$from == "B" & ed$to == "D"))
  # The v-structures are exactly the ones the PDAG had.
  expect_setequal(parents(dag, "B"), c("A", "C"))
  expect_setequal(parents(dag, "D"), "B")
})

test_that("dag_from_pdag orients each undirected edge exactly once", {
  pdag <- caugi(
    A %-->% C,
    B %-->% C,
    A %---% D,
    class = "PDAG"
  )

  dag <- dag_from_pdag(pdag)
  ed <- edges(dag)
  has_dir <- function(from, to) {
    any(ed$from == from & ed$to == to & ed$edge == "-->")
  }

  expect_true(has_dir("A", "C"))
  expect_true(has_dir("B", "C"))
  expect_true(xor(has_dir("A", "D"), has_dir("D", "A")))
  expect_false(any(ed$edge == "---"))
  expect_equal(nrow(ed), 3L)
})

test_that("pgmpy Meek fixtures: meek_closure on PDAG applies Meek closure", {
  cpdag_from_pdag <- function(pdag) {
    cp_session <- rs_meek_closure(pdag@session)
    .session_to_caugi(cp_session, node_names = nodes(pdag)$name)
  }
  edge_set <- function(g) {
    ed <- edges(g)
    sort(paste(ed$from, ed$edge, ed$to))
  }

  # Case 1 (pgmpy test_pdag_to_cpdag): A->B and B-C => B->C
  pdag1 <- caugi(A %-->% B, B %---% C, class = "PDAG")
  cpdag1 <- cpdag_from_pdag(pdag1)
  expect_setequal(
    edge_set(cpdag1),
    edge_set(caugi(A %-->% B, B %-->% C, class = "PDAG"))
  )

  # Case 2: orientation propagates along B-C-D.
  pdag2 <- caugi(A %-->% B, B %---% C, C %---% D, class = "PDAG")
  cpdag2 <- cpdag_from_pdag(pdag2)
  expect_setequal(
    edge_set(cpdag2),
    edge_set(caugi(A %-->% B, B %-->% C, C %-->% D, class = "PDAG"))
  )

  # Case 3 (pgmpy fixture): keep B-C undirected.
  pdag3 <- caugi(A %-->% B, D %-->% C, B %---% C, class = "PDAG")
  cpdag3 <- cpdag_from_pdag(pdag3)
  expect_setequal(
    edge_set(cpdag3),
    edge_set(pdag3)
  )

  # Case 4 (pgmpy fixture): extra parent evidence orients B->C.
  pdag4 <- caugi(
    A %-->% B,
    D %-->% C,
    D %-->% B,
    B %---% C,
    class = "PDAG"
  )
  cpdag4 <- cpdag_from_pdag(pdag4)
  expect_setequal(
    edge_set(cpdag4),
    edge_set(caugi(A %-->% B, D %-->% C, D %-->% B, B %-->% C, class = "PDAG"))
  )

  # Case 5 (rule-2 style): A->B->C and A-C => A->C.
  pdag5 <- caugi(A %-->% B, B %-->% C, A %---% C, class = "PDAG")
  cpdag5 <- cpdag_from_pdag(pdag5)
  expect_setequal(
    edge_set(cpdag5),
    edge_set(caugi(A %-->% B, B %-->% C, A %-->% C, class = "PDAG"))
  )
})

test_that("condition_marginalize and helper branches are covered", {
  cg_ug <- caugi(A %o->% B, class = "UNKNOWN")
  expect_error(
    condition_marginalize(cg_ug, cond_vars = "A"),
    "`cg` must be an AG for `condition_marginalize\\(\\)`"
  )

  cg <- caugi(A %-->% B, class = "DAG")
  out_trivial <- condition_marginalize(cg, marg_vars = "B")
  expect_equal(out_trivial@graph_class, "AG")
  expect_setequal(out_trivial@nodes$name, "A")

  # Force non-list anteriors branch (line coverage via expected error).
  cg2 <- caugi(A %-->% B %-->% C, class = "DAG")
  expect_error(
    testthat::with_mocked_bindings(
      condition_marginalize(cg2, cond_vars = "B"),
      anteriors = function(...) "A",
      .package = "caugi"
    ),
    "must be the same length"
  )

  expect_true(
    caugi:::.not_m_separated_for_all_subsets(
      cg = cg2,
      node_a = "A",
      node_b = "C",
      other_nodes = character(0),
      cond_vars = character(0)
    )
  )

  expect_false(
    caugi:::.not_m_separated_for_all_subsets(
      cg = cg2,
      node_a = "A",
      node_b = "C",
      other_nodes = "B",
      cond_vars = character(0)
    )
  )

  edge_rev <- caugi:::.edge_type_from_anteriors(
    node_a = "A",
    node_b = "B",
    cond_vars = character(0),
    anteriors_list = list(A = "B", B = NULL)
  )
  expect_identical(edge_rev$from, "B")
  expect_identical(edge_rev$edge, "-->")
  expect_identical(edge_rev$to, "A")
})

test_that("not_m_separated helper matches R reference across subset scenarios", {
  helper_reference <- function(cg, node_a, node_b, other_nodes, cond_vars) {
    n_other <- length(other_nodes)
    subsets <- if (n_other == 0L) {
      list(NULL)
    } else {
      c(
        list(other_nodes),
        if (n_other > 1L) {
          unlist(
            lapply(
              seq_len(n_other - 1L),
              function(k) combn(other_nodes, n_other - k, simplify = FALSE)
            ),
            recursive = FALSE
          )
        },
        list(NULL)
      )
    }

    for (subset in subsets) {
      z <- c(cond_vars, subset)
      if (length(z) == 0L) {
        z <- NULL
      }
      if (m_separated(cg, X = node_a, Y = node_b, Z = z)) {
        return(FALSE)
      }
    }
    TRUE
  }

  # 1) Separation found only for full subset.
  # A <- B -> C and A <- D -> C: both B and D are required to block all paths.
  g_full <- caugi(B %-->% A + C, D %-->% A + C, class = "DAG")
  expect_true(m_separated(g_full, X = "A", Y = "C", Z = c("B", "D")))
  expect_false(m_separated(g_full, X = "A", Y = "C", Z = "B"))
  expect_false(m_separated(g_full, X = "A", Y = "C", Z = "D"))

  rust_full <- caugi:::.not_m_separated_for_all_subsets(
    g_full,
    node_a = "A",
    node_b = "C",
    other_nodes = c("B", "D"),
    cond_vars = character(0)
  )
  ref_full <- helper_reference(
    g_full,
    node_a = "A",
    node_b = "C",
    other_nodes = c("B", "D"),
    cond_vars = character(0)
  )
  expect_identical(rust_full, ref_full)

  # 2) Separation found only for an intermediate subset size.
  # A <- B -> C is blocked by conditioning on B, but conditioning on D opens A->E<-C via E->D.
  g_mid <- caugi(B %-->% A + C, A %-->% E, C %-->% E, E %-->% D, class = "DAG")
  expect_true(m_separated(g_mid, X = "A", Y = "C", Z = "B"))
  expect_false(m_separated(g_mid, X = "A", Y = "C", Z = c("B", "D")))

  rust_mid <- caugi:::.not_m_separated_for_all_subsets(
    g_mid,
    node_a = "A",
    node_b = "C",
    other_nodes = c("B", "D"),
    cond_vars = character(0)
  )
  ref_mid <- helper_reference(
    g_mid,
    node_a = "A",
    node_b = "C",
    other_nodes = c("B", "D"),
    cond_vars = character(0)
  )
  expect_identical(rust_mid, ref_mid)

  # 3) Separation found for cond_vars alone.
  g_cond <- caugi(S %-->% A + C, T %-->% A, class = "DAG")
  expect_true(m_separated(g_cond, X = "A", Y = "C", Z = "S"))

  rust_cond <- caugi:::.not_m_separated_for_all_subsets(
    g_cond,
    node_a = "A",
    node_b = "C",
    other_nodes = "T",
    cond_vars = "S"
  )
  ref_cond <- helper_reference(
    g_cond,
    node_a = "A",
    node_b = "C",
    other_nodes = "T",
    cond_vars = "S"
  )
  expect_identical(rust_cond, ref_cond)

  # Additional parity checks on small random DAGs (all subsets feasible).
  set.seed(123)
  for (i in seq_len(20L)) {
    g <- generate_graph(n = 8, m = 8, class = "DAG")
    nn <- nodes(g)$name
    pair <- sample(nn, size = 2, replace = FALSE)
    rest <- setdiff(nn, pair)
    other <- sample(rest, size = sample.int(3L, 1), replace = FALSE)
    cond_pool <- setdiff(rest, other)
    cond <- if (length(cond_pool) > 0L) {
      sample(
        cond_pool,
        size = sample.int(min(2L, length(cond_pool)), 1) - 1L,
        replace = FALSE
      )
    } else {
      character(0)
    }

    rust <- caugi:::.not_m_separated_for_all_subsets(
      g,
      node_a = pair[[1]],
      node_b = pair[[2]],
      other_nodes = other,
      cond_vars = cond
    )
    ref <- helper_reference(
      g,
      node_a = pair[[1]],
      node_b = pair[[2]],
      other_nodes = other,
      cond_vars = cond
    )
    expect_identical(rust, ref)
  }
})

# ──────────────────────────────────────────────────────────────────────────────
# ───────────────────────── Enumerate / count DAGs in MEC ──────────────────────
# ──────────────────────────────────────────────────────────────────────────────

# Canonical signature: edge set as a sorted character vector "from->to".
.dag_edge_signature <- function(cg) {
  e <- edges(cg)
  if (nrow(e) == 0L) {
    return(character(0))
  }
  sort(paste0(as.character(e$from), "->", as.character(e$to)))
}

test_that("enumerate_dags returns 3 DAGs for an A--B--C chain", {
  pdag <- caugi(
    A %---% B,
    B %---% C,
    class = "PDAG"
  )

  dags <- enumerate_dags(pdag)
  expect_length(dags, 3L)
  expect_equal(count_dags(pdag), 3)

  for (d in dags) {
    expect_equal(d@graph_class, "DAG")
    expect_true(all(edges(d)$edge == "-->"))
    sig <- .dag_edge_signature(d)
    # No DAG may contain the v-structure A->B<-C.
    expect_false(all(c("A->B", "C->B") %in% sig))
  }

  sigs <- vapply(
    dags,
    function(d) paste(.dag_edge_signature(d), collapse = "|"),
    character(1)
  )
  expect_length(unique(sigs), 3L)
})

test_that("enumerate_dags returns the singleton for a v-structure CPDAG", {
  pdag <- caugi(
    A %-->% C,
    B %-->% C,
    class = "PDAG"
  )

  dags <- enumerate_dags(pdag)
  expect_length(dags, 1L)
  expect_equal(count_dags(pdag), 1)
  expect_setequal(.dag_edge_signature(dags[[1]]), c("A->C", "B->C"))
})

test_that("enumerate_dags returns 6 DAGs for an undirected triangle", {
  pdag <- caugi(
    A %---% B,
    A %---% C,
    B %---% C,
    class = "PDAG"
  )

  dags <- enumerate_dags(pdag)
  expect_length(dags, 6L)
  expect_equal(count_dags(pdag), 6)

  sigs <- vapply(
    dags,
    function(d) paste(.dag_edge_signature(d), collapse = "|"),
    character(1)
  )
  expect_length(unique(sigs), 6L)
})

test_that("enumerate_dags propagates background-knowledge orientations via Meek", {
  # A->C<-B (fixed v-structure) plus C--D: Meek R1 forces C->D, so MEC = {one DAG}.
  pdag <- caugi(
    A %-->% C,
    B %-->% C,
    C %---% D,
    class = "PDAG"
  )

  dags <- enumerate_dags(pdag)
  expect_length(dags, 1L)
  expect_setequal(
    .dag_edge_signature(dags[[1]]),
    c("A->C", "B->C", "C->D")
  )
})

test_that("enumerate_dags includes the dag_from_pdag extension", {
  pdag <- caugi(
    A %---% B,
    B %---% C,
    class = "PDAG"
  )
  one <- dag_from_pdag(pdag)
  all <- enumerate_dags(pdag)
  target <- paste(.dag_edge_signature(one), collapse = "|")
  sigs <- vapply(
    all,
    function(d) paste(.dag_edge_signature(d), collapse = "|"),
    character(1)
  )
  expect_true(target %in% sigs)
})

test_that("count_dags agrees with length(enumerate_dags()) on mixed fixtures", {
  fixtures <- list(
    caugi(A %---% B, class = "PDAG"),
    caugi(A %---% B, B %---% C, class = "PDAG"),
    caugi(A %---% B, A %---% C, B %---% C, class = "PDAG"),
    caugi(A %-->% C, B %-->% C, C %---% D, class = "PDAG"),
    caugi(A %---% B, C %---% D, D %---% E, class = "PDAG")
  )
  for (p in fixtures) {
    expect_equal(count_dags(p), length(enumerate_dags(p)))
  }
})

test_that("enumerate_dags and count_dags reject non-PDAG input", {
  admg <- caugi(A %<->% B, class = "ADMG")
  expect_error(enumerate_dags(admg), "can only be applied to PDAGs")
  expect_error(count_dags(admg), "can only be applied to PDAGs")
})

test_that("enumerate_dags rejects v-structures introduced by Meek R2/R3/R4 propagation", {
  # Trial-8 fixture from random-CPDAG cross-validation against pcalg.
  # Skeleton is dense enough that R2/R3/R4 propagation during enumeration
  # can introduce new v-structures (3->2<-5 with 3 not adjacent to 5);
  # those branches must be rejected. pcalg::pdag2allDags returns 0 here.
  pdag <- caugi(
    `7` %-->% `1`,
    `7` %-->% `3`,
    `7` %-->% `5`,
    `7` %-->% `6`,
    `1` %---% `2`,
    `1` %---% `3`,
    `1` %---% `4`,
    `1` %---% `5`,
    `1` %---% `6`,
    `2` %---% `3`,
    `2` %---% `4`,
    `2` %---% `5`,
    `2` %---% `6`,
    `3` %---% `4`,
    `4` %---% `5`,
    `4` %---% `6`,
    `5` %---% `6`,
    nodes = c("1", "2", "3", "4", "5", "6", "7"),
    class = "PDAG"
  )
  expect_equal(count_dags(pdag), 0)
  expect_length(enumerate_dags(pdag), 0L)
})

test_that("enumerate_dags matches pcalg::pdag2allDags on random CPDAGs", {
  skip_if_not_installed("pcalg")

  # pcalg amat.cpdag convention: m[j, i] = 1, m[i, j] = 0  iff  i -> j.
  # caugi's as_adjacency() uses the opposite convention, so transpose.
  pcalg_sig <- function(row, p, nn) {
    m <- matrix(row, p, p, byrow = TRUE, dimnames = list(nn, nn))
    es <- character(0)
    for (i in seq_len(p)) {
      for (j in seq_len(p)) {
        if (i != j && m[j, i] == 1 && m[i, j] == 0) {
          es <- c(es, paste0(nn[i], "->", nn[j]))
        }
      }
    }
    paste(sort(es), collapse = "|")
  }
  caugi_sig <- function(cg) {
    e <- edges(cg)
    if (nrow(e) == 0L) {
      return("")
    }
    paste(
      sort(paste0(as.character(e$from), "->", as.character(e$to))),
      collapse = "|"
    )
  }

  for (seed in 1:8) {
    cg <- generate_graph(n = 5, p = 0.4, class = "CPDAG", seed = seed)
    adj <- t(as_adjacency(cg))
    nn <- rownames(adj)

    caugi_dags <- sort(vapply(enumerate_dags(cg), caugi_sig, character(1)))
    pc <- pcalg::pdag2allDags(adj)
    pcalg_dags <- if (
      is.null(pc$dags) || !is.matrix(pc$dags) || nrow(pc$dags) == 0L
    ) {
      character(0)
    } else {
      sort(vapply(
        seq_len(nrow(pc$dags)),
        function(i) pcalg_sig(pc$dags[i, ], length(nn), nn),
        character(1)
      ))
    }
    expect_identical(caugi_dags, pcalg_dags)
    expect_equal(count_dags(cg), length(pcalg_dags))
  }
})