# Using column generation in constraint-and-column generation for adjustable robust optimization > FLP

```{r echo=FALSE}
library(rmarkdown)
library(kableExtra)
library(tidyr)
library(ggplot2)
suppressPackageStartupMessages(library(dplyr))
```

## Loading the data

Our results can be found in the `results.flp.csv` file with the following columns:

- "tag": a tag always equal to "result" used grep the result line in our execution log file.
- "instance": the path to the instance.
- "standard_phase_time_limit": the time limit for the standard phase (i.e., without using CG).
- "master_solver": the solver used for solving the CCG master problem: STD for standard, i.e., Gurobi, CG for column generation.
- "status": the final status.
- "reason": the final status reason.
- "has_large_scaled": true if the CG phase has been started, false otherwise.
- "n_iterations": the number of iterations.
- "total_time": the total time to solve the problem.
- "master_time": the time spent solving the master problem.
- "adversarial_time": the time spent solving the adversarial problem.
- "best_bound": the best bound found.
- "best_obj": the best feasible point value.
- "relative_gap": the final relative gap.
- "absolute_gap": the final absolute gap.
- "adversarial_unexpected_status": the status of the adversarial problem solver if it is not Optimal.
-  "with_heuristic": true if the CG-based heuristic is used.
- "n_facilities": the number of facilities in the instance.
- "n_customers":  the number of customers in the instance.
- "Gamma": the value for the uncertainty budget $\Gamma$.
- "blank": this column is left blank.

```{r}
data = rbind(
  read.csv("results.csv", header = FALSE) 
)
colnames(data) = c("slurm_file", "tag", "instance", "standard_phase_time_limit_raw", "Gamma", "with_heuristic", "standard_phase_time_limit", "master_solver", "status", "reason", "has_large_scaled", "n_iterations", "total_time", "master_time", "adversarial_time", "best_bound", "best_obj", "relative_gap", "absolute_gap", "second_stage.mean", "second_stage.std_dev", "adversarial_unexpected_status", "memory_used", "memory_limit")
data = data %>% mutate(instance = basename(instance),
                n_facilities = as.numeric(sub(".*F(\\d+)_C.*", "\\1", instance)),
                n_customers = as.numeric(sub(".*C(\\d+)_.*", "\\1", instance))
                )
```

We start by removing the "tag" and the "blank" columns.

```{r}
data = data[, !(names(data) %in% c("tag", "blank"))]
```

We then remove instances which were not considered in our experimental results because
none of the approaches could solve sufficiently many instances.

```{r}
#data = data[data$n_customers != 40,]
#data = data[data$n_facilities == 10,]
```

For homogeneity, we fix the total_time of unsolved instances to the time limit. 

```{r}
if (sum(data$tag == "iteration") > 0 ) {
  data[data$tag == "iteration",]$total_time = 10800 
}

if (sum(data$total_time > 10800) > 0 ) {
  data[data$total_time > 10800,]$total_time = 10800 
}
```

Then, we create a column named "method" which gives a specific name to each method, comprising the approach for solving the CCG master problem,
the time limit of the standard phase and a flag indicating if the CG-based heuristic was used.

```{r}
data$method = paste0(data$master_solver, "_TL", data$standard_phase_time_limit, "_H", data$with_heuristic)
unique(data$method)
data = data[data$method != "STD_TLInf_H1" & data$method != "CG_TL120_H0" & data$method != "CG_TL120_H1",]
```

Our final data reads.

```{r, echo = FALSE}
paged_table(data)
```

## Sanity Check

```{r}
# Define the relative tolerance
tolerance <- 10^-2

# Filter the data where time < 10800 and group by 'instance'
validation <- data %>%
  filter(total_time < 10800) %>%
  group_by(paste0(instance, "_", Gamma)) %>%
  summarise(min_best_obj = min(best_obj), max_best_obj = max(best_obj)) %>%
  mutate(valid = abs(max_best_obj - min_best_obj) / max(abs(min_best_obj), 1) <= tolerance)

# Check if all instances are valid
if (all(validation$valid)) {
  print("All methods find the same best_obj value within the relative tolerance for all instances.")
} else {
  print("Methods do not find the same best_obj value within the relative tolerance for some instances.")
  print(validation %>% filter(!valid))  # Show the instances that failed validation
}
```

