Computational Times

# Step 0: Compute size per unique instance and assign quantile bins as factors
instance_sizes <- all_results %>%
  distinct(instance, n_vars, n_ctrs) %>%
  mutate(
    size = n_vars + n_ctrs,
    size_group = ntile(size, 5),
    size_group = factor(size_group)  # convert to factor here
  )

# Step 1: Create readable interval labels for each bin, keeping size_group as factor
group_labels <- instance_sizes %>%
  group_by(size_group) %>%
  summarise(
    min_size = min(size),
    max_size = max(size),
    .groups = "drop"
  ) %>%
  mutate(
    label = paste0(min_size, "--", max_size),
    size_group = factor(size_group)  # factor again for safety
  )

# Step 2: Join size and group info into main data (filter methods first)
all_results_with_size <- all_results %>%
  filter(method %in% c("MibS", "MILP", "MILP-first", "PADM-1000")) %>%
  left_join(instance_sizes %>% select(instance, size_group, size), by = "instance") %>%
  left_join(group_labels, by = "size_group") %>%
  mutate(size_group_label = factor(label, levels = group_labels$label))  # preserve order

# Step 3: Count total instances per group (unique problems)
total_counts <- all_results_with_size %>%
  distinct(instance, size_group_label) %>%
  count(size_group_label, name = "count") %>%
  mutate(type = "Total")

# Step 4: Count solved instances per method
solved_counts <- all_results_with_size %>%
  filter(time < 3600) %>%
  group_by(size_group_label, method) %>%
  summarise(count = n(), .groups = "drop") %>%
  rename(type = method)

# Step 4.1: Fill in missing combinations (show 0 count bars)
all_methods <- unique(all_results_with_size$method)
all_groups <- unique(all_results_with_size$size_group_label)

solved_counts <- solved_counts %>%
  complete(size_group_label = all_groups, type = all_methods, fill = list(count = 0))

# Step 5: Combine and plot
plot_data <- bind_rows(total_counts, solved_counts)

