demo.R

# The path of "R-proj"
proj_path <- getwd()
src_path <- paste(proj_path, '/src', sep = "")
base_res_path <- paste(proj_path, '/result', sep = "") # Base results directory
data_path <- paste(proj_path, '/data', sep = "")

source(file.path(src_path, 'CITest.R'))
source(file.path(src_path, 'Evaluation.R')) # This contains test_adj
source(file.path(src_path, 'MissingValuePC.R'))
source(file.path(src_path, 'SyntheticDataGeneration.R'))

# --- Define Parameters for Data Generation ---
num_var_val <- 20
num_sample_val <- 100
mode_val <- "mnar"
num_extra_e_val <- 5
num_m_val <- 10
seed_val <- 777

folder_name <- paste0(
  "s", num_sample_val,
  "_", mode_val,
  "_v", num_var_val,
  "_e", num_extra_e_val,
  "_m", num_m_val,
  "_seed", seed_val
)
current_res_path <- file.path(base_res_path, folder_name) # Path for current run's results

if (!dir.exists(current_res_path)) {
  dir.create(current_res_path, recursive = TRUE)
  cat("Created results directory for this run:", current_res_path, "\n")
} else {
  cat("Directory for this run already exists:", current_res_path, "\n")
}

# --- Synthethic data generation ---
cat("Generating synthetic data...\n")
gen_result_list <- gen_data(
  num_sample = num_sample_val,
  mode = mode_val,
  num_var = num_var_val,
  num_extra_e = num_extra_e_val,
  num_m = num_m_val,
  seed = seed_val
)

data_complete <- gen_result_list$data_complete
data_m <- gen_result_list$data_m
myCPDAG <- gen_result_list$ground_truth$cpdag # Ground truth graph object (e.g., graphNEL)

suffStat_m <- list(data = data_m)
suffStat <- list(C = cor(data_complete), n = num_sample_val)

cat("Running causal discovery algorithms...\n")

# PC on the complete dataset
res_com_pc <- pc(suffStat,
                 gaussCItest,
                 alpha = 0.01,
                 p = num_var_val)

# Test-wise deletion PC
res_tw <- pc(suffStat_m,
             gaussCItest.td,
             alpha = 0.01,
             p = num_var_val)

# MVPC with the PermC correction method
res.mvpc.permc <- mvpc(suffStat_m,
                       gaussCItest.td,
                       gaussCItest.permc,
                       alpha = 0.01,
                       p = num_var_val)

# MVPC with the DRW correction method
res.mvpc.drw <- mvpc(suffStat_m,
                     gaussCItest.td,
                     gaussCItest.drw,
                     alpha = 0.01,
                     p = num_var_val)

# --- Convert Graph Objects to Adjacency Matrices ---
cat("Converting graphs to adjacency matrices and saving them...\n")

# 1. Ground Truth CPDAG
if (is(myCPDAG, "graphNEL")) {
  ground_truth_adj_matrix <- as(myCPDAG, "matrix")
} else if (is.matrix(myCPDAG)) {
  ground_truth_adj_matrix <- myCPDAG
} else {
  warning("myCPDAG type not directly recognized as graphNEL or matrix, attempting as(myCPDAG, 'matrix').")
  ground_truth_adj_matrix <- try(as(myCPDAG, "matrix"))
  if (inherits(ground_truth_adj_matrix, "try-error")) {
    stop("Failed to convert myCPDAG to matrix. Please check its structure.")
  }
}
write.csv(ground_truth_adj_matrix, file.path(current_res_path, "ground_truth_cpdag_adj.csv"), row.names = TRUE)
cat("Saved ground_truth_cpdag_adj.csv\n")

# 2. PC on complete data
res_com_pc_adj_matrix <- as(res_com_pc@graph, "matrix")
write.csv(res_com_pc_adj_matrix, file.path(current_res_path, "res_com_pc_adj.csv"), row.names = TRUE)
cat("Saved res_com_pc_adj.csv\n")

# 3. Test-wise deletion PC
res_tw_adj_matrix <- as(res_tw@graph, "matrix")
write.csv(res_tw_adj_matrix, file.path(current_res_path, "res_tw_adj.csv"), row.names = TRUE)
cat("Saved res_tw_adj.csv\n")

