plotting.Rmd

---
title: "Figure making"
output: 
  github_document:
    toc: true
    toc_depth: 5
date: "2025-10-21"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Load in functions and packages
```{r eval=F, message = F}
packages = c("openxlsx", "here") 
installed = packages %in% installed.packages()[, "Package"]
if (any(!installed)) {
  install.packages(packages[!installed])
}
lapply(packages, library, character.only = TRUE)

source(here("plotting/plotting-functions.R")) # this calls the other scripts
```

## Load in data summaries and other global parameters
```{r eval=F, message = F}
N = 100000

cb.data = read.xlsx(here("data/summaries/critical-bubble-underdominance-alleles.xlsx")) %>% 
  as_tibble()

under.data = read.xlsx(here("data/summaries/rectangle1D-underdominance-alleles.xlsx")) %>% 
  as_tibble() 

gd.data = read.xlsx(here("data/summaries/rectangle1D-homing-gene-drive.xlsx")) %>%
  as_tibble() 

tade.mod1D.data = read.xlsx(here("data/summaries/rectangle1D-tade-modification.xlsx")) %>%
  as_tibble()

tade.mod2D.data = read.xlsx(here("data/summaries/circle2D-tade-modification.xlsx")) %>%
  as_tibble()

D.labels = c(TeX(r"($10^{-5}$)"),TeX(r"($10^{-4}$)"),TeX(r"($10^{-3}$)"), TeX(r"($10^{-2}$)"))

text_size = 20
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))


linewidth = 0.75 # default
min_pt_size = 3 # for D sizing
max_pt_size = min_pt_size*4
```

## 1. Factors that influence the critical bubble