## Empirical Cumulative Distribution Function (ECDF)

We plot the ECDF of computation time over our set of instances for all approaches.

```{r}
ggplot(data, aes(x = total_time, col = method)) + stat_ecdf(pad = FALSE) +
  coord_cartesian(xlim = c(0,10800)) +
  # scale_x_log10() +
  theme_minimal()


ggplot(data, aes(x = memory_used, col = method)) +
  stat_ecdf(pad = FALSE) +
  theme_minimal()
```

```{r}
ggplot(data, aes(x = total_time / n_iterations, col = method)) + stat_ecdf(pad = FALSE) +
  # scale_x_log10() +
  theme_minimal()
```
We export these results in csv to print them in tikz.

```{r}
data_with_ecdf = data %>%
  group_by(method) %>%
  arrange(total_time) %>%
  mutate(ecdf_value = ecdf(total_time)(total_time)) %>%
  ungroup()

for (method in unique(data_with_ecdf$method)) {
  output = data_with_ecdf[data_with_ecdf$method == method,]
  output = output[,c("total_time", "ecdf_value")]
  output$log_total_time = log10(output$total_time)
  output = output[output$total_time < 10800,]
  write.csv(output, file = paste0(method, ".csv"), row.names = FALSE)
}
```

## Summary table

In this section, we create a table summarizing the main outcome of our computational experiments.

We first focus on the solved instances.

```{r}
summary_data_lt_10800 <- data %>%
  filter(total_time < 10800) %>%
  group_by(n_facilities, n_customers, Gamma, method) %>%
  summarize(
    avg_total_time = mean(total_time, na.rm = TRUE),
    avg_master_time = mean(master_time, na.rm = TRUE),
    avg_adversarial_time = mean(adversarial_time, na.rm = TRUE),
    avg_n_iterations = mean(n_iterations, na.rm = TRUE),
    sum_has_large_scaled = sum(has_large_scaled),
    num_lines = n(),
    .groups = "drop"
  ) %>%
  ungroup() %>%
  arrange(n_facilities, n_customers, Gamma, method)
```

Then, we compute averages over the unsolved instances.

```{r}
summary_data_ge_10800 <- data %>%
  filter(total_time >= 10800) %>%
  group_by(n_facilities, n_customers, Gamma, method) %>%
  summarize(
    avg_n_iterations_unsolved = mean(n_iterations, na.rm = TRUE),
    num_lines_unsolved = n(),
    .groups = "drop"
  ) %>%
  ungroup() %>%
  arrange(n_facilities, n_customers, Gamma, method)
```

Finally, we merge our results.

```{r}
transposed_data_lt_10800 <- summary_data_lt_10800 %>%
  pivot_wider(names_from = method, values_from = avg_total_time:num_lines)

transposed_data_ge_10800 <- summary_data_ge_10800 %>%
  pivot_wider(names_from = method, values_from = avg_n_iterations_unsolved:num_lines_unsolved) %>%
  select(-n_facilities, -n_customers, -Gamma) 

transposed_data_lt_10800 %>% kable()
transposed_data_ge_10800 %>%  kable()

#cbind(
#  transposed_data_lt_10800,
#  transposed_data_ge_10800
#) %>%
#  kable() %>%
#  kable_styling(full_width = FALSE, position = "center")

```


## Second-stage Deviations

```{r}
ggplot(data, aes(x = method, y = second_stage.std_dev / abs(second_stage.mean))) +
  geom_boxplot() +
  labs(title = "Boxplot of std.dev",
       x = "Method",
       y = "Number of Iterations") +
  theme_minimal()
```