# 4. MVPC with PermC
# Assuming mvpc output has @graph slot like pcAlgo. Adjust if mvpc returns graph object differently.
res_mvpc_permc_adj_matrix <- as(res.mvpc.permc@graph, "matrix")
write.csv(res_mvpc_permc_adj_matrix, file.path(current_res_path, "res_mvpc_permc_adj.csv"), row.names = TRUE)
cat("Saved res_mvpc_permc_adj.csv\n")

# 5. MVPC with DRW
# Assuming mvpc output has @graph slot like pcAlgo. Adjust if mvpc returns graph object differently.
res_mvpc_drw_adj_matrix <- as(res.mvpc.drw@graph, "matrix")
write.csv(res_mvpc_drw_adj_matrix, file.path(current_res_path, "res_mvpc_drw_adj.csv"), row.names = TRUE)
cat("Saved res_mvpc_drw_adj.csv\n")

cat("All adjacency matrices saved.\n")


# --- Evaluation: SHD and F1 Score ---
cat("Calculating SHD and F1 scores...\n")

# Helper function to calculate F1 from precision and recall
calculate_f1_from_pr <- function(precision, recall) {
  if (is.na(precision) || is.na(recall) || (precision + recall) == 0) {
    return(0) # Or NA, depending on how you want to handle it
  }
  return(2 * (precision * recall) / (precision + recall))
}

methods_results <- list()

# 1. Complete Data PC
shd_com_pc <- shd(res_com_pc, myCPDAG)
pr_com_pc <- test_adj(myCPDAG, res_com_pc@graph) # test_adj takes graph objects
f1_com_pc <- calculate_f1_from_pr(pr_com_pc$precision, pr_com_pc$recall)
methods_results[["CompleteData_PC"]] <- list(shd = shd_com_pc, precision = pr_com_pc$precision, recall = pr_com_pc$recall, f1 = f1_com_pc)

# 2. Test-wise Deletion PC
shd_tw <- shd(res_tw, myCPDAG)
pr_tw <- test_adj(myCPDAG, res_tw@graph)
f1_tw <- calculate_f1_from_pr(pr_tw$precision, pr_tw$recall)
methods_results[["TestWise_PC"]] <- list(shd = shd_tw, precision = pr_tw$precision, recall = pr_tw$recall, f1 = f1_tw)

# 3. MVPC with PermC
shd_permc <- shd(res.mvpc.permc, myCPDAG)
pr_permc <- test_adj(myCPDAG, res.mvpc.permc@graph)
f1_permc <- calculate_f1_from_pr(pr_permc$precision, pr_permc$recall)
methods_results[["MVPC_PermC"]] <- list(shd = shd_permc, precision = pr_permc$precision, recall = pr_permc$recall, f1 = f1_permc)

# 4. MVPC with DRW
shd_drw <- shd(res.mvpc.drw, myCPDAG)
pr_drw <- test_adj(myCPDAG, res.mvpc.drw@graph)
f1_drw <- calculate_f1_from_pr(pr_drw$precision, pr_drw$recall)
methods_results[["MVPC_DRW"]] <- list(shd = shd_drw, precision = pr_drw$precision, recall = pr_drw$recall, f1 = f1_drw)

evaluation_df_methods <- do.call(rbind, lapply(names(methods_results), function(name) {
  data.frame(
    Method = name,
    SHD = methods_results[[name]]$shd,
    Precision = methods_results[[name]]$precision,
    Recall = methods_results[[name]]$recall,
    F1_Score = methods_results[[name]]$f1
  )
}))

evaluation_df <- data.frame(
  Num_Samples = num_sample_val,
  Mode = mode_val,
  Num_Variables = num_var_val,
  Num_Extra_Edges = num_extra_e_val,
  Num_Missing_Vars_Affected = num_m_val, # Clarified name for num_m
  Seed = seed_val
)

evaluation_df <- cbind(evaluation_df, evaluation_df_methods)

write.csv(evaluation_df, file.path(current_res_path, "evaluation_metrics.csv"), row.names = FALSE)
cat("Evaluation metrics (including dataset parameters, SHD, Precision, Recall, F1) saved to evaluation_metrics.csv in", current_res_path, "\n")

cat("Demo script finished.\n")