README_old.Rmd

---
output: github_document
editor_options: 
  markdown: 
    wrap: sentence
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  dpi = 300
)
```

# Zoning for Autonomous Vehicles

<!-- badges: start -->

<!-- badges: end -->

The goal of zav (Zoning for Autonomous Vehicles) is to show the implementation of the problem instance generation and simulation experiment conducted for the project.

## Installation

You can install the development version of zav from [GitHub](https://github.com/) with:

``` r
# install.packages("devtools")
# note that there is a GitHub only dependency
# devtools::install_github("cuhklinlab/SWKM")
devtools::install_github("Rosenkrands/zav")
```

## Instance generation

To generate a problem instance we can utilize the `generate_2d_instance` function.

```{r example}
library(zav)
library(ggplot2)
library(dplyr)

instance <- generate_2d_instance(no_of_points = 100)
# plot_point(instance = instance)
# ggplot(instance$data) +
#     geom_point(aes(x,y,size=`Arrival rate`),
#                         shape = 21, fill = alpha("black", .2)) +
#     theme_void()

ggplot(instance$data) +
    geom_point(aes(x,y),
                        shape = 21, fill = alpha("black", .2)) +
    theme_void() + coord_fixed()

ggsave("instance_overview_uniform.pdf", width = 6, height = 3)
```

## Zoning solutions

### Weighted K-Means

The below code chunk shows how we can utilize the `solve_wkmeans` function to generate a solution for our problem instance.

```{r solution}
solution_wkm <- solve_wkmeans(instance, no_of_centers = 5, type = "swkm")
plot_bases(solution = solution_wkm)
```

### Genetic Algorithm

The below chunk shows how we can utilize the `solve_ga` function to generate a solution for our problem instance.

First we need to precalculate centroids to use as input for the GA.

```{r, cache=T}
centroids <- grid_centroids(instance, dimension = 5)
```

Having now the centroids and distances between demand points and centroids, we are able to give this as input for the GA.
As a default the GA will have a maximum of 10 iterations for demonstration purposes, in reality we would have a much higher number of iterations.

```{r, cache=T}
solution <- solve_ga(instance, centroids, no_of_centers = 5, obj = "SAFE")
```

```{r}
plot_bases(solution)
```

## Simulation

Here it is shown how you can make a simulation based on the solution just found:

```{r simulation, cache = T}
#' The below line is for testing purposes
#' solution = solution_wkm; seed = 1;n_replications = 1;flight = "free";max_dist = 1000000;LOS = 600;warmup = 0;speed_agent = .25;verbose = F; dist_haversine = T
simulation_result <- simulation(
  solution = solution_wkm,
  seed = 1,
  n_replications = 1,
  flight = "zoned",
  queue = T,
  max_dist = 1000000,
  LOS = 600,
  warmup = 0,
  speed_agent = .25,
  verbose = F
)
```

```{r}
# simulation_result$metrics[[1]]$response_time_performance
```

```{r distances, fig.dim=c(8,3)*1.1}
simulation_result$metrics[[1]]$distances %>%
  ggplot(aes(x = time, 
             y = distance,
             color = paste0(id1,'-',id2))) +
  geom_line() +
  labs(color = "Pair") +
  theme(legend.position = "none")
```

### Free-flight/Soft zoning

As an argument to the simulation function we have `flight` that together with the `max_dist` argument allow the user to obtain simulation results from a "Free-flight" setup.
Here we will give a brief introduction to the concept.

If we set `fligth = "free"` and `max_dist` arbitrarily large any UAV would be able to service any demand.
This is what we would refer to as free-flight.

However when any UAV can service any demand then mean response will suffer from long travel distances.
Therefore we propose using the following scheme to decide on a maximum distance for the service area

$\\text{r} = \\text{dist}\_\\text{max} + \\alpha \\cdot \\text{number of UAVs}$

where $\\alpha$ denotes the scaling factor and $\\text{dist}\_\\text{max}$ denotes the maximum distance between a demand point and its base location.

To further illustrate let us consider the following situation.

```{r}
instance <- generate_2d_instance(no_of_points = 40)
solution <- solve_wkmeans(instance, no_of_centers = 4, type = "swkm")