```{r eval=F}
t.sizes = calculate_text_size(artboard_width = 921, overleaf_scaling = 0.6)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))

linewidth = 1

#### panel A: sqrt(D/s) 
sqrt.terms = c(0.001, 0.01, 0.1) # 3 magnitudes of sqrt(D/s)
D = 1e-5 
baseline.phat = 0.2 
s.is = D/sqrt.terms^2

bubble1 = critical_bubble_vector_underdominance(phat = baseline.phat, D = D, 
                                        s = s.is[1],
                                        center_0 = T) # sqrt.term = 0.001
bubble2 = critical_bubble_vector_underdominance(phat = baseline.phat, D = D, s = s.is[2], 
                                        center_0 = T) # sqrt.term = 0.01
bubble3 = critical_bubble_vector_underdominance(phat = baseline.phat, D = D, s = s.is[3], 
                                        center_0 = T) # sqrt.term = 0.1
xlow = -1; xhigh = 1 

#palette = brewer.pal(n = 5, name = "PuBuGn")[c(3:5)]
palette = carto_pal(4, "Emrld")[2:4]
sqrt.labs = c(TeX(r"($10^{-3}$)"), TeX(r"($10^{-2}$)"), TeX(r"($10^{-1}$)"))

subplot.D.s.legend =ggplot(data = bubble1$data, aes(x = x, y = p)) + 
  geom_line(aes(color = "a"), size = linewidth) +
  geom_line(data = bubble2$data, aes(x=x, y=p, color = "b"), size=linewidth) +
  geom_line(data = bubble3$data, aes(x = x, y = p, color = "c"), size = linewidth) +
  xlab(TeX(r"(\textit{$x$})")) + ylab(TeX(r"(\textit{$p$})")) + xlim(xlow, xhigh) +
  ylim(0,1) + 
  theme_minimal() +
  scale_x_continuous(limits = c(-1,1), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0)) +
  scale_colour_manual(name = TeX(r"(\textit{$\sqrt{\frac{D}{s_u}}$})"), values = c("a" = palette[1], "b" = palette[2], "c" = palette[3]), 
                      labels = sqrt.labs) +
  sizes.theme 


panelA = subplot.D.s.legend + 
  theme(legend.position = "none")  +
  theme(axis.ticks.length = unit(0.3, "cm"), axis.ticks.margin = unit(0.3, "cm"))

# legend
# ggsave(plot = subplot.D.s.legend, 
#        filename = paste0(fig.dir,"cb-recapitulation/critical-bubble-equation-factors-subplotA-withlegend.pdf"),
#        width= 6, height = 6)

#### panel B: phat
baseline.sqrtD.over.s=0.02
D = 1e-5
s = D/baseline.sqrtD.over.s^2
phat.vals = c(0.05, 0.2, 0.4)

bubble1 = critical_bubble_vector_underdominance(phat = phat.vals[1], D = D, s = s,
                                        center_0 = T) # phat=0.05
bubble2 = critical_bubble_vector_underdominance(phat = phat.vals[2], D = D, s = s, 
                                        center_0 = T) # phat=0.2
bubble3 = critical_bubble_vector_underdominance(phat = phat.vals[3], D = D, s = s, 
                                        center_0 = T) # phat=0.4
xlow = -1; xhigh = 1 # based on plot 3

subplot.phat.option3.withlegend = ggplot(data = bubble1$data, aes(x = x, y = p)) + 
  geom_line(aes(color = "a"), size = linewidth) +
  geom_line(data = bubble2$data, aes(x=x, y=p, color = "b"), size=linewidth) +
  geom_line(data = bubble3$data, aes(x = x, y = p, color = "c"), size = linewidth) +
  xlab(TeX(r"(\textit{$x$})")) + ylab(TeX(r"(\textit{$p$})")) + xlim(xlow, xhigh) +
  ylim(0,1) + 
  scale_colour_manual(name = TeX(r"(\textit{$\hat{p}$})"), values = c("a" = palette[1], "b" = palette[2], "c" = palette[3]), 
                      labels = c('0.05','0.2','0.4')) +
  theme_minimal() +
  scale_x_continuous(limits = c(-1,1), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0))  +
  sizes.theme


panelB = subplot.phat.option3.withlegend + theme(legend.position = "none") +
  theme(legend.position = "none")  +
  theme(axis.ticks.length = unit(0.3, "cm"), axis.ticks.margin = unit(0.3, "cm"))

# legend
# ggsave(plot = subplot.phat.option3.withlegend,
#        filename = paste0(fig.dir,"cb-recapitulation/critical-bubble-equation-factors-subplotB-withlegend.pdf"),
#        width = 6, height = 6)

#### panel C: phat vs M(phat)
phat.seq = seq(0.0, 0.5, length.out = 1000)
log.term = log_term_in_auc(phat.seq)

p.info = tibble(phat = phat.seq, log.term = log.term)

panelC = p.info %>%
  ggplot(aes(x = phat, y = log.term)) + geom_line(size = linewidth) + 
  ylab(TeX(r"(\textit{$M(\hat{p})$})")) + xlab(TeX(r"(\textit{$\hat{p}$})")) +
  theme_minimal() +
  scale_x_continuous(limits = c(0,0.5), expand = c(0,0),breaks = c(0.0,0.1,0.2,0.3,0.4,0.5)) +
  scale_y_continuous(limits = c(0,5), expand = c(0,0), breaks = c(1.00,2.00,3.00,4.00,5.00),
                     labels = c("1","2","3","4","5")) +
  sizes.theme +
  theme(legend.position = "none")  +
  theme(axis.ticks.length = unit(0.3, "cm"), axis.ticks.margin = unit(0.3, "cm"))

#### panel D: 
ptilda.seq = ptilda_underdominance(phat.seq)

ptilda.info = tibble(phat = phat.seq, ptilda = ptilda.seq)

panelD = ptilda.info %>%
  ggplot(aes(x = phat, y = ptilda)) +
  xlab(TeX(r"(\textit{$\hat{p}})")) +ylab(TeX(r"(\textit{$p_c(\hat{p})$})")) +
  geom_line(size = linewidth) + theme_minimal() +
  scale_x_continuous(limits = c(0,0.5), expand = c(0.0,0)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0),  breaks = c(0.25, 0.5, 0.75, 1.0)) +
  sizes.theme +
  theme(legend.position = "none")  +
  theme(axis.ticks.length = unit(0.3, "cm"), axis.ticks.margin = unit(0.3, "cm"))

dist_width = -0.05
dist_height = -0.18

combined = plot_grid(panelA, NULL, panelB, 
                     NULL, NULL, NULL,
                     panelC, NULL, panelD,
                     nrow = 3, ncol = 3,
                     align = "hv",
                     rel_widths = c(1, dist_width, 1),
                     rel_heights = c(1, dist_height, 1))

# ggsave(plot = combined, filename = paste0(fig.dir,"cb-recapitulation/combined-critical-bubble-factors.pdf"),
#        width = 12, height = 12)
```

## 2. Results confirming the critical bubble threshold in SLiM

