---
title: "Supplementary Report"
output:
  pdf_document:
    highlight: zenburn
    number_sections: yes
    toc: yes
  html_document:
    code_folding: hide
    toc: yes
    toc_depth: 1
    number_sections: true
params:
  csv_path: "data.csv"
  models_csv_path: "models.csv"
  timeout_ms: 180000
---

This is a supplementary report for the paper **"A Preliminary Analysis on the Effect of Randomness in a CEGAR Framework"** by **Akos Hajdu and Zoltan Micskei**, presented at the **25th PhD Mini-Symposium** (2018), organized by the _Department of Measurement and Information Systems_ at the _Budapest University of Technology and Economics_.

# Initialization

```{r includes, message = F}
library(tidyverse)
library(stringr)
library(gridExtra)
library(GGally)

theme_rotate_x <- theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
theme_noticks <- theme(axis.ticks = element_blank())
fill_heat <- scale_fill_gradient(low = "green", high = "red", na.value = "white")
percent <- function(actual, total, round = 0) { round(actual / total * 100, digits = round) }
```

## Load data

```{r}
d <- read_csv(
  params$csv_path,
  col_types = cols(
    Model = col_character(),
    Domain = col_character(),
    Refinement = col_character(),
    Search = col_character(),
    RandomVars = col_character(),
    Result = col_character(),
    TimeMs = col_integer(),
    Iterations = col_integer(),
    ArgSize = col_integer(),
    ArgDepth = col_integer(),
    ArgMeanBranchFactor = col_double(),
    CexLen = col_integer()
  )
)

d_model_attrs <- read_csv(
  params$models_csv_path,
  col_types = cols(
    Model = col_character(),
    Inputs = col_integer(),
    Latches = col_integer(),
    AndGates = col_integer(),
    Vars = col_integer(),
    Locations = col_integer(),
    Edges = col_integer(),
    Assigns = col_integer(),
    Assumes = col_integer(),
    Havocs = col_integer()
  )
)
```

## Clean data

```{r}
d$Domain <- factor(d$Domain)
d$Refinement <- gsub("BW_BIN_ITP", "BIN", d$Refinement)
d$Refinement <- gsub("SEQ_ITP", "SEQ", d$Refinement)
d$Refinement <- factor(d$Refinement)
d$RandomVars <- as.logical(d$RandomVars)

# Combine RandomVars and Search into a single column
d <- d %>% mutate(Randomized = ifelse(RandomVars, "VARS", ifelse(Search == "RND", "SEARCH", "DET")))
d <- d  %>% select(-Search, -RandomVars)
d$Randomized <- factor(d$Randomized)

# Create 'Config' column
d <- d %>% mutate(Config = paste(
  substr(Domain, 1, 1),
  substr(Refinement, 1, 1),
  "-",
  substr(Randomized, 1, 1),
  sep = ""
))
d$Config <- factor(d$Config)

all_params <- c("Randomized", "Domain", "Refinement")

# Convert 'Result' to logical (and remove exceptions)
n_rows_empty <- sum(is.na(d$Result))
n_rows_exception <- nrow(d %>% filter(str_detect(d$Result, "\\[EX\\]")))
d$Result <- as.logical(d$Result)

# Shorten model names by removing irrelevant details
simplify_name <- function(name) {
  name <- gsub("../models/cfa/", "", name)
  name <- gsub("../models/sts/", "", name)
  name <- gsub("\\.cfa", "", name)
  name <- gsub("\\.aag", "", name)
  name <- gsub("hwmcc", "hw", name)
  name <- gsub("-unreach-call", "", name)
  name <- gsub("_true", "", name)
  name <- gsub("_false", "", name)
  name <- gsub("\\.c", "", name)
  name <- gsub("Problem0", "prob", name)
  name <- gsub("label", "lab", name)
  name <- gsub("(ssh/.*)_0", "\\1", name)
  name <- gsub("(eca/.*)_0", "\\1", name)
  name
}
d$Model <- simplify_name(d$Model)
d_model_attrs$Model <- simplify_name(d_model_attrs$Model)

# Determine category
d <- d %>% mutate(Category = sub("/([^/]*)$","", d$Model))
d$Model <- as.factor(d$Model)

str(d)
```