plot_bases(solution) + coord_fixed()
```

We can the illustrate the service area for different values of $\\alpha$ as follows.

```{r, echo = F, message=F}
centroids <- solution$instance %>%
    select(`Centroid id`, x.centroid, y.centroid) %>%
    distinct() %>%
    mutate(r = (solution$instance %>%
      mutate(distance = sqrt((x - x.centroid)^2 + (y - y.centroid)^2)) %>%
      summarise(distance = max(distance)))$distance
    )

ggplot(solution$instance) +
  ggforce::geom_circle(
    data = centroids %>% 
      left_join(tibble(n = c(0)), by = character()) %>%
      mutate(r = r + n*4,
             `Scaling factor` = factor(n)),
    aes(x0 = x.centroid, y0 = y.centroid, r = r, fill = `Centroid id`, linetype = `Scaling factor`),
    alpha = .1, color = "grey"
  ) +
  geom_segment(aes(x = x, y = y, xend = x.centroid, yend = y.centroid),
               color = "gray") +
  geom_point(aes(x,y, fill = `Centroid id`, size = `Arrival rate`, alpha = .2),
                      shape = 21, color = "black") +
  geom_point(
    data = centroids,
    aes(x.centroid, y.centroid), shape = 15, size = 2
  ) +
  theme_void() + coord_fixed() +
  scale_alpha(guide = "none") +
  scale_fill_discrete(guide = "none") +
  scale_linetype(guide = "none") +
  scale_size(guide = "none")

# ggsave("free_flight_distance_constraints.pdf", width = 6,height = 5)
```

## Experiment Results

To get analyze the results from the experiment conducted for the project we do the following.

```{r}
results <- experiment_results()

# Renaming the labels
results$simulation <- results$simulation %>%
  mutate(`Flight configuration` =
           factor(as.character(`Flight configuration`), 
                  levels = c("Zoned", "Free-flight: 0", "Free-flight: .2", "Free-flight: .5", "Free-flight: no constraint"),
                  labels = c("Zoning", "Maximum flight range", "Free-flight: .2", "Free-flight: .5", "Free-flight")))
```

The `results` variable is a list of 3 `tibbles` containing information about instances, solutions and simulations respectively.

### Instances

If we take a look at the instances first we can see that there is a total of `r nrow(results$instance)` instances.
These are distributed across the two arrival rate variances as seen below.

```{r, message = F}
results$instance |> 
  group_by(`Arrival rate distribution`, `Arrival rate variance`) |> 
  summarise(n = n())
```

### Solutions

Taking a look at the solutions we see that there are a total of `r nrow(results$solution)` solutions.
These are distributed across:

-   `Solution method` that have `r nrow(unique(results$solution[, "Solution method"]))` levels: `r unique(results$solution[, "Solution method"][[1]])`. <!-- * `Number of UAVs` that have `r nrow(unique(results$solution[, "Number of UAVs"]))` levels: `r unique(results$solution[, "Number of UAVs"][[1]])`. -->

#### Solution comparsion

If we start by comparing ga-tot and wkm-swkm we see that the latter is able to outperform the first.

```{r, echo=F, message=F, fig.dim=c(8,3)*1.1}
results$solution |> 
  filter(`Solution method` != "ga-safe") |> 
  group_by(`Solution method`, `Arrival rate variance`, `Arrival rate distribution`, `Number of UAVs`) |> 
  summarise(sd = sd(TOT), TOT = mean(TOT)) |> 
  ggplot(aes(x = `Arrival rate variance`, y = TOT, group = `Solution method`, color = `Solution method`)) +
  geom_point() +
  geom_line() +
  geom_errorbar(aes(ymin = TOT - sd, ymax = TOT + sd), width = .05)