```{r eval=F, message=F}
t.sizes = calculate_text_size(artboard_width = 1306.62, overleaf_scaling = 1)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))

options(scipen=999)

# Panel A
param.no = 3 
file = paste0("paramset", param.no, "-cb-probabilistic-compiled.csv")
data = read_csv(here(paste0("data/underdominance-alleles/bubble/", file)))

pred.size = cb.data$full_cb_auc[which(cb.data$param_no == param.no)]*N
s = data$s[1]; phat = data$phat[1]; D = data$D[1]

title = TeX(r"(\textit{$D = 0.00001, \hat{p} = 0.05, s_u = 0.8$})")

panelA = fit.L2.logistic.wrapper(csv = data, predicted_size = pred.size, 
                                        last_size = 3000, line_thickness = linewidth, add_fitting_range_lines = F)$plot +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0), limits = c(0,1),breaks = c(0.25, 0.5, 0.75, 1.0)) +
  theme_minimal() +
  sizes.theme + 
  xlab("Average release size") + ylab("P(increase)")
  
# panel B
param.no = 11
file = paste0("paramset", param.no, "-cb-probabilistic-compiled.csv")
data = read_csv(here(paste0("data/underdominance-alleles/bubble/", file)))

pred.size = cb.data$full_cb_auc[which(cb.data$param_no == param.no)]*N
s = data$s[1]; phat = data$phat[1]; D = data$D[1]

title = TeX(r"(\textit{$D = 0.00001, \hat{p} = 0.05, s_u = 0.05$})")

panelB = fit.L2.logistic.wrapper(csv = data, predicted_size = pred.size, last_size = 3000,
                                        line_thickness = linewidth, add_fitting_range_lines = F)$plot +
  scale_y_continuous(expand = c(0,0), limits = c(0,1), breaks = c(0.25, 0.5, 0.75, 1.0)) +
  scale_x_continuous(expand = c(0,0)) +
  theme_minimal()+ sizes.theme+
  xlab("Average release size") + ylab("P(increase)")

dist = -0.15
combined = plot_grid(panelA, 
                     NULL,
                     panelB,
                     ncol = 1, nrow = 3,
                     align = "hv",
                     rel_heights = c(1, dist, 1))

# ggsave(plot = combined, filename = paste0(fig.dir,"critical-bubble-in-slim/cb-slim-subplotsA-and-B.pdf"),
#        width = 6, height = 12)

# Add panel C
D5 = TeX(r"($10^{-5}$)")

min.alpha = min(cb.data$slim_bubble_div_full_bubble)
max.alpha = max(cb.data$slim_bubble_div_full_bubble)

panelC_legend = cb.data %>% 
  mutate(phat_f = as.factor(phat),
         s_f = as.factor(s),
         D_f = as.factor(D)) %>%
  ggplot(aes(x = full_cb_size, y = slim_cb_critical_release_size, ymin = slim_cb_q0.05, 
             ymax = slim_cb_q0.95)) +
  geom_abline(slope = 1, intercept = 0, color = "darkgrey", linewidth = 0.75) +
  geom_point(aes(color = phat_f, shape = s_f, size = D_f, alpha = slim_bubble_div_full_bubble)) +
  geom_errorbar(aes(color = phat_f,
                    alpha = slim_bubble_div_full_bubble)) +
  labs(alpha = "slim auc / \nfull auc", 
       color = TeX(r"(\textit{$\hat{p}})"), shape = TeX(r"(\textit{$s_u$})")) +
  scale_size_discrete(name = TeX(r"(\textit{$D$})"),
                      breaks = c(1e-5, 1e-4, 1e-3, 1e-2),
             labels = c(TeX(r"($10^{-5}$)"),TeX(r"($10^{-4}$)"),TeX(r"($10^{-3}$)"), TeX(r"($10^{-2}$)")),
             range = c(min_pt_size, max_pt_size))+
  scale_alpha_continuous(name = "fraction of critical\nbubble represented",
              range = c(min.alpha,max.alpha),
              breaks = c(min.alpha, 0.5, max.alpha),
              labels = c("0.17","0.5","1.0"))+
  xlab("Critical bubble release size") +
  ylab("Individual-based model critical release size") +
  theme_minimal() +
  scale_x_log10(breaks = c(1,10,100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,175000)) + 
  scale_y_log10(breaks = c(10,100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,175000)) +
  annotation_logticks() +
  theme(panel.grid.minor = element_blank())+
  sizes.theme + guides(shape = guide_legend(override.aes = list(size = min_pt_size)),
                      color = guide_legend(override.aes = list(size = min_pt_size)),
                      alpha = guide_legend(override.aes = list(size = min_pt_size)))


panelC_nolegend = panelC_legend + theme(legend.position = "none") +
  theme(axis.title.y = element_text(margin = margin(r = 10)))

# ggsave(plot = panelC_legend + theme(axis.title.y = element_text(margin = margin(r = 10))),
#        width = 12, height = 12,
#        filename = paste0(fig.dir,"critical-bubble-in-slim/critical-bubble-slim-results-subplotC-with-legend.pdf"))
# 
# ggsave(plot = panelC_nolegend,
#        width = 12, height = 12,
#        filename = paste0(fig.dir,"critical-bubble-in-slim/cb-panelC.pdf"))


### panels for the legend ######
param.no = 34
i = which(cb.data$param_no == param.no)
param.no = cb.data$param_no[i]
D = cb.data$D[i]
s = cb.data$s[i]
phat = cb.data$phat[i]
bw = cb.data$full_cb_auc[i]

legend_in = 20
linewidth = 1.0
p34 = plot_bubble_chopped_off(phat = phat, D = D, s = s, line_thickness = linewidth, line_alpha = 1) +
  theme_minimal() +
  scale_x_continuous(limits = c(0.5 - 7.5, 0.5 + 7.5),expand = c(0,0),
                     breaks = seq(-7, 8, by = 3)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0)) +
  ylab(TeX(r"(\textit{p})")) + xlab("") + 
  sizes.theme + xlim(-2, 3) +
  theme(axis.text.x=element_blank(), 
      axis.ticks.x=element_blank(), 
      axis.text.y=element_blank(), 
      axis.ticks.y=element_blank()) +
  theme(axis.title.y = element_text(margin = margin(r = legend_in)))

# ggsave(plot = p34, width = 2.5, height = 2.5,
#        filename = paste0(fig.dir,"critical-bubble-in-slim/cb-panelC-legend-low-end.pdf"))


param.no = 23
i = which(cb.data$param_no == param.no)
param.no = cb.data$param_no[i]
D = cb.data$D[i]
s = cb.data$s[i]
phat = cb.data$phat[i]
bw = cb.data$full_cb_auc[i]

p23 = plot_bubble_chopped_off(phat = phat, D = D, s = s, line_thickness = linewidth, line_alpha = 1) +
  theme_minimal() +
  scale_x_continuous(limits = c(-1,2),expand = c(0,0),
                     breaks = seq(0.5 - 1.5, 0.5 + 1.5, by = 1)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0)) +
  xlim(-2, 3) +
  ylab(TeX(r"(\textit{p})")) + xlab(TeX(r"(\textit{x})")) + 
  sizes.theme +
  theme(axis.text.x=element_blank(), 
      axis.ticks.x=element_blank(), 
      axis.text.y=element_blank(), 
      axis.ticks.y=element_blank()) +
    theme(axis.title.y = element_text(margin = margin(r = legend_in))) 

 
# ggsave(plot = p23, width = 2.5, height = 2.5,
#        filename = paste0(fig.dir,"critical-bubble-in-slim/cb-panelC-legend-high-end.pdf"))
```