# Overview

## Summary

```{r}
n_meas_rep_total <- nrow(d)
n_meas_rep_succ <- nrow(d %>% filter(!is.na(Result)))
n_models_total <- length(unique(d$Model))
n_configs_total <- length(unique(d$Config))
categories <- unique(d$Category)
n_categories <- length(categories)
n_rows_problem <- n_rows_empty + n_rows_exception

d_model_config <- d %>%
  group_by_(.dots = c("Config", "Model", "Category", all_params)) %>%
  summarise(
    Reps = n(),
    SuccCount = sum(!is.na(Result)),
    SuccSD = sd(!is.na(Result)),
    Succ = SuccCount > 0,
    TimeMsAvg = mean(TimeMs, na.rm = T),
    TimeMsRSD = ifelse(SuccCount == 1, 0, sd(TimeMs, na.rm = T) / TimeMsAvg),
    IterationsAvg = mean(Iterations, na.rm = T),
    IterationsRSD = ifelse(SuccCount == 1, 0, sd(Iterations, na.rm = T) / IterationsAvg),
    ArgSizeAvg = mean(ArgSize, na.rm = T),
    ArgSizeRSD = ifelse(SuccCount == 1, 0, sd(ArgSize, na.rm = T) / ArgSizeAvg),
    ArgDepthAvg = mean(ArgDepth, na.rm = T),
    ArgDepthRSD = ifelse(SuccCount == 1, 0, sd(ArgDepth, na.rm = T) / ArgDepthAvg),
    CexLenRSD = ifelse(sum(!is.na(CexLen)) == 0, NA, ifelse(sum(!is.na(CexLen)) == 1, 0, sd(CexLen, na.rm = T) / mean(CexLen, na.rm = T)))
    )

n_meas_total <- nrow(d_model_config)
n_meas_succ <- nrow(d_model_config %>% filter(Succ))

if (n_meas_total !=  n_models_total * n_configs_total) {
  warning(sprintf("%d models * %d configs = %d expected measurements, but there are %d actual measurements",
                  n_models_total, n_configs_total, n_models_total * n_configs_total, n_meas_total))
}
d_model <- d_model_config %>% group_by_("Model", "Category") %>%
  summarise(
    SuccRunCount = sum(SuccCount),
    SuccCount = sum(Succ),
    Succ = SuccCount > 0)

n_models_succ <- nrow(d_model %>% filter(Succ))

d_config <- d_model_config %>%
  group_by_(.dots = c("Config", all_params)) %>%
  summarise(
    SuccRunCount = sum(SuccCount),
    SuccCount = sum(Succ))
```

**Terminology**

- A _run_ is a single execution of the algorithm _configuration_ on an input _model_.
    - A run is _successful_ if it terminates within the time limit.
- A _measurement_ is the collection of all repeated runs of a configuration on a model.
    - A measurement is _successful_ if it includes at least one successful run. In this case we also say that the configuration _verified_ the model.

**Results**

- There are **`r n_models_total` models** (from `r n_categories` categories) and **`r n_configs_total` configurations**, giving **`r n_meas_total` measurements**.
    - Each measurement was **repeated `r paste(unique(d_model_config$Reps), sep = "/")` times**, giving a total number of **`r n_meas_rep_total` runs**.
- **`r n_meas_rep_succ` / `r n_meas_rep_total` runs (`r percent(n_meas_rep_succ, n_meas_rep_total)`%) are successful** with a timeout of **timeout of `r params$timeout_ms/1000` s**.
    - **`r n_meas_succ` / `r n_meas_total` measurements (`r percent(n_meas_succ, n_meas_total)`%) are successful**.
    - There are **`r n_rows_problem` problems** (see details later).