```

### Simulations

Looking at the simulation results we can see how the utilization is affected by the queuing strategy.
First we can take a look at the zoned solution approach.

```{r, echo=F, message=F, utilization, fig.dim=c(8,3)*1.1}
results$simulation %>%
  select(`Instance id`, `Arrival rate distribution`, `Arrival rate variance`, `Solution method`, `Number of UAVs`, `Flight configuration`, `Queue strategy`, utilization) %>%
  filter(`Solution method` == "wkm-swkm", `Flight configuration` == "Zoned",
                `Arrival rate distribution` == "uniform") %>%
  transmute(group = paste0(`Instance id`, "-", `Queue strategy`), `Instance id`, `Queue strategy`, `Arrival rate variance`, utilization) %>%
  tidyr::unnest(cols = utilization) %>%
  group_by(hours, `Queue strategy`, `Arrival rate variance`) %>%
  summarise(min_util = min(utilization), max_util = max(utilization), mean_util = mean(utilization)) %>%
  ggplot(aes(x= hours, group = `Queue strategy`)) +
  geom_line(aes(y = mean_util, color = `Queue strategy`)) +
  geom_ribbon(aes(ymin = min_util, ymax = max_util, fill = `Queue strategy`), alpha = .4) +
  scale_y_continuous(limits = c(0,1)) +
  facet_wrap(~`Arrival rate variance`, labeller = label_both)
```

If we compare with the zoned solution approach with a free-flight approach we see that the utilization go to 1 with the FCFS queue.

Looking at the simulation performance measures mean response and ploss, given the no queue strategy, it would seem that we are giving up some response time to lower the ploss for $\\alpha = 0$.
But for other values of $\\alpha$ the performance are worse of overall.

```{r queue_no, echo=F, message=F, fig.dim=c(8,3)*1.1}
results$simulation |> 
  filter(`Arrival rate distribution` == "uniform", 
                `Queue strategy` == "No queue",
                `Number of UAVs` == "high",
         `Solution method` == "wkm-swkm") |> 
  tidyr::pivot_longer(cols = c(`Mean response`, Ploss), names_to = "Measure") |> 
  mutate(Measure = factor(Measure, levels = c("Mean response", "Ploss"))) |>
  group_by(Measure,`Arrival rate variance`,`Flight configuration`, `Number of UAVs`) |> 
  summarise(error = qt(.975, df = n() - 1)*sd(value)/sqrt(n()),
            value = mean(value)) |> 
  ggplot(aes(x = `Arrival rate variance`, y = value, fill = `Flight configuration`)) +
  geom_col(position = position_dodge(.9)) +
  geom_errorbar(aes(ymin = value - error, ymax = value + error), width = .2, position = position_dodge(.9)) +
  facet_wrap(~Measure, scales = "free", labeller = label_both)

# ggsave("no_queue_metrics_flight_configuration.pdf", width = 8, height = 3)
```

Results are exaggerated for the FCFS queue strategy.

```{r queue_fcfs, echo=F, message=F, fig.dim=c(8,3)*1.1}
results$simulation |> 
  filter(`Queue strategy` == "FCFS", `Solution method` == "wkm-swkm") |> 
  ggplot(aes(x = `Demands in queue`, fill = `Flight configuration`, color = `Flight configuration`)) +
  geom_density(alpha = .2)

# ggsave("demands_in_queue_flight_configuration.pdf", width = 8, height = 2)
```

Looking at just the numbers we get for mean response the following numbers for zoned and free-flight: no constraint.

```{r, echo = F, message = F}
results$simulation |> 
  filter(`Solution method` %in% c("wkm-swkm"),
         `Flight configuration` %in% c("Zoned", "Free-flight: no constraint")) |> 
  tidyr::pivot_longer(cols = c(`Mean response`, Ploss), names_to = "Measure") |>
  filter(Measure == "Mean response") |> 
  group_by(`Solution method`, `Queue strategy`, `Flight configuration`, Measure) |>
  summarise(value = mean(value)) |> 
  tidyr::pivot_wider(names_from = `Flight configuration`, values_from = value) |> 
  dplyr::select(-`Solution method`) |> 
  knitr::kable(format = "pipe")