## 3. Results of a rectangle release in SLiM compared to release under the critical bubble of diffusion theory

```{r eval=F}
options(scipen=999) 

t.sizes = calculate_text_size(artboard_width = 1306.62, overleaf_scaling = 1)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))

line_thickness = 1.0


### panel A and B in same column
param = 14
i = which(under.data$param_no==param)
param.no = under.data$param_no[i]
cb_width_here = under.data$full_cb_size[i]/100000
phat = under.data$phat[i]
D = under.data$D[i]
s = under.data$s[i]

title = TeX(r"(\textit{$D = 0.0001, \hat{p} = 0.05, s_u = 0.3$})")

panelA = plot_bubble_and_square(phat = phat, D=D, s = s, square.color = "black",
                                         line_thickness = 1.3, alpha = 0.2) + 
  scale_x_continuous(limits = c(0,1), expand = c(0,0)) + 
  scale_y_continuous(limits = c(0,1), expand = c(0,0),
                     breaks = c(0.25, 0.5, 0.75, 1.0))  +
  theme_minimal()+
  xlab(TeX(r"(\textit{$x$})")) + ylab(TeX(r"(\textit{$p$})"))+
 sizes.theme 

param = 18
i = which(under.data$param_no==param)
param.no = under.data$param_no[i]
cb_width_here = under.data$full_cb_size[i]/100000
phat = under.data$phat[i]
D = under.data$D[i]
s = under.data$s[i]

title = TeX(r"(\textit{$D = 0.00001, \hat{p} = 0.4, s_u = 0.8$})")

panelB = plot_bubble_and_square(phat = phat, D = D, s = s, square.color = "black",
                                         line_thickness = 1.3, alpha = 0.2)+ 
  scale_x_continuous(limits = c(0,1), expand = c(0,0)) + 
  scale_y_continuous(limits = c(0,1), expand = c(0,0), 
                     breaks = c(0.25, 0.5, 0.75, 1.0))  +
  theme_minimal()+
  xlab(TeX(r"(\textit{$x$})")) + ylab(TeX(r"(\textit{$p$})"))+
  sizes.theme

dist = -0.15
combined = plot_grid(panelA, 
                     NULL,
                     panelB,
                     ncol = 1, nrow = 3,
                     align = "hv",
                     rel_heights = c(1, dist, 1))

# ggsave(plot = combined,
#        width = 6, height = 12,
#        filename = paste0(fig.dir,"/underdominance-rectangle/compared-to-CB-size/bubble-vs-rectangle-panelA-and-B.pdf"))
# 

# panel C
D.labels = c(TeX(r"($10^{-5}$)"),TeX(r"($10^{-4}$)"),TeX(r"($10^{-3}$)"), TeX(r"($10^{-2}$)"))

yaxis.in = 10

panelC = under.data %>% 
  mutate(s_f = as.factor(s), D_f = as.factor(D), phat_f = as.factor(phat)) %>%
  ggplot(aes(x =full_cb_size, y = slim_rectangle_critical_release_size,
             ymin = slim_rectangle_q0.05, ymax = slim_rectangle_q0.95)) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linewidth = 0.75) +
  geom_point(aes(size = D, shape = s_f, color = phat_f)) +
  geom_errorbar(aes(color = phat_f)) +
  labs(shape = TeX(r"(\textit{$s_u$})"), color = TeX(r"(\textit{$\hat{p}$})"), size = TeX(r"(\textit{$D$})")) +
  xlab("Critical bubble release size") +
  ylab("Individual-based model critical release size\n(from rectangle)") +
  scale_x_log10(breaks = c(1, 10, 100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,175000)) + 
  scale_y_log10(breaks = c(10, 100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,175000)) +
  theme_minimal() +
  scale_size(name = TeX(r"(\textit{$D$})"), range = c(min_pt_size, max_pt_size),
             breaks = c(1e-5, 1e-4, 1e-3, 1e-2), 
             labels = D.labels) +
  annotation_logticks() +
  theme(panel.grid.minor = element_blank()) +
  sizes.theme + guides(shape = guide_legend(override.aes = list(size = min_pt_size))) +
  theme(axis.title.y = element_text(margin = margin(r = yaxis.in)))


# ggsave(plot = panelC,
#        width = 12, height = 12,
#        filename = paste0(fig.dir,"/underdominance-rectangle/compared-to-CB-size/bubble-vs-rectangle-panelC-legend.pdf"))


panelC.no.legend = panelC +
  theme(legend.position = "none")

# ggsave(plot = panelC.no.legend,
#        width = 12, height = 12,
#        filename = paste0(fig.dir,"/underdominance-rectangle/compared-to-CB-size/bubble-vs-rectangle-panelC-nolegend.pdf"))
```