ggplot(plot_data, aes(x = size_group_label, y = count, fill = type)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Instance and Solved Counts by Size Interval",
    x = "Size Interval",
    y = "Count",
    fill = "Type"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

wide_data = plot_data %>%
  pivot_wider(
    names_from = type,
    values_from = count,
    values_fill = 0  # Fill missing combinations with 0
  )

#write.csv(wide_data, "barplot_solved.csv", row.names = FALSE, quote = FALSE)

plot_ecdf_of_time = function(data, log = FALSE) {
  p = ggplot(data, aes(x = time, color = method)) +
    stat_ecdf(geom = "step", linewidth = 1) +
    labs(
      title = "ECDF of Time",
      x = "Time (s)",
      y = "Instances (%)"
    ) +
    theme_minimal() +
    scale_color_brewer(palette = "Set1")  
  
  if (log) {
    p = p + scale_x_log10()
  }
  
  return(p)
}

p = plot_ecdf_of_time(all_results %>% filter(!startsWith(method, "PADM-interdiction")))
p

export_data_from_ggplot_ecdf_csv(p, "ecdf_time_all.csv")

padm_raw = get_results_from_hpc("PADM") %>% filter(method == "PADM-100")
ggplot(padm_raw, aes(x = lambda, color = method)) +
    stat_ecdf(geom = "step", linewidth = 1) +
    scale_x_log10() +
    labs(
      title = "ECDF of Lambda bound",
      x = "Value",
      y = "Instances (%)"
    ) +
    theme_minimal() +
    scale_color_brewer(palette = "Set1")  

summary(padm_raw$lambda)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.000e+00 2.582e+03 1.727e+04 4.650e+18 1.196e+05 1.000e+20

# padm_raw %>% filter(lambda > 1e10)

Optimality Gap

plot_ecdf_of_gap_wrt_best_known = function(data) {
  results_with_gaps = data %>%
  filter(!is.na(best_known_obj)) %>% 
  filter(bobilib_status == "optimal") %>%
  filter(time < 3600) %>%
  mutate(gap = abs(best_obj - best_known_obj) / (1e-10  + abs(best_known_obj))) %>%
  mutate(gap = ifelse(gap < 1e-4, 1e-4, gap))
  
  #print(results_with_gaps %>% filter(method == "MILP-first" & gap > .8 & gap < 1.2 ))

  p = ggplot(results_with_gaps, aes(x = gap, color = method)) +
    stat_ecdf(geom = "step", linewidth = 1) +
    #scale_x_log10() +
    coord_cartesian(xlim = c(1e-4, 2e1)) +
    labs(
      title = "ECDF of Gap",
      x = "Gap (%)",
      y = "Instances (%)"
    ) +
    theme_minimal() +
    scale_color_brewer(palette = "Set1")
  
  return(p)
}

plot_ecdf_of_gap_wrt_best_known(all_results %>% filter(method == "PADM-100" | method == "PADM-1000"))

solved_by_both = solved_by_both = all_results %>% filter(method == "PADM-100" | method == "PADM-1000" | method == "MILP-first") %>%
  group_by(instance) %>%
  filter(all(time < 3600 & time >= 0)) %>%
  ungroup()

p = plot_ecdf_of_gap_wrt_best_known(solved_by_both)
p

export_data_from_ggplot_ecdf_csv(p, "ecdf_gap_all.csv")

Status change in BOBILib

all_results_with_status = all_results %>%
  filter(method == "PADM-1000" & time < 3600) %>%
  mutate(
    new_bobilib_status = ifelse(bobilib_status == "open with feasible point" & best_obj < best_known_obj - 1e-6, "was open with feasible point\n and is improving", bobilib_status),
    new_bobilib_status = ifelse(bobilib_status == "optimal" & abs(best_obj - best_known_obj) < 1e-6, "was optimal\n and found optimal point", new_bobilib_status),
    new_bobilib_status = ifelse(bobilib_status == "open", "was open", new_bobilib_status),
    new_bobilib_status = ifelse(bobilib_status == "optimal" & abs(best_obj - best_known_obj) > 1e-6, "was optimal\n and found feasible point", new_bobilib_status),
    new_bobilib_status = ifelse(new_bobilib_status == "open with feasible point" & best_obj >= best_known_obj, "was open with feasible point\n and found non-improving\n feasible point", new_bobilib_status)
  )

padm_found_feasible_point = all_results %>%
  filter(method == "PADM-1000" & time < 3600)

was_open_and_is_improving = length(padm_found_feasible_point %>% filter(bobilib_status == "open with feasible point" & best_obj < best_known_obj - 1e-6))
was_optimal_and_found_optimal = length(padm_found_feasible_point %>% filter(bobilib_status == "optimal" & abs(best_obj - best_known_obj) < 1e-6))
was_open_and_found_feasible = length(padm_found_feasible_point %>% filter( bobilib_status == "open" ))
was_optimal_and_found_feasible = length(padm_found_feasible_point %>% filter(bobilib_status == "optimal" & abs(best_obj - best_known_obj) > 1e-6))

ggplot(all_results_with_status, aes(x = new_bobilib_status)) +
  geom_bar(fill = "steelblue", width = .9) +
  geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.5) +
  labs(title = "Status change in BOBILib",
       x = "Status",
       y = "Count") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 0, hjust = .5, vjust = 1, lineheight = 1)
  )

all_results %>% 
  filter(method == "PADM-1000" | method == "MILP-first") %>%
  filter(time < 3600) %>%
  filter(  (bobilib_status == "open") ) %>%
  select(instance) %>%
  unique() %>% 
  count()

##    n
## 1 66

Computational time

all_results_common_to_padm_and_padm_interdiction = 
solved_by_both = all_results %>% filter(method == "PADM-1000" | method == "PADM-interdiction-1000") %>%
  group_by(instance) %>%
  filter(all(time >= 0)) %>%
  ungroup()

solved_by_both = solved_by_both = all_results %>% filter(method == "PADM-1000" | method == "PADM-interdiction-1000") %>%
  group_by(instance) %>%
  filter(all(time < 3600 & time >= 0)) %>%
  ungroup()

p = plot_ecdf_of_time(all_results_common_to_padm_and_padm_interdiction)
p

export_data_from_ggplot_ecdf_csv(p, "ecdf_time_interdiction.csv")

Optimality gap

plot_ecdf_of_gap_wrt_best_known(all_results %>% filter(method == "PADM-interdiction-1000"))

p = plot_ecdf_of_gap_wrt_best_known(solved_by_both)
p

export_data_from_ggplot_ecdf_csv(p, "ecdf_gap_interdiction.csv")

#to_check = read.csv("to_check.csv", col.names = c("tag", "instance", "check", "reason"))
#to_check = to_check %>% mutate(instance = basename(instance), instance = gsub("\\.aux$", "", instance))

#all_results %>% filter(method == "MILP-first" & instance %in% to_check$instance & time < 3600)