```

## Showing mean response

```{r}
results$simulation |> 
  filter(
    `Arrival rate distribution` == "uniform",
    `Arrival rate variance` == "low",
    `Solution method` == "wkm-swkm",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
  ) |> 
  ggplot() +
  geom_boxplot(aes(y = `Mean response`/60, color = `Flight configuration`)) +
  facet_wrap(
     ~ `Queue strategy`, 
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  labs(y = "Response time (minutes)") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

ggsave("mean_response.pdf", width = 8, height = 3)
```

with all response times not just means:

```{r}
results$simulation |> 
  filter(
    `Arrival rate distribution` == "uniform",
    `Arrival rate variance` == "low",
    `Solution method` == "wkm-swkm",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
  ) |> 
  tidyr::unnest(cols = response_times) |> 
  ggplot() +
  geom_boxplot(aes(y = response_time/60, color = `Flight configuration`)) +
  facet_wrap(
     ~ `Queue strategy`, 
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  labs(y = "Response time (minutes)") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
```


## Showing Ploss

```{r}
results$simulation |> 
  filter(
    `Arrival rate distribution` == "uniform",
    `Arrival rate variance` == "low",
    `Solution method` == "wkm-swkm",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight"),
    `Queue strategy` == "No queue"
  ) |> 
  ggplot() +
  geom_boxplot(aes(y = `Ploss`, color = `Flight configuration`)) +
  facet_grid(
     ~ `Queue strategy`, 
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  labs(y = "% of demands lost") +
  theme(legend.position = "none") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

ggsave("Ploss.pdf", width = 3, height = 2.5)
```

## Showing distance

```{r}
# results$simulation |> 
#   filter(
#     `Arrival rate distribution` == "uniform",
#     `Solution method` == "wkm-swkm",
#     `Number of UAVs` == "high",
#     `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
#   ) |> 
#   tidyr::unnest(cols = min_distances) |> 
#   group_by(`Flight configuration`, `Queue strategy`, `Arrival rate variance`) |>
#   filter(distance == min(distance)) |> 
#   ggplot() +
#   geom_boxplot(aes(y = distance, color = `Flight configuration`)) +
#   facet_grid(
#     `Queue strategy` ~ `Arrival rate variance`, 
#     labeller = label_both, 
#     scales = "free"
#   ) +
#   #scale_color_grey() +
#   theme_bw() +
#   labs(y = "Minimum distance")

results$simulation |> 
  filter(
    `Arrival rate variance` == "low",
    `Arrival rate distribution` == "uniform",
    `Solution method` == "wkm-swkm",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
  ) |> 
  ggplot() +
  geom_boxplot(aes(y = `Minimum distance`, color = `Flight configuration`)) +
  facet_grid(
     ~ `Queue strategy`,
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  labs(y = "Minimum distance") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

ggsave("minimum_distance.pdf", width = 8, height = 3)
```

```{r}
# results$simulation |> 
#   filter(
#     `Arrival rate distribution` == "uniform",
#     `Solution method` == "wkm-swkm",
#     `Number of UAVs` == "high",
#     `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
#   ) |> 
#   tidyr::unnest(cols = min_distances) |> 
#   group_by(`Flight configuration`, `Queue strategy`, `Arrival rate variance`) |>
#   filter(distance < quantile(distance, probs = .05)) |> 
#   ggplot() +
#   geom_boxplot(aes(y = distance, color = `Flight configuration`)) +
#   facet_grid(
#     `Queue strategy` ~ `Arrival rate variance`, 
#     labeller = label_both, 
#     scales = "free"
#   ) +
#   #scale_color_grey() +
#   theme_bw() +
#   labs(y = "5th percentile distance")

results$simulation |> 
  filter(
    `Arrival rate distribution` == "uniform",
    `Arrival rate variance`== "low",
    `Solution method` == "wkm-swkm",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
  ) |>
  ggplot() +
  geom_boxplot(aes(y = `5th percentile distance`, color = `Flight configuration`)) +
  facet_grid(
     ~ `Queue strategy`, 
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  labs(y = "5th percentile distance") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
  

ggsave("5th_minimum_distance.pdf", width = 8, height = 3)
```


## Example with real data

To see how the methods presented apply to real data we include the following section.
From the R package `maxcovr` we found the York crime data set.
The distribution of points from the data set looks as follows.

```{r, echo=F, message=F, fig.dim=c(8,3)*1.1}
data(york, package = "maxcovr")

ggplot2::ggplot(york %>% select(long, lat) %>% distinct()) +
  ggplot2::geom_point(ggplot2::aes(long,lat),
                      shape = 21, fill = ggplot2::alpha("black", .2)) +
  ggplot2::theme_void() + coord_fixed()

# ggsave("location_overview_v2.pdf", width = 6, height = 3)
```

We have sampled demand points from this data set, to mimic the original distribution, and made solutions and simulations.
The sample with assigned arrival rates is displayed in the following figure.

```{r}
# Read precalculated york solutions
clean_york <- function(york_solution) {
  if (grepl("GA", york_solution)) {
    solution <- readRDS(york_solution)

    solution$no_of_centers <- nrow(solution$centroids)

    solution$clusters <- solution$centroids %>%
      dplyr::mutate(`Cluster id` = dplyr::row_number())

    solution$instance <- solution$instance %>%
      dplyr::inner_join(solution$clusters %>% dplyr::select(-x,-y), by = "Centroid id") %>%
      dplyr::select(-`Centroid id`) %>%
      dplyr::mutate(`Centroid id` = as.character(`Cluster id`)) %>%
      dplyr::select(-`Cluster id`)

    return(solution)
  } else {
    return(readRDS(york_solution))
  }
}

york_solutions <- lapply(
  list.files("inst/extdata/york", full.names = T) %>% as.list(),
  clean_york
)

ggplot2::ggplot(york_solutions[[1]]$instance) +
  ggplot2::geom_point(ggplot2::aes(x,y, size = `Arrival rate`),
                      shape = 21, fill = ggplot2::alpha("black", .2)) +
  ggplot2::theme_void() + coord_fixed()

# ggsave("york_sample_with_arrival_rate.pdf", width = 8, height = 3)
```

We obtain metadata in the following way.

```{r}
# The results from the york data set are available here
york <- readRDS(system.file("extdata", "york_simulation_metadata.rds", package = "zav"))
```

When comparing the solution we can see that `wkm-swkm` achieve a lower TOT objective as compared to `ga-tot`.

```{r, echo = F, message=F, fig.dim=c(8,3)*1.1}
york |> 
  group_by(`Solution method`, `Arrival rate variance`, `Number of UAVs`) |> 
  summarise(sd = sd(TOT), TOT = mean(TOT)) |> 
  ggplot(aes(x = `Arrival rate variance`, y = TOT, group = `Solution method`, color = `Solution method`)) +
  geom_point() +
  geom_line() +
  facet_wrap(~`Number of UAVs`, labeller = label_both)
```

Comparing now utilization of UAVs in the simulations we see that with low number of UAVs, low arrival rate variance and FCFS queue strategy we have higher utilization with `ga-tot` as compared to `wkm-swkm`.

```{r, echo=F, message=F, fig.dim=c(8,6)*1.1}
york %>%
  select(`Solution method`,`Number of UAVs`, queue, `Arrival rate variance`, utilization) %>%
  tidyr::unnest(cols = utilization) %>%
  group_by(hours, `Solution method`, queue, `Arrival rate variance`, `Number of UAVs`) %>%
  summarise(min_util = min(utilization), max_util = max(utilization), mean_util = mean(utilization)) %>%
  ggplot(aes(x= hours, group = paste0(`Solution method`, queue))) +
  geom_line(aes(y = mean_util, linetype = queue, color = `Solution method`)) +
  # geom_ribbon(aes(ymin = min_util, ymax = max_util, fill = `Queue strategy`), alpha = .4) +
  scale_y_continuous(limits = c(0,1)) +
  facet_wrap(`Number of UAVs`~`Arrival rate variance`, labeller = label_both)

```

We believe this could be due to `ga-tot` having to base locations sharing the dense part of the service area, where as `wkm-swkm` only have one base location covering the dense part of the area.

```{r, echo=F, message=F, fig.dim=c(8,3.5)*1.1}
clean_york <- function(york_solution) {
  if (grepl("GA", york_solution)) {
  solution <- readRDS(york_solution)

  solution$no_of_centers <- nrow(solution$centroids)

  solution$clusters <- solution$centroids %>%
    dplyr::mutate(`Cluster id` = dplyr::row_number())

  solution$instance <- solution$instance %>%
    dplyr::inner_join(solution$clusters %>% dplyr::select(-x,-y), by = "Centroid id") %>%
    dplyr::select(-`Centroid id`) %>%
    dplyr::mutate(`Centroid id` = as.character(`Cluster id`)) %>%
    dplyr::select(-`Cluster id`)

    return(solution)
  } else {
    return(readRDS(york_solution))
  }
}

york_solutions <- lapply(
  list.files("inst/extdata/york", full.names = T) %>% as.list(),
  clean_york
)

plot_bases(york_solutions[[11]]) + theme(legend.position = "none")
ggsave("york_wkm_low_uav.pdf", width = 3, height = 3)

plot_bases(york_solutions[[9]]) + theme(legend.position = "none")
ggsave("york_wkm_high_uav.pdf", width = 3, height = 3)

plot_bases(york_solutions[[4]]) + theme(legend.position = "none")
ggsave("york_tot_low_uav.pdf", width = 3, height = 3)

plot_bases(york_solutions[[8]]) + theme(legend.position = "none")
ggsave("york_vot_low_uav.pdf", width = 3, height = 3)
```

To see how this affect the mean response for the zoned approach with no queue we can take a look at the following figure.

```{r york_simulation_metrics, echo=F, message=F, fig.dim=c(8,6)*1.1}
york |>
  filter(`Queue strategy` == "No queue",
         `Flight configuration` == "Zoned") |> 
  tidyr::pivot_longer(cols = c(`Mean response`, Ploss), names_to = "Measure") |> 
  mutate(Measure = factor(Measure, levels = c("Mean response", "Ploss"))) |>
  group_by(Measure,`Arrival rate variance`, `Solution method`, `Number of UAVs`) |> 
  summarise(value = mean(value)) |> 
  ggplot(aes(x = `Arrival rate variance`, y = value, color = `Solution method`, group = paste0(`Measure`,`Solution method`))) +
  geom_point(alpha = .75) +
  geom_line(alpha = .75) +
  facet_grid(Measure~`Number of UAVs`, scales = "free", labeller = label_both)
```

Looking at the figure we conclude that in general we see shorter response times for `wkm-swkm` but also a high Ploss when comparing to `ga-tot`.

Next we can take a look at how introducing a queue affects the findings.

```{r york_simulation_metrics_queue, echo=F, message=F, fig.dim=c(8,6)*1.1}
york |>
  filter(`Queue strategy` == "FCFS",
         `Flight configuration` == "Zoned") |> 
  tidyr::pivot_longer(cols = c(`Mean response`, `Demands in queue`), names_to = "Measure") |> 
  mutate(Measure = factor(Measure, levels = c("Mean response", "Demands in queue"))) |>
  group_by(Measure,`Arrival rate variance`, `Solution method`, `Number of UAVs`) |> 
  summarise(value = mean(value)) |> 
  ggplot(aes(x = `Arrival rate variance`, y = value, color = `Solution method`, group = paste0(`Measure`,`Solution method`))) +
  geom_point(alpha = .75) +
  geom_line(alpha = .75) +
  facet_wrap(Measure~`Number of UAVs`, scales = "free", labeller = label_both)
```

First we notice that the winner in terms of mean response time have shifted and we now see better performance for `ga-tot` as compared to `wkm-swkm`.
Looking at the demands in queue measure we also see better performance for `ga-tot` as compared to `wkm-swkm`.

### Flight configurations

Having looked at the zoned approach with different queuing strategies, this section will deal with the soft zoning strategies to see if they can improve performance.

In the next figure we show how the mean response is affected by solution method, queue strategy and arrival rate variance.
We restrict to only looking at high number of UAVs, the same pattern is true for low number of UAVs but exaggerated in case of the FCFS queue.

```{r, echo=F, message=F, fig.dim=c(8,6)*1.1}
york |> 
  filter(`Number of UAVs` == "high") |>
  ggplot(aes(x = `Arrival rate variance`, y = `Mean response`, fill = `Flight configuration`)) +
  geom_col(position = position_dodge(.9)) +
  facet_grid(`Queue strategy`~`Solution method`, scales = "free", labeller = label_both)
```

In the case of having no queue, we see the same kind of pattern as from the uniformly generated instances, here the zoned approach can achieve the best mean response as it, by design, only chooses to service the closest demands.

However this is not the case when introducing the queue, we suspect this is because now the even distribution of demand between UAVs become important as demands pile up in the queue a certain base is overworked.
This seems to be the case, in particular for the `wkm-swkm` solution method where mean response for the zoned approach are much higher than any of the free-flight configurations.

To further see the importance of balancing total arrival rate between zones we can take a look at the demands in queue metric for the high UAV case.
Interestingly we have here that the free-flight: no constraint in all situation are able to serve all demand, leaving no demand in queue at the end of the simulation, where as for the `wkm-swkm` solution method leaves demand in queue for all other flight configurations.

```{r, echo=F, message=F, fig.dim=c(8,3)*1.1}
york |> 
  filter(`Number of UAVs` == "high",
         `Queue strategy` == "FCFS") |>
  ggplot(aes(x = `Arrival rate variance`, y = `Demands in queue`, fill = `Flight configuration`)) +
  geom_col(position = position_dodge(.9)) +
  facet_grid(`Queue strategy`~`Solution method`, scales = "free", labeller = label_both)
```

Lastly we can look at the Ploss for the no queue strategy.

```{r, echo=F, message=F, fig.dim=c(8,3)*1.1}
york |> 
  filter(`Number of UAVs` == "high",
         `Queue strategy` == "No queue") |>
  ggplot(aes(x = `Arrival rate variance`, y = Ploss, fill = `Flight configuration`)) +
  geom_col(position = position_dodge(.9)) +
  facet_grid(`Queue strategy`~`Solution method`, scales = "free", labeller = label_both)
```

This also shows that the only flight configuration that is able to serve all demand is the free-flight: no constraint.

Even though we have seen here the good performance for free-flight: no constraint, we would expect that this is also because of the low utilization of UAVs overall in the system.
Had there been more strain on the system, then free-flight: no constraint would likely not have had the time to fly the long distances and thus have served all demand.
In fact the ability to fly long distances to serve demand could be detrimental to overall performance if the UAV does not take the distance into consideration as part of the assignment policy.

## Showing mean response

```{r}
york <- readRDS("inst/extdata/york_simulation_metadata.rds")

york <- york %>%
  mutate(`Flight configuration` =
           factor(as.character(`Flight configuration`), 
                  levels = c("Zoned", "Free-flight: 0", "Free-flight: no constraint"),
                  labels = c("Zoning", "Maximum flight range", "Free-flight")))

york |> 
  filter(
    `Solution method` == "wkm-swkm",
    `Arrival rate variance` == "low",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
  ) |> 
  # tidyr::unnest(cols = response_times) |>
  ggplot() +
  geom_boxplot(aes(y = `Mean response`/60, color = `Flight configuration`)) +
  facet_wrap(
     ~ `Queue strategy`, 
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  labs(y = "Response time (minutes)") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

ggsave("mean_response_york.pdf", width = 8, height = 3)
```

## Showing Ploss

```{r}
york |> 
  filter(
    `Arrival rate variance` == "low",
    `Solution method` == "ga-tot",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight"),
    `Queue strategy`== "No queue"
  ) |> 
  ggplot() +
  geom_boxplot(aes(y = `Ploss`, color = `Flight configuration`)) +
  facet_grid(
     ~ `Queue strategy`, 
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  labs(y = "")

ggsave("Ploss_york.pdf", width = 5, height = 2.5)
```

## Showing distance

```{r}
york |> 
  filter(
    `Arrival rate variance` == "low",
    `Solution method` == "wkm-swkm",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
  ) |> 
  ggplot() +
  geom_boxplot(aes(y = `Minimum distance`*110, color = `Flight configuration`)) +
  facet_grid(
     ~ `Queue strategy`, 
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  labs(y = "Minimum distance") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

ggsave("minimum_distance_york.pdf", width = 8, height = 3)
```


```{r}
york |> 
  filter(
    `Solution method` == "wkm-swkm",
    `Arrival rate variance` == "high",
    `Number of UAVs` == "high",
    `Flight configuration` %in% c("Zoning", "Maximum flight range", "Free-flight")
  ) |> 
  ggplot() +
  geom_boxplot(aes(y = `5th percentile distance`*110, color = `Flight configuration`)) +
  facet_grid(
     ~ `Queue strategy`, 
    labeller = label_both, 
    scales = "free"
  ) +
  #scale_color_grey() +
  theme_bw() +
  labs(y = "5th percentile distance") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

ggsave("5th_minimum_distance_york.pdf", width = 8, height = 3)
```