## 4. Results of a rectangle release in SLiM compared to *deSolve* predictions

```{r eval=F}
t.sizes = calculate_text_size(artboard_width = 1306.62, overleaf_scaling = 1)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))


### Panel A: Evolution of wave for paramset13 (D = 1e-5, s = 0.3, phat = 0.2, predicted width = 0.0075)
# paramset13 (D = 1e-5, s = 0.3, phat = 0.2, predicted width = 0.0075)

d30 = under.data %>% filter(param_no==13)
phat = d30$phat; D = d30$D; s = d30$s; pred.width = d30$deSolve_rectangle_critical_size/100000

times2 = c(0, 20, 60, 100) 

line.thickness = 1.0
plot.list = plot_deSolve_wave_separate_times(phat = phat, s = s, D = D,
                                             release_width = pred.width, 
                                             times = times2, text_size = text_size,
                                             critical_bubble_color = "red",
                                             shared_x_axis = T, 
                                             critical_bubble_linetype = "solid",
                                             line.thickness = linewidth, titles = F)

### combine plots
plot1 = plot.list$`t = 0` + theme(axis.text.x = element_blank(), axis.title.x = element_blank())
plot2 = plot.list$`t = 20` + theme(axis.text.x = element_blank(), axis.title.x = element_blank())
plot3 = plot.list$`t = 60` + theme(axis.text.x = element_blank(), axis.title.x = element_blank())
plot4 = plot.list$`t = 100`

dist = -0.448
combined = plot_grid(plot1,
                     NULL,
                     plot2,
                     NULL,
                     plot3, 
                     NULL,
                     plot4, 
                     nrow = 7, ncol = 1, align = "hv",
                     rel_heights = c(1, dist, 1, dist, 1, dist, 1))

# ggsave(plot = combined,
#        filename = paste0(fig.dir,"/underdominance-rectangle/deSolve/wave-evolution.pdf"),
#        width = 6, height = 12)

### Panel C: deSolve predicted release size versus critical release size under SLiM rectangle

options(scipen=999)
yaxis.in = 10

panelC = under.data %>% 
  mutate(s_f = as.factor(s), D_f = as.factor(D), phat_f = as.factor(phat)) %>%
  ggplot(aes(x = deSolve_rectangle_critical_size, y = slim_rectangle_critical_release_size,
             ymin = slim_rectangle_q0.05, ymax = slim_rectangle_q0.95)) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linewidth = 0.75) +
  geom_point(aes(size = D, shape = s_f, color = phat_f)) +
  geom_errorbar(aes(color = phat_f)) +
  labs(shape =  TeX(r"(\textit{$s_u$})"), color = TeX(r"(\textit{$\hat{p}$})"), size =  TeX(r"(\textit{$D$})")) +
  xlab("Reaction-diffusion model critical release size") +
  ylab("Individual-based model critical release size") +
  scale_x_log10(breaks = c(1, 10, 100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,175000)) + 
  scale_y_log10(breaks = c(10, 100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,175000)) +
  theme_minimal() +
  scale_size(name = TeX(r"(\textit{$D$})"), range = c(min_pt_size, max_pt_size),
             breaks = c(1e-5, 1e-4, 1e-3, 1e-2), 
             labels = D.labels) +
  annotation_logticks() +
  theme(panel.grid.minor = element_blank()) +
  sizes.theme + guides(shape = guide_legend(override.aes = list(size = min_pt_size))) +
  theme(axis.title.y = element_text(margin = margin(r = yaxis.in)))

panelC_no_legend = panelC +
  theme(legend.position = "none")

# ggsave(plot = panelC,
#        width = 12, height = 12,
#        filename = paste0(fig.dir,"/underdominance-rectangle/deSolve/deSolve-panelC-with-legend.pdf"))
# 
# ggsave(plot = panelC_no_legend,
#        width = 12, height = 12,
#        filename = paste0(fig.dir,"/underdominance-rectangle/deSolve/deSolve-panelC-no-legend.pdf"))
```

## 5. TADE diagram

* This was made in PowerPoint