- **`r n_models_succ` / `r n_models_total` models (`r percent(n_models_succ, n_models_total)`%) were verified** by at least one configuration.

## Problems

- There are **`r n_rows_empty` empty** results (no timeout, no exception). This usually happens when the SMT solver crashes.

- There are **`r n_rows_exception` exceptions**.

```{r}
if (n_rows_exception > 0) {
  knitr::kable(d_exceptions %>% group_by(Type = Result) %>% summarise(Configs = toString(unique(Config)), Count = n()))
}
```

## Models

The **`r n_models_total` models** are distributed in the **`r n_categories` categories** in the following way.

```{r, fig.height = 1.5, fig.width = 4.5}
ggplot(d_model) +
  geom_bar(aes(x = Category, fill = Category)) +
  coord_flip() + guides(fill = F) +
  scale_fill_brewer(palette = "Set1") +
  ggtitle("Distribution of models among the categories")
```

The following plots show for each model the number of successful runs and the number of configurations that verified them.

```{r, fig.height = 3.5, fig.width = 10}
grid.arrange(
  ggplot(d_model) +
    geom_point(aes(reorder(Model, -SuccRunCount), SuccRunCount, color = Category), size = 2) +
    theme_rotate_x + scale_color_brewer(palette = "Set1") +
    labs(x = "Model", y = "Successful runs", title = "Number of successful runs for the models"),
  ggplot(d_model) +
    geom_point(aes(reorder(Model, -SuccCount), SuccCount, color = Category), size = 2) +
    theme_rotate_x + scale_color_brewer(palette = "Set1") +
    coord_cartesian(ylim = c(0, n_configs_total)) +
    labs(x = "Model", y = "Verified by", title = "Number of configs that verified the models"),
  ncol = 2)
```

The following table summarizes the attributes of the models. Hardware models can be described by the number of hardware elements, while software and PLC models can be described by the number of control locations and statements.

```{r}
knitr::kable(d_model_attrs %>% select(Model, Inputs, Latches, AndGates) %>% filter(!is.na(Inputs)))
knitr::kable(d_model_attrs %>% select(Model, Locations, Edges, Assigns, Assumes, Havocs) %>% filter(!is.na(Locations)))
```


## Configurations

Configurations are listed in the following table with their abbreviation, parameters, successful runs and verified models.

```{r}
knitr::kable(d_config)
```

The following plots show the number of models verified and the successful runs by each config. The number of models verified ranges between `r min(d_config$SuccCount)` and `r max(d_config$SuccCount)`, while the number of successful runs is between `r min(d_config$SuccRunCount)` and `r max(d_config$SuccRunCount)`.

```{r, fig.height = 3, fig.width = 10}
grid.arrange(
  ggplot(d_config) +
    geom_bar(aes(reorder(Config, SuccCount), SuccCount, fill = Randomized), stat = "identity") +
    theme_rotate_x + scale_fill_brewer(palette = "Set1") +
    labs(x = "Config", y = "Models verified", title = "Number of models verified by each config"),
  ggplot(d_config) +
    geom_bar(aes(reorder(Config, SuccRunCount), SuccRunCount, fill = Randomized), stat = "identity") +
    theme_rotate_x + scale_fill_brewer(palette = "Set1") +
    labs(x = "Config", y = "Succ. runs", title = "Number of successful runs by each config"),
  ncol = 2)
```

The following plot compares the number of verified models to the number of successful runs for each configuration.

```{r, fig.height = 3, fig.width = 4}
ggplot(d_config) +
  geom_point(aes(SuccCount, SuccRunCount, color = Randomized)) +
  scale_color_brewer(palette = "Set1") +
  labs(x = "Number of models verified", y = "Number of succ. runs")
```


## Output metrics

The following histograms show an overview on the range of the output metrics colored by the category.

```{r, fig.width = 10, fig.height = 7}
grid.arrange(
  ggplot(d) + geom_bar(aes(!is.na(Result), fill = Category)) + xlab("Succ") + scale_fill_brewer(palette = "Set1"),
  ggplot(d) + geom_histogram(aes(TimeMs, fill = Category), na.rm = T) + scale_x_log10() + scale_fill_brewer(palette = "Set1"),
  ggplot(d) + geom_histogram(aes(Iterations, fill = Category), na.rm = T) + scale_x_log10() + scale_fill_brewer(palette = "Set1"),
  ggplot(d) + geom_histogram(aes(ArgSize, fill = Category), na.rm = T) + scale_x_log10() + scale_fill_brewer(palette = "Set1"),
  ggplot(d) + geom_histogram(aes(ArgDepth, fill = Category), na.rm = T) + scale_x_log10() + scale_fill_brewer(palette = "Set1"),
  ggplot(d) + geom_histogram(aes(CexLen, fill = Category), na.rm = T) + scale_x_log10() + scale_fill_brewer(palette = "Set1"),
  ncol = 2)
```

The following scatterplot matrix shows the numerical output metrics. 

```{r, fig.width = 10, fig.height = 10}
ggpairs(d %>% filter(!is.na(Result)) %>% select(TimeMs, Iterations, ArgSize, ArgDepth, CexLen))
```


# Success rates (RQ1)

The following heatmaps show for each measurement whether it was successful (at least one successful run) and the number of successful executions (if there were any). In the latter plot, those measurements are marked (with circles for randomized search and triangles for randomized variable names) where a randomized configuration verified a model that its deterministic counterpart could not.

```{r, fig.height = 4, fig.width = 8}
d_model_config2 <- d_model_config %>% ungroup() %>%
  select(Config, Model, Randomized, Domain, Refinement, SuccCount)
d_model_config2 <- inner_join(
  d_model_config2,
  d_model_config2 %>% filter(Randomized == "DET"),
  by = c("Model", "Domain", "Refinement")) %>%
  mutate(Marked = ifelse((SuccCount.x > 0) & (SuccCount.y == 0), ifelse(Randomized.x == "VARS", "V", "S"), NA))

d_model_config$Marked <- d_model_config2$Marked

grid.arrange(
  ggplot(d_model_config, aes(Config, Model)) +
    geom_tile(aes(fill = Succ), color = "black") +
    theme_noticks + theme_rotate_x +
    ggtitle("At least one successful run"),
  ggplot(d_model_config, aes(Config, Model)) +
    geom_tile(aes(fill = ifelse(SuccCount > 0, SuccCount, NA)), color = "black") +
    geom_point(aes(shape = Marked), na.rm = T) + guides(shape=F) +
    theme_noticks + theme_rotate_x +
    scale_fill_gradient(low = "red", high = "green", na.value = "white", limits = c(1, NA), breaks = c(1, 10, 20, 30)) +
    labs(title = "Number of successful runs", fill = "Count"),
  ncol = 2)
```

```{r}
det_succ <- unique((d %>% filter(!is.na(Result)) %>% filter(Randomized == "DET"))$Model)
rnd_succ <- unique((d %>% filter(!is.na(Result)) %>% filter(Randomized != "DET"))$Model)
det_only <- setdiff(det_succ, rnd_succ)
rnd_only <- setdiff(rnd_succ, det_succ)
```

There are `r length(det_only)` models that were only verified by deterministic configs: `r paste(det_only)`.

There are `r length(rnd_only)` models that were only verified by randomized configs: `r paste(rnd_only)`.

# Variance of the output metrics (RQ2)

For each measurement, there are multiple runs, forming a distribution over the output metrics.

## Individual distributions

The following plots show the distributions of the output metrics for each measurement individually.

```{r, fig.height = 15, fig.width = 12}
ggplot(d) +
  geom_bar(aes(!is.na(Result), fill = Randomized)) +
  facet_grid(Model ~ Config) +
  theme_bw() + theme(strip.text.y = element_text(angle = 0)) + theme_rotate_x +
  scale_fill_brewer(palette = "Set1") +
  labs(x = "Succ", title = "Distributions of successful executions")
```

```{r, fig.height = 15, fig.width = 12}
ggplot(d) +
  geom_histogram(aes(TimeMs, fill = Randomized), bins = 15, na.rm = T) +
  facet_grid(Model ~ Config) +
  theme_bw() + theme(strip.text.y = element_text(angle = 0)) + theme_rotate_x +
  scale_x_log10() + scale_fill_brewer(palette = "Set1") +
  ggtitle("Distributions of execution time")
```

```{r, fig.height = 15, fig.width = 12}
ggplot(d) +
  geom_histogram(aes(Iterations, fill = Randomized), bins = 15, na.rm = T) +
  facet_grid(Model ~ Config) +
  theme_bw() + theme(strip.text.y = element_text(angle = 0)) + theme_rotate_x +
  scale_x_log10() + scale_fill_brewer(palette = "Set1") +
  ggtitle("Distributions of iterations")
```

```{r, fig.height = 15, fig.width = 12}
ggplot(d) +
  geom_histogram(aes(ArgSize, fill = Randomized), bins = 15, na.rm = T) +
  facet_grid(Model ~ Config) +
  theme_bw() + theme(strip.text.y = element_text(angle = 0)) + theme_rotate_x +
  scale_x_log10() + scale_fill_brewer(palette = "Set1") +
  ggtitle("Distributions of ARG size")
```

```{r, fig.height = 15, fig.width = 12}
ggplot(d) +
  geom_histogram(aes(ArgDepth, fill = Randomized), bins = 15, na.rm = T) +
  facet_grid(Model ~ Config) +
  theme_bw() + theme(strip.text.y = element_text(angle = 0)) + theme_rotate_x +
  scale_x_log10() + scale_fill_brewer(palette = "Set1") +
  ggtitle("Distributions of ARG depth")
```

```{r, fig.height = 15, fig.width = 12}
ggplot(d) +
  geom_histogram(aes(CexLen, fill = Randomized), bins = 15, na.rm = T) +
  facet_grid(Model ~ Config) +
  theme_bw() + theme(strip.text.y = element_text(angle = 0)) + theme_rotate_x +
  scale_x_log10() + scale_fill_brewer(palette = "Set1") +
  ggtitle("Distributions of cex. length")
```

## Individual deviations

In this research question we are mainly interested in the **deviations** of the distributions. The following heatmaps summarize the deviations for each measurement and each output metric. The _SuccSD_ variable shows the _standard deviation_ of the success of the verification (0: false, 1: true). Other metrics correspont to the _relative standard deviation_ (RSD = SD / mean) of the output metrics.

```{r, fig.height = 10, fig.width = 8}
grid.arrange(
  ggplot(d_model_config, aes(Config, Model)) +
    geom_tile(aes(fill = SuccSD), color = "black") +
    theme_noticks + theme_rotate_x +
    fill_heat + ggtitle("Success SD"),
  
  ggplot(d_model_config, aes(Config, Model)) +
    geom_tile(aes(fill = TimeMsRSD), color = "black") +
    theme_noticks + theme_rotate_x +
    fill_heat + ggtitle("Time RSD"),
  
  ggplot(d_model_config, aes(Config, Model)) +
    geom_tile(aes(fill = IterationsRSD), color = "black") +
    theme_noticks + theme_rotate_x +
    fill_heat + ggtitle("Iterations RSD"),
  
  ggplot(d_model_config, aes(Config, Model)) +
    geom_tile(aes(fill = ArgSizeRSD), color = "black") +
    theme_noticks + theme_rotate_x +
    fill_heat + ggtitle("ARG size RSD"),
  
  ggplot(d_model_config, aes(Config, Model)) +
    geom_tile(aes(fill = ArgDepthRSD), color = "black") +
    theme_noticks + theme_rotate_x +
    fill_heat + ggtitle("ARG depth RSD"),
  
  ggplot(d_model_config, aes(Config, Model)) +
    geom_tile(aes(fill = CexLenRSD), color = "black") +
    theme_noticks + theme_rotate_x +
    fill_heat + ggtitle("Counterexample length RSD"),
ncol = 2)
```

## Summarized deviations

The following plots summarize the distribution of the deviations, grouped by the parameter that is randomized.

```{r, fig.height = 3, fig.width = 10}
grid.arrange(
  ggplot(d_model_config) + geom_boxplot(aes(x = Randomized, y = SuccSD), na.rm = T) + theme_rotate_x,
  ggplot(d_model_config) + geom_boxplot(aes(x = Randomized, y = TimeMsRSD), na.rm = T) + theme_rotate_x,
  ggplot(d_model_config) + geom_boxplot(aes(x = Randomized, y = IterationsRSD), na.rm = T) + theme_rotate_x,
  ggplot(d_model_config) + geom_boxplot(aes(x = Randomized, y = ArgSizeRSD), na.rm = T) + theme_rotate_x,
  ggplot(d_model_config) + geom_boxplot(aes(x = Randomized, y = ArgDepthRSD), na.rm = T) + theme_rotate_x,
  ggplot(d_model_config) + geom_boxplot(aes(x = Randomized, y = CexLenRSD), na.rm = T) + theme_rotate_x,
  ncol = 6)
```

The following plot shows the same data, but using a uniform y axis scale.

```{r, fig.height = 4}
d_model_config2 <- d_model_config %>%
  gather("Metric", "Deviation", SuccSD, TimeMsRSD, IterationsRSD, ArgSizeRSD, ArgDepthRSD, CexLenRSD)
ggplot(d_model_config2) +
  geom_boxplot(aes(Metric, Deviation, color = Randomized), na.rm = T)
```

The following plot shows the same data, but outliers above 1.2 are trimmed (there are `r nrow(d_model_config2 %>% filter(Deviation >= 1.2))` of them).

```{r, fig.height = 4}
ggplot(d_model_config2 %>% filter(Deviation < 1.2)) +
  geom_boxplot(aes(Metric, Deviation, color = Randomized), na.rm = T)
```

The following parallel coordinate plots show the deviations of each output metric, where each line corresponds to a measurement.

```{r, fig.width = 10, fig.height = 3}
grid.arrange(
  ggparcoord(d_model_config %>% ungroup() %>% filter(Randomized == "SEARCH") %>%
               select(TimeMsRSD, IterationsRSD, ArgSizeRSD, ArgDepthRSD, CexLenRSD)) +
    ggtitle("RSDs for randomized search"),
  ggparcoord(d_model_config %>% ungroup() %>% filter(Randomized == "VARS") %>%
               select(TimeMsRSD, IterationsRSD, ArgSizeRSD, ArgDepthRSD, CexLenRSD)) +
    ggtitle("RSDs for randomized variables"),
  ncol = 2)
```


The following table shows the mean values of the deviations for each output metric grouped by the parameter that is randomized.

```{r}
knitr::kable(d_model_config %>% group_by(Randomized) %>% summarize(
  Succ = mean(SuccSD, na.rm = T),
  TimeMs = mean(TimeMsRSD, na.rm = T),
  Iterations = mean(IterationsRSD, na.rm = T),
  ArgSize = mean(ArgSizeRSD, na.rm = T),
  ArgDepth = mean(ArgDepthRSD, na.rm = T),
  CexLen = mean(CexLenRSD, na.rm = T)
), caption = "Mean value of the deviations")
```

The following table shows the mean values of the deviations for each output metric for each configuration.

```{r}
knitr::kable(d_model_config %>% group_by(Config) %>% summarize(
  Succ = mean(SuccSD, na.rm = T),
  TimeMs = mean(TimeMsRSD, na.rm = T),
  Iterations = mean(IterationsRSD, na.rm = T),
  ArgSize = mean(ArgSizeRSD, na.rm = T),
  ArgDepth = mean(ArgDepthRSD, na.rm = T),
  CexLen = mean(CexLenRSD, na.rm = T)
), caption = "Mean value of the deviations")
```