library(ggplot2)
library(patchwork)
library(dplyr)
library(tidyverse)
library(scales)
library(stringr)
library(tidymodels)

formula_stats <- read_delim('./reference-values.csv') %>% mutate(instance = sub('\\.cnf', '', instance),
                                                             mus_instance = str_detect(instance, 'mus'),
                                                             instance_class = sub('-.*', '', instance),
                                                             clique = str_detect(instance, 'clique'),
                                                             full_instance_class = ifelse(
                                                               instance_class == 'coloring',
                                                               ifelse(clique, 'coloring-clique', 'coloring'),
                                                               instance_class
                                                             ))

p <- ggplot(formula_stats %>% filter(mus_instance, instance_class == '3cnf'), aes(x = formula_size)) +
  geom_histogram(binwidth = 1, color = 'black', fill = 'white') +
  labs(x = 'Formula size', y = 'Number of samples')
p
ggsave('full-distribution-3cnf-mus.svg', plot = p, width = 3, height = 3)

p <- ggplot(formula_stats %>% filter(!mus_instance, instance_class == '3cnf'), aes(x = formula_size)) +
  geom_histogram(binwidth = 4, color = 'black', fill = 'white') +
  labs(x = 'Formula size', y = 'Number of samples')
p
ggsave('full-distribution-3cnf-non-mus.svg', plot = p, width = 3, height = 3)

df_mus <- read_delim('./results.mus.csv') %>%
  mutate(solver = sub('_.*', '', run_name),
         instance = sub('.*_', '', run_name),
         optimal = proof_length == proof_bound,
         mus_instance = TRUE)

df_orig <- read_delim('./results.orig.csv') %>%
  mutate(solver = sub('_.*', '', run_name),
         instance = sub('.*_', '', run_name),
         optimal = proof_length == proof_bound,
         gap_to_mus_bound = NA,
         mus_instance = FALSE)

df <- bind_rows(df_mus, df_orig)

df %>% filter(!tracecheck_result) %>% any()

bugs <- df %>% group_by(instance) %>% filter(optimal) %>%
  summarise(max_proof_length = max(proof_length),
            min_proof_length = min(proof_length)) %>%
  filter(max_proof_length != min_proof_length)

df %>% group_by() %>% summarize(n = length(unique(instance)))
bugs %>% count

df_bugs <- df |> inner_join(bugs) %>% filter(proof_length > min_proof_length)
df_bugs

df_dpll <- df %>% filter(optimal) %>% filter(solver == 'proofmin-dpll') |> inner_join(formula_stats)
df_cadical <- df %>% filter(optimal) %>% filter(solver == 'proofmin-cadical') |> inner_join(formula_stats)

df_delta <- df_dpll |> full_join(df_cadical, by = c('instance', 'instance_class', 'formula_size'))
p <- ggplot(df_delta, aes(x = proof_length.x, y = proof_length.y, color = instance_class)) +
  geom_point(shape = 4) +
  geom_abline(color = 'red', linetype = 'dashed') +
  scale_x_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x)))
p

p <- ggplot(df_delta, aes(x = formula_size, y = proof_length.y / proof_length.x, color = instance_class)) +
  geom_point(shape = 4) +
  scale_x_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x)))
p

df_with_refs <- df %>% filter(solver == 'proofmin') |> inner_join(formula_stats)
df_with_refs %>% filter(ref_proof_len < proof_length & optimal) %>% count

df_with_refs %>% filter(instance_class == 'php') %>% View
df_with_refs %>% filter(instance_class == 'parity') %>% View
df_with_refs %>% filter(instance_class == 'op') %>% View

p <- ggplot(df_with_refs %>% filter(!is.na(proof_length) & optimal),
            aes(x = formula_size, y = ref_proof_len / proof_length)) +
  geom_point(shape = 4)
p

p <- ggplot(df_with_refs %>% filter(!is.na(proof_length) & optimal),
             aes(x = ref_proof_len, y = proof_length, color = formula_size)) +
  geom_point() +
  geom_abline(color = 'red', linetype = 'dashed') +
  geom_abline(slope = 1/1.5, color = 'red', linetype = 'dashed') +
  annotate('text', x=25, y=30, label = 'CaDiCaL proof length = optimal proof length',
           angle = 65, color = 'navy') +
  annotate('text', x=42, y=25, label = 'CaDiCaL proof length = 1.5 × optimal proof length',
           angle = 55, color = 'navy') +
  scale_color_gradient(low = "#FFDC00", high = "#FF0000",) + 
  scale_x_continuous(expand = c(0, NA)) + scale_y_continuous(expand = c(0, NA)) + 
  theme(legend.position = c(.15, .85), legend.background=element_blank()) +
  labs(color = 'Formula size',
       x = 'Proof length from CaDiCaL', y = 'Optimal proof length', shape = 'Minimally UNSAT?') +
  coord_cartesian(xlim = c(0, 150), ylim = c(0, 75))
p
ggsave('cadical-vs-optimal.svg', plot = p, width = 6, height = 3)

p <- ggplot(df_with_refs %>% filter(!is.na(proof_length) & instance_class == '3cnf'),
            aes(x = ref_proof_len, y = proof_length, color = formula_size)) +
  geom_point() +
  geom_abline(color = 'red', linetype = 'dashed') +
  geom_abline(slope = 1/1.5, color = 'red', linetype = 'dashed') +
  annotate('text', x=60, y=70, label = 'CaDiCaL proof length = optimal proof length', angle = 57.5) +
  annotate('text', x=100, y=60, label = 'CaDiCaL proof length = 1.5 × optimal proof length', angle = 45) +
  scale_color_gradient(low = "#56B1F7", high = "#132B43",) + 
  theme(legend.position = c(.15, .55), legend.background=element_blank()) +
  labs(color = 'Formula size',
       x = 'Proof length from CaDiCaL', y = 'Best proof length', shape = 'Minimally UNSAT?')
p
ggsave('cadical-vs-best-3cnf.svg', plot = p, width = 8, height = 8)

p <- ggplot(df_with_refs %>% filter(!is.na(proof_length) & proof_length < 1e5 & ref_proof_len < 1e5),
            aes(x = ref_proof_len, y = proof_length, color = instance_class)) +
  geom_point() +
  geom_abline(color = 'red', linetype = 'dashed') +
  scale_x_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))) +
  theme(legend.position = c(.15, .55), legend.background=element_blank()) +
  labs(color = 'Formula type',
       x = 'Proof length from Kissat', y = 'Best proof length', shape = 'Minimally UNSAT?')
p
ggsave('kissat-vs-best.svg', plot = p, width = 8, height = 8)

p <- ggplot(df_with_refs %>% filter(!is.na(proof_length)),
       aes(x = proof_length / ref_proof_len, y = full_instance_class, color = mus_instance)) +
  geom_boxplot() +
  geom_vline(xintercept=1, linetype="dashed", color = "red") +
  scale_x_log10() +
  theme(legend.position = c(.8, .2),
        legend.background=element_blank(),
        axis.text.y = element_text(angle = 45, vjust = 0.5, hjust=1)) +
  labs(x = 'Proof length relative to Kissat', y = 'Instance class', color = 'MUS instance')
p
ggsave('bnb-vs-kissat-ratios.svg', plot = p, width = 5, height = 2)


df_mus %>% filter(optimal) %>% group_by(instance) %>% slice(1) %>% group_by(gap_to_mus_bound) %>% count()

df_baseline <- df_mus %>% filter(optimal) %>% filter(solver == 'baseline') |> inner_join(formula_stats)
df_our <- df_mus %>% filter(optimal) %>% filter(solver == 'proofmin') |> inner_join(formula_stats)

df_our <- df_baseline |> full_join(df_our, by = c('instance', 'gap_to_mus_bound', 'formula_size')) %>%
  mutate(rel_time = duration.y / duration.x)

df_our |> inner_join(formula_stats %>% select(instance, instance_class)) %>%
  group_by(instance_class, is.na(duration.x), is.na(duration.y)) %>% count()

df_our %>% filter(is.na(duration.y) & !is.na(duration.x)) %>% View

df_our %>% filter(rel_time > 1) %>% View

p <- ggplot(df_our %>% filter(!is.na(rel_time)),
            aes(y = rel_time, x = formula_size)) +
  geom_point(shape = 4) +
  geom_hline(yintercept=1, linetype="dashed", color = "red") +
  labs(x = 'Formula size') +
  scale_y_log10("Relative solving time",
                breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) + 
  theme(legend.position = c(.8, .9),
        legend.background=element_blank())
p
ggsave('relative-time.svg', plot = p, width = 3, height = 4)

p <- ggplot(df_our, aes(x = factor(gap_to_mus_bound), y = duration.y)) +
  geom_boxplot(show.legend = FALSE, color='darkgray') +
  geom_jitter(size = 0.75, shape = 4, color = 'black') +
  scale_y_log10("Time to solve with branch-and-bound, s",
                breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) + 
  labs(
    x = 'Gap to MUS bound',
    y = 'Time to solve with branch-and-bound, s'
  ) +
  coord_flip()
p
ggsave('bnb-time-gap.svg', plot = p, width = 5, height = 3)

p <- ggplot(df_our %>% filter(gap_to_mus_bound > 0),
            aes(x = factor(formula_size), y = duration.y, color = factor(gap_to_mus_bound))) +
  geom_boxplot(show.legend = FALSE, color='darkgray', outliers = FALSE) +
  geom_jitter(size = 0.75, shape = 4) +
  scale_y_log10("Time to solve with branch-and-bound, s",
                breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) + 
  labs(
    x = 'Formula size',
    y = 'Time to solve with branch-and-bound, s',
    color = 'Gap to MUS bound'
  ) +   theme(legend.position='top')
p
ggsave('bnb-time-clauses-nonzero.svg', plot = p, width = 5, height = 3)

p <- ggplot(df_our %>% filter(gap_to_mus_bound == 0),
            aes(x = factor(formula_size), y = duration.y)) +
  geom_boxplot(show.legend = FALSE, color='darkgray', outliers = FALSE) +
  geom_jitter(size = 0.75, shape = 4) +
  scale_y_log10("Time to solve with branch-and-bound, s",
                breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) + 
  labs(
    x = 'Formula size',
    y = 'Time to solve with branch-and-bound, s',
  ) +   theme(legend.position='top')
p
ggsave('bnb-time-clauses-zero.svg', plot = p, width = 5, height = 3)