## 6. TADE modification 2D
```{r eval=F}
t.sizes = calculate_text_size(artboard_width = 900, overleaf_scaling = 0.6)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))

### 2D results (main)
options(scipen=999)
plt2D = tade.mod2D.data %>% 
  mutate(s_phat = paste0(s, ", ", round(phat, 3))) %>%
  mutate(D_f = as.factor(D)) %>%
  ggplot(aes(x = deSolve_circle_critical_size, y = slim_circle_critical_release_size, 
             ymin = slim_circle_q0.05, ymax = slim_circle_q0.95)) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linewidth = 0.75) +
  geom_point(aes(color = s_phat, size = D)) +
  geom_errorbar(aes(color = s_phat)) +
  labs(color = TeX(r"(\textit{$s_g, \hat{p}$})"), size = TeX(r"(\textit{$D$})")) +
  xlab("Reaction-diffusion model critical release size") +
  ylab("Individual-based model critical release size") + 
  scale_x_log10(breaks = c(1, 10, 100, 1000, 10000, 100000),  expand = expansion(mult = 0.05), limits = c(1,100000)) + 
  scale_y_log10(breaks = c(10, 100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,100000)) +
  theme_minimal() +
  scale_size(name = TeX(r"(\textit{$D$})"), range = c(min_pt_size, max_pt_size),
             breaks = c(1e-5, 1e-4, 1e-3, 1e-2), 
             labels = D.labels) +
  annotation_logticks() +
  theme(panel.grid.minor = element_blank()) +
  sizes.theme +  guides(shape = guide_legend(override.aes = list(size = min_pt_size))) +
  theme(axis.title.y = element_text(margin = margin(r = yaxis.in)))

# ggsave(plot = plt2D,
#        filename = paste0(fig.dir,"tade-modification/twoD/tade-mod-2D-legend.pdf"),
#        width = 12, height = 12)
# 
# ggsave(plot = plt2D+ theme(legend.position = "none"),
#        filename = paste0(fig.dir,"tade-modification/twoD/tade-mod-2D-no-legend.pdf"),
#        width = 12, height = 12)
```


## S1. How well the critical bubble is approximated in SLiM

```{r eval=F}
t.sizes = calculate_text_size(artboard_width = 458, overleaf_scaling = 0.3)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))

linewidth = 1.2

### Panel A: paramset 11 - D = 1e-5, s = 0.05, phat = 0.05 
param.no = 11
i = which(cb.data$param_no == param.no)
D = cb.data$D[i]
s = cb.data$s[i]
phat = cb.data$phat[i]
bw = cb.data$full_cb_auc[i]

title = TeX(r"(\textit{$D = 0.00001, \hat{p} = 0.05, s_u = 0.05$})")

panelA = plot_bubble_chopped_off(phat = phat, D = D, s = s, line_thickness = linewidth, line_alpha = 1) +
  scale_x_continuous(limits = c(-1, 2), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0)) +
  theme_minimal()+ ylab(TeX(r"(\textit{$p$})"))+ xlab(TeX(r"(\textit{$x$})")) +
  sizes.theme

### Panel B - paramset41: D = 0.001, phat = 0.2, s = 0.05 
param.no = 41
i = which(cb.data$param_no == param.no)
param.no = cb.data$param_no[i]
D = cb.data$D[i]
s = cb.data$s[i]
phat = cb.data$phat[i]
bw = cb.data$full_cb_auc[i]

title = TeX(r"(\textit{$D = 0.001, \hat{p} = 0.2, s_u = 0.05$})")

panelB = plot_bubble_chopped_off(phat = phat, D = D, s = s, line_thickness = linewidth, line_alpha = 1) +
  theme_minimal() +
  scale_x_continuous(limits = c(-1, 2), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0)) +
  xlab(TeX(r"(\textit{$x$})")) + ylab(TeX(r"(\textit{$p$})"))+
  sizes.theme

dist = -0.25
together = plot_grid(panelA,
                     NULL,
                     panelB,
                     ncol = 1, nrow = 3, align = "hv",
                     rel_heights = c(1,dist,1))


# ggsave(plot = together, width = 6, height = 12,
#        filename = paste0(fig.dir,"critical-bubble-in-slim/together-approx-of-critical-bubble.pdf"))
```

## S2. The critical local frequency for different release widths for paramset8 (via deSolve algorithm)
```{r eval=F}
t.sizes = calculate_text_size(artboard_width = 900, overleaf_scaling = 0.5)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))

data_D1e_5 = read.xlsx(here("data/summaries/critical-p0-underdominance/critical-p0-D1e-5.xlsx")) %>% as_tibble()
data_D1e_4 = read.xlsx(here("data/summaries/critical-p0-underdominance/critical-p0-D1e-4.xlsx")) %>% as_tibble()
data_D1e_3 = read.xlsx(here("data/summaries/critical-p0-underdominance/critical-p0-D1e-3.xlsx")) %>% as_tibble()

title = TeX(r"(\textit{$\hat{p} = 0.2, s_u = 0.8$})")

palette = carto_pal(4, "Emrld")[2:4]

this.D = c(TeX(r"($10^{-5}$)"),TeX(r"($10^{-4}$)"),TeX(r"($10^{-3}$)"))

linewidth = 1.2
p.comb = ggplot(data_D1e_5, aes(x = release_width, y = critical_p0)) +
  geom_line(aes(color = "a"), size = linewidth) +
  ylab(TeX(r"(Critical \textit{$p_0$})")) +
  xlab("Release width") +
  geom_hline(yintercept = 0.2, color = "darkgrey",  size = linewidth) +
  theme_minimal()+
  sizes.theme +
  geom_line(data = data_D1e_4, aes(x = release_width, y = critical_p0, color = "b"), 
            size = linewidth) +
  geom_line(data = data_D1e_3, aes(x = release_width, y = critical_p0, color = "c"), 
            size = linewidth) +
  scale_colour_manual(name = TeX(r"(\textit{$D$})"), values = c("a" = palette[1], "b" = palette[2], "c" = palette[3]), 
                      labels =this.D) +
  scale_y_continuous(expand = c(0,0), limits = c(0,1), breaks = c(0.2, 0.4, 0.6, 0.8,1.0)) 

p.comb.noleg = p.comb + theme(legend.position = "none")

# ggsave(plot = p.comb, 
#        filename = paste0(fig.dir,"underdominance-rectangle/deSolve/critical-p0-withlegend.pdf"),
#        height = 12, width = 12)
# 
# ggsave(plot = p.comb.noleg, 
#        filename = paste0(fig.dir,"underdominance-rectangle/deSolve/critical-p0-nolegend.pdf"),
#        height = 12, width = 12)

# scratch
data5 = data_D1e_5 %>% mutate(minimum_area = critical_p0*release_width)
data4 = data_D1e_4  %>% mutate(minimum_area = critical_p0*release_width)
data3 = data_D1e_3  %>% mutate(minimum_area = critical_p0*release_width)
```


## S3. 1D homing gene drive results

```{r eval=F}
t.sizes = calculate_text_size(artboard_width = 900, overleaf_scaling = 0.6)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))

yaxis.in = 12

options(scipen=999)
gd.1D = gd.data %>% 
  filter(s!=0.05)%>%
  mutate(s_f = as.factor(s), D_f = as.factor(D), pg_hat_f = as.factor(phat)) %>%
  ggplot(aes(x = deSolve_rectangle_critical_size, y = slim_rectangle_critical_release_size, 
             ymin = slim_rectangle_q0.05, ymax = slim_rectangle_q0.95)) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linewidth = 0.75) +
  geom_point(aes(color = pg_hat_f, shape = s_f, size = D)) +
  geom_errorbar(aes(color = pg_hat_f)) +
  labs(shape = TeX(r"(\textit{$s_g$})"), color = TeX(r"(\textit{$\hat{p}$})"), size = TeX(r"(\textit{$D$})")) +
  xlab("Reaction-diffusion model critical release size") + 
  ylab("Individual-based model critical release size") + 
  scale_x_log10(breaks = c(1, 10, 100, 1000, 10000, 100000),  expand = expansion(mult = 0.05), limits = c(1,100000)) + 
  scale_y_log10(breaks = c(10, 100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,100000)) +
  theme_minimal() +
  scale_size(name = TeX(r"(\textit{$D$})"), range = c(min_pt_size, max_pt_size),
             breaks = c(1e-5, 1e-4, 1e-3, 1e-2), 
             labels = D.labels) +
  annotation_logticks() +
  theme(panel.grid.minor = element_blank()) +
  sizes.theme +  guides(shape = guide_legend(override.aes = list(size = min_pt_size))) +
  theme(axis.title.y = element_text(margin = margin(r = yaxis.in)))

gd.1D.nolegend = gd.1D +
  theme(legend.position = "none")

# ggsave(plot = gd.1D,
#        filename = paste0(fig.dir,"gene-drive-1D/gene-drive-1D-with-legend.pdf"),
#        width = 12, height = 12)
# 
# ggsave(plot = gd.1D.nolegend,
#        filename = paste0(fig.dir,"gene-drive-1D/gene-drive-1D-with-nolegend.pdf"),
#        width = 12, height = 12)
```

## S4. 1D TADE modification
```{r eval=F}
t.sizes = calculate_text_size(artboard_width = 900, overleaf_scaling = 0.6)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))

print(paste0("subplot label text size = ", t.sizes$subplot_text_size))

### 1D results (supplemental)
yaxis.in = 12
options(scipen=999)
plt = tade.mod1D.data %>%
  mutate(s_phat = paste0(s, ", ", round(phat, 3))) %>%
  mutate(D_f = as.factor(D)) %>%
  ggplot(aes(x = deSolve_rectangle_critical_size, y = slim_rectangle_critical_release_size, 
             ymin = slim_rectangle_q0.05, ymax = slim_rectangle_q0.95)) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linewidth = 0.75) +
  geom_point(aes(color = s_phat, size = D)) +
  geom_errorbar(aes(color = s_phat)) +
  labs(color = TeX(r"(\textit{$s_g, \hat{p}$})"), size = TeX(r"(\textit{$D$})")) +
  xlab("Reaction-diffusion model critical release size") +
  ylab("Individual-based model critical release size") +
  scale_x_log10(breaks = c(1, 10, 100, 1000, 10000, 100000),  expand = expansion(mult = 0.05), limits = c(1,100000)) + 
  scale_y_log10(breaks = c(10, 100, 1000, 10000, 100000), expand = c(0,0), limits = c(1,100000)) +
  theme_minimal() +
  scale_size(name = TeX(r"(\textit{$D$})"), range = c(min_pt_size, max_pt_size),
             breaks = c(1e-5, 1e-4, 1e-3, 1e-2), 
             labels = D.labels) +
  annotation_logticks() +
  theme(panel.grid.minor = element_blank()) +
  sizes.theme +  guides(shape = guide_legend(override.aes = list(size = min_pt_size))) +
  theme(axis.title.y = element_text(margin = margin(r = yaxis.in)))


# ggsave(plot = plt,
#        filename = paste0(fig.dir,"tade-modification/oneD/tade-mod-1D-legend.pdf"),
#          width = 12, height = 12)
# 
# ggsave(plot = plt+ theme(legend.position = "none") ,
#        filename = paste0(fig.dir,"tade-modification/oneD/tade-mod-1D-no-legend.pdf"),
#        width = 12, height = 12)
```


## S5. TADE modification 1D versus 2D SLiM results
```{r eval=F}
t.sizes = calculate_text_size(artboard_width = 900, overleaf_scaling = 0.6)
text_size = t.sizes$text_size
sizes.theme =  theme(
    axis.title.x = element_text(size = text_size,
                                margin = margin(t = 10, r = 0, b = 0, l = 0)), 
    axis.title.y = element_text(size = text_size,
                                margin = margin(t = 0, r = 10, b = 0, l = 0)),  
    axis.text.x = element_text(size = text_size),     
    axis.text.y = element_text(size = text_size),
    legend.text = element_text(size = text_size),
    legend.title = element_text(size = text_size),
    plot.margin = unit(c(1,1,1,1), "cm"),
    panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 1),
    aspect.ratio = 1,
    plot.title = element_text(hjust=0.5,size = text_size))


tade.slim = bind_cols(tade.mod1D.data, tade.mod2D.data %>% 
                        dplyr::select(slim_circle_q0.05, slim_circle_q0.95, slim_circle_critical_release_size))

tade2 = tade.slim %>%
  mutate(width = slim_rectangle_critical_release_size/100000,
         width.low = slim_rectangle_q0.05/100000,
         width.high = slim_rectangle_q0.95/100000,
         diameter = 2*sqrt(slim_circle_critical_release_size/(100000*pi)),
         diameter.low = 2*sqrt(slim_circle_q0.05/(100000*pi)),
         diameter.high = 2*sqrt(slim_circle_q0.95/(100000*pi)))

# x-axis is 1D, y-axis is 2D
ex = 0.1
yaxis.in = 12
options(scipen=999)
plt = tade2 %>% 
  mutate(s_phat = paste0(s, ", ", round(phat, 3))) %>%
  mutate(D_f = as.factor(D)) %>%
  ggplot(aes(x = width, xmin = width.low, xmax = width.high,
             y = diameter, ymin = diameter.low, ymax = diameter.high)) +
  geom_abline(slope = 1, intercept = 0, color = "grey", linewidth = 0.75) +
  geom_point(aes(color = s_phat, size = D)) +
  geom_errorbar(aes(color = s_phat)) +
  geom_errorbarh(aes(color = s_phat)) +
  labs(color = TeX(r"(\textit{$s_g, \hat{p}$})"), size = TeX(r"(\textit{$D$})")) +
  xlab("Individual-based model critical width") + 
  ylab("Individual-based model critical diameter") + 
  theme_minimal() +
  annotation_logticks(base = 10) +
  scale_x_log10(limits = c(0.00001,1), breaks = c(0.00001, 0.0001, 0.001, 0.01, 0.1, 1),
                labels = c("0.00001", "0.0001", "0.001", "0.01", "0.1", "1"), expand = expansion(mult = ex)) +
  scale_y_log10(limits = c(0.00001,1), breaks = c(0.00001, 0.0001, 0.001, 0.01, 0.1, 1),
                labels = c("0.00001", "0.0001", "0.001", "0.01", "0.1", "1"), expand = expansion(mult = ex)) +
  scale_size(name = TeX(r"(\textit{$D$})"), range = c(min_pt_size, max_pt_size),
             breaks = c(1e-5, 1e-4, 1e-3, 1e-2), 
             labels = D.labels) +
  theme(panel.grid.minor = element_blank()) +
  sizes.theme +  guides(shape = guide_legend(override.aes = list(size = min_pt_size))) +
  theme(axis.title.y = element_text(margin = margin(r = yaxis.in)))


# ggsave(plot = plt,
#        filename = paste0(fig.dir,"tade-modification/tade-mod-1D-vs-2D-legend.pdf"),
#          width = 12, height = 12)
# 
# ggsave(plot = plt+ theme(legend.position = "none") ,
#        filename = paste0(fig.dir,"tade-modification/tade-mod-1D-vs-2D-no-legend.pdf"),
#        width = 12, height = 12)
```