Created
May 27, 2025 11:06
-
-
Save sachsmc/7c53a4a487303a13f58688ded5d3d07e to your computer and use it in GitHub Desktop.
Using mplrs to do the vertex enumeration in parallel from causaloptim
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| library(causaloptim) | |
| library(rcdd) | |
| write_problem <- function(opt, obj) { | |
| # The Primal LP | |
| c0 <- obj$c0 | |
| n <- nrow(c0) | |
| A_e <- obj$R | |
| p <- obj$parameters | |
| m_e <- nrow(A_e) | |
| A_l <- obj$iqR | |
| if (is.null(A_l)) A_l <- matrix(data = 0, nrow = 0, ncol = n) | |
| b_l <- obj$iqb | |
| if (is.null(b_l)) b_l <- matrix(data = 0, nrow = 0, ncol = 1) | |
| if (!is.numeric(A_l) || !is.numeric(b_l)) stop("Inequality entries must be numeric.") | |
| if (ncol(A_l) != n) stop("Dimension mismatch of inequality constaints.") | |
| m_l <- nrow(A_l) | |
| if (nrow(b_l) != m_l) stop("Dimension mismatch of inequality constaints.") | |
| m <- m_l + m_e | |
| # The Dual LP | |
| a1 <- rbind(cbind(t(A_l), t(A_e)), | |
| cbind(diag(x = 1, nrow = m_l, ncol = m_l), matrix(data = 0, nrow = m_l, ncol = m_e))) | |
| b1 <- rbind(c0, | |
| matrix(data = 0, nrow = m_l, ncol = 1)) | |
| if (opt == "max") { | |
| a1 <- -a1 | |
| b1 <- -b1 | |
| } | |
| hrep <- rcdd::makeH(a1 = a1, b1 = b1) | |
| hrep <- rcdd::redundant(hrep)$output | |
| # vrep <- scdd(input = hrep, adjacency = TRUE, | |
| # inputadjacency = TRUE, incidence = TRUE, inputincidence = TRUE) | |
| # matrix_of_vrep <- vrep$output | |
| c(sprintf("Test.%s", opt), "H-representation", "begin", sprintf("%.0f %.0f rational", nrow(hrep), ncol(hrep)), | |
| apply(hrep, 1, paste, collapse = " "), "end") | |
| } | |
| read_problem <- function(file, opt, obj) { | |
| # The Primal LP | |
| c0 <- obj$c0 | |
| n <- nrow(c0) | |
| A_e <- obj$R | |
| p <- obj$parameters | |
| m_e <- nrow(A_e) | |
| A_l <- obj$iqR | |
| if (is.null(A_l)) A_l <- matrix(data = 0, nrow = 0, ncol = n) | |
| b_l <- obj$iqb | |
| if (is.null(b_l)) b_l <- matrix(data = 0, nrow = 0, ncol = 1) | |
| if (!is.numeric(A_l) || !is.numeric(b_l)) stop("Inequality entries must be numeric.") | |
| if (ncol(A_l) != n) stop("Dimension mismatch of inequality constaints.") | |
| m_l <- nrow(A_l) | |
| if (nrow(b_l) != m_l) stop("Dimension mismatch of inequality constaints.") | |
| m <- m_l + m_e | |
| # The Dual LP | |
| a1 <- rbind(cbind(t(A_l), t(A_e)), | |
| cbind(diag(x = 1, nrow = m_l, ncol = m_l), matrix(data = 0, nrow = m_l, ncol = m_e))) | |
| b1 <- rbind(c0, | |
| matrix(data = 0, nrow = m_l, ncol = 1)) | |
| if (opt == "max") { | |
| a1 <- -a1 | |
| b1 <- -b1 | |
| } | |
| text <- readLines(file) | |
| mrows <- text[(which(text == "begin") + 2):(which(text == "end") - 1)] | |
| matrix_of_vrep <- do.call(rbind, lapply(strsplit(mrows, " ", fixed = TRUE), \(x) as.numeric(x[x != ""]))) | |
| indices_of_vertices <- matrix_of_vrep[ , 1] == 0 & matrix_of_vrep[ , 2] == 1 | |
| vertices <- matrix_of_vrep[indices_of_vertices, -c(1, 2), drop = FALSE] # the rows of this matrix are the vertices of the convex polytope | |
| # The Bound | |
| c1_num <- rbind(b_l, 1) | |
| expressions <- apply(vertices, 1, function(y) causaloptim:::evaluate_objective(c1_num = c1_num, p = p, y = y)) | |
| elements <- paste(expressions, sep = ",", collapse = ",\n") | |
| opt_bound <- paste0(if (opt == "min") "\nMAX {\n" else "\nMIN {\n", elements, "\n}\n") | |
| opt_bound <- structure(list(expr = opt_bound, | |
| type = if (opt == "min") "lower" else "upper", | |
| dual_vertices = vertices, | |
| dual_vrep = matrix_of_vrep, | |
| expressions = expressions), | |
| class = "optbound") | |
| opt_bound | |
| } | |
| get_bounds <- function(obj, filemin, filemax) { | |
| lower_bound <- read_problem(filemin, opt = "min", obj) | |
| upper_bound <- read_problem(filemax, opt = "max", obj) | |
| bounds <- c(lower = lower_bound$expr, upper = upper_bound$expr) | |
| vreps_of_duals <- list(lower = lower_bound$dual_vrep, upper = upper_bound$dual_vrep) | |
| structure(list(bounds = bounds, | |
| expressions = list(lower = lower_bound$expressions, | |
| upper = upper_bound$expressions), | |
| logs = NULL, | |
| bounds_function = interpret_bounds(bounds, obj$parameters)), class = "balkebound") | |
| } | |
| b <- initialize_graph(graph_from_literal(Z-+ X, X -+ Y, Ur -+ X, Ur -+ Y)) | |
| obj <- analyze_graph(b, constraints = NULL, effectt = "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}") | |
| coptim_bounds <- optimize_effect_2(obj) | |
| writeLines(write_problem("min", obj), "ivbpmin.ine") | |
| writeLines(write_problem("max", obj), "ivbpmax.ine") | |
| ## | |
| # mpirun -np 4 mplrs ivbpmin.ine ivbpmin.ext | |
| ## mpirun -np 4 mplrs ivbpmax.ine ivbpmax.ext | |
| lrsbounds <- get_bounds(obj, "ivbpmin.ext", "ivbpmax.ext") | |
| testdat <- as.list(sample_distribution(obj$causal_model)) | |
| do.call(lrsbounds$bounds_function, testdat) | |
| do.call(coptim_bounds$bounds_function, testdat) | |
| ## johannes' example | |
| # # define reduced graph | |
| graph_reduced <- initialize_graph(igraph::graph_from_literal(Z -+ A, A-+Y, G -+ A, | |
| Ur -+ G, Ur -+ F, Ur-+A, Ur -+ Y)) | |
| plot(graph_reduced) | |
| V(graph_reduced)$name | |
| #define parameters | |
| N_Z <- 2 | |
| N_A <- 3 | |
| N_X <- 2 | |
| N_G <- 3 | |
| N_F <- 3 | |
| N_Y <- 2 | |
| # # Assuming V(graph)$name = "Z" "A" "Y" "G" "Ur" "F" | |
| V(graph_reduced)$leftside <- c(1, 0, 0, 0, 0, 0) | |
| V(graph_reduced)$latent <- c(0, 0, 0, 0, 1, 0) | |
| V(graph_reduced)$nvals <- c(N_Z, N_A, N_Y, N_G, 2, N_F) | |
| objective_F <- paste0("p{Y(A = ", 0: (N_A - 1) , ") = 1; F = ", 0: (N_G - 1),"}", collapse = " + ") #nolint | |
| objective_G <- paste0("p{Y(A = ", 0: (N_G - 1 ), ") = 1; G = ", 0: (N_G - 1),"}", collapse = " - ") #nolint | |
| bound_reduced <- paste0(objective_F, " - ", objective_G) | |
| obj_reduced <- analyze_graph(graph_reduced, constraints = NULL, effectt = bound_reduced) | |
| # bounds_reduced <- optimize_effect_2(obj_reduced) | |
| writeLines(write_problem("min", obj_reduced), "johanmin.ine") | |
| writeLines(write_problem("max", obj_reduced), "johanmax.ine") | |
| ## mpirun -np 12 mplrs johanmin.ine johanmin.ext | |
| ## mpirun -np 12 mplrs johanmax.ine johanmax.ext | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment