Instance Statistics

instances = read.csv("instances.csv", header = FALSE)
colnames(instances) = c("tag", "instance", "mps", "n_vars", "n_ctrs")

# Remove tag
instances$tag = NULL

# Clean instance name
instances$instance = sub("_.*", "", basename(instances$instance))
instances$mps = basename(instances$mps)

instances = instances[!duplicated(instances$instance), ]

paged_table(instances)

Results

time_limit = 1800

data = read.csv("results.csv", header = FALSE)
stats = c("time", "status", "reason", "objective", "n_outer_iterations", "n_inner_iterations","l1_norm", "l2_norm")
colnames(data) = c("tag", "instance", "n_vars", "n_ctrs", "desired_space_dim", "n_mutable_coefficients", "n_mutable_costs", "n_mutable_rhs", "method", "update_rule", "update_rule_parameter", "norm", "initial_penalty",  "warm_start", "sol_file", paste0("warm_start_", stats), stats, "solution_ok", "unconstrained_obj", "constrained_obj", "gap")

# Remove tag
data$tag = data$n_vars = data$n_ctrs = NULL

# Add a one-word solver description
data$solver = paste0(data$method, " - ", data$update_rule, " ", data$update_rule_parameter, " - ", data$norm, " - init ", data$initial_penalty, " - warm start ", data$warm_start)

# Add 'n_mutable_columns' based on instance name
data$n_mutable_columns = gsub(".*_(\\d+)$", "\\1", data$instance)

# Clean instance name
data$full_instance = basename(data$instance)

data$instance = sub("_.*", "", data$full_instance)

data$total_time = ifelse(is.na(data$warm_start_time), 0, data$warm_start_time) + data$time

data$solved = data$status == "Feasible" & data$total_time < time_limit

n_unsolved = sum(!data$solved)
if (n_unsolved > 0) {
  data[!data$solved,]$total_time = time_limit
}

data$n_mutable_columns <- as.numeric(as.character(data$n_mutable_columns))  # Convert to numeric

Merge Instances and Results

data = merge(data, instances[, c("instance", "n_vars", "n_ctrs")], 
                     by = "instance",
                     all.x = TRUE)

Summary à la Kurtz

This is the same table as Table 6 in “Counterfactual Explanations for Linear Optimization”, J. Kurtz (2024).

We focus on l1-norm with warm-start.

sub_data = data %>%
  filter(method == "PADM",
         update_rule == "adapt",
         norm == "l1",
         initial_penalty == 5e2,
         warm_start == 1)

var_bounds <- list(c(0, 534), c(534, 2167), c(2167, 22275))
ctr_bounds <- list(c(0, 351), c(351, 906), c(906, 16675))
labels <- c("small", "medium", "large")

summary_table = NULL

for (i in seq_along(labels)) {
  for (j in seq_along(labels)) {
    
    sub_summary_table <- sub_data %>%
      filter(
        n_vars >= var_bounds[[i]][1] & n_vars <= var_bounds[[i]][2],
        n_ctrs >= ctr_bounds[[j]][1] & n_ctrs <= ctr_bounds[[j]][2]
      ) %>%
      mutate(
        var_cat = labels[i],
        ctr_cat = labels[j]
      ) %>%
      group_by(var_cat, ctr_cat, n_mutable_columns) %>%
      summarise(
        `# inst.` = n_distinct(instance),  # Count number of instances
        `feasible (in %)` = mean(solved) * 100,  # Percentage of feasible instances
        `# mutable objective param.` = mean(n_mutable_costs, na.rm = TRUE),  # Avg mutable objective parameters
        `# mutable constraint param.` = mean(n_mutable_coefficients, na.rm = TRUE),  # Avg mutable constraint parameters
        .groups = "drop"
      ) %>%
      arrange(var_cat, ctr_cat, n_mutable_columns) %>%  # Sort by var_cat, ctr_cat, and n_mutable_columns
      ungroup()
    
    summary_table = rbind(summary_table, sub_summary_table)    
    
  }
}

summary_table %>%
  kable("html", align = "c", col.names = c()) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE, position = "center") %>%
  add_header_above(c("n" = 1, "m" = 1, "# mut. columns" = 1, "# inst." = 1, "feasible (in %)" = 1, 
                     "# mutable objective param." = 1, "# mutable constraint param." = 1)) %>%
  group_rows("small", 1, 6) %>%
  group_rows("medium", 7, 15) %>%
  group_rows("large", 16, 24)

Solved Instances by Number of Mutable Columns

bar_data = sub_data %>%
  filter(solved == TRUE) %>%  # Only consider solved instances
  group_by(instance, n_mutable_columns) %>%
  summarise(solved_count = n(), .groups = "drop")  # Count solved instances

# Create bar plot
ggplot(bar_data, aes(x = instance, y = solved_count, fill = factor(n_mutable_columns))) +
  geom_bar(stat = "identity", position = "dodge") +  # Bar plot with bars side-by-side
  coord_flip() +
  labs(
    title = "Number of Solved Instances by n_mutable_columns for Each Instance",
    x = "Instance",
    y = "Number of Solved Instances",
    fill = "n_mutable_columns"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position = "top")  # Rotate x-axis labels for better readability

Computational Times

for (norm in unique(data$norm)) {
  
  sub_data = data[data$norm == norm, ]
  
  plot = ggplot(sub_data, aes(x = total_time, group = solver, color = solver)) + 
    stat_ecdf(geom = "step") + 
    labs(title = paste0("ECDF of Time for Each Solver using ", norm, " norm"),
         x = "Time",
         y = "ECDF",
         color = "Solver") +
    theme_minimal() +
    theme(legend.position = "bottom") +
    scale_x_continuous(breaks = seq(0, max(sub_data$total_time), by = 60), limits = c(0, time_limit)) +
    scale_y_continuous(breaks = seq(0, 1, by = 0.1))

  print(plot)
}

for (norm in unique(data$norm)) {
  
  sub_data = data[data$norm == norm,]
  
  sub_data = sub_data %>%
    group_by( sub("^(.+_[0-9]+)_.*", "\\1", full_instance)) %>%
    filter(any(solved == TRUE)) %>%
    ungroup()

  plot = ggplot(sub_data, aes(x = total_time, group = solver, color = solver)) + 
    stat_ecdf(geom = "step") + 
    labs(title = paste0("ECDF of Time for Each Solver using ", norm, " norm"),
         x = "Time",
         y = "ECDF",
         color = "Solver") +
    theme_minimal() +
    theme(legend.position = "bottom") +
    scale_x_continuous(breaks = seq(0, max(sub_data$total_time), by = 300), limits = c(0, time_limit)) +
    scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +
    facet_wrap(~ n_mutable_columns, ncol = 3)  # Facet by mutable_columns, 3 columns
  
  ########## SAVING DATA TO FILE FOR TIKZ ##########
  ecdf_data = data.frame(time = sub_data$total_time)
  
  for (n_mutable_columns_val in unique(sub_data$n_mutable_columns)) {
    
    for (warm_start_val in unique(sub_data$warm_start)) {
      
      facet_data = sub_data %>% filter(n_mutable_columns == n_mutable_columns_val & warm_start == warm_start_val)
      
      ecdf_values = ecdf(facet_data$total_time)(sub_data$total_time)
      
      ecdf_data = cbind(ecdf_data, y = ecdf_values)
      colnames(ecdf_data)[ncol(ecdf_data)] = paste0(n_mutable_columns_val, "_w", warm_start_val)
      
    }
  }
  
  ecdf_data <- as.data.frame(ecdf_data)
  ecdf_data <- ecdf_data[order(ecdf_data$time), ]
  ecdf_data = ecdf_data %>% filter(time < 1800)
  
  indices <- seq(1, length(ecdf_data$time), length.out = 100)
  ecdf_data <- ecdf_data[indices, ]
  
  file_name <- paste0("ecdf_", norm, ".csv")
  write.csv(ecdf_data, file_name, row.names = FALSE)
  ###################################################
  
  print(plot)
   
}

Scatter Plot of Solved and Unsolved Instances

sub_data = data[data$update_rule == "adapt" & data$norm == "l1" & data$initial_penalty == 5e2,]

# Create a scatter plot with colors based on the normalized time
ggplot(sub_data, aes(x = n_vars, y = n_ctrs, color = as.numeric(solved))) +
  geom_point(alpha = .5, size = 1) +  # Transparency and point size
  scale_x_log10() +  # Logarithmic scale for n_vars
  scale_y_log10() +  # Logarithmic scale for n_ctrs
  scale_color_gradient(low = "red", high = "green") +  # Color gradient
  labs(title = "Scatter Plot of (n_vars, n_ctrs) with Log Scale and Transparency",
       x = "Number of Variables (n_vars)",
       y = "Number of Constraints (n_ctrs)",
       color = "Solved") +
  theme_minimal()  # Minimal theme for clean look