library(tidyverse)
library(samplr)
library(abc)
select <- dplyr::select # abc package loads MASS which has a function called select
# but so does dplyr! So if you load tidyverse first then select will be MASS::select 
# (because you loaded it last). You can solve this by loading tidyverse last or by
# doing what we're doing in line 4. Just thought to flag it here as it's an error
# we've encountered often! 


# Setup
sim_summaries <- readRDS("simulations/sim_summaries.rds")
obs_summaries <- readRDS("data/obs_summaries.rds")

obs_summaries[1:2,]
sim_summaries[1:2,]

# Practice: Find out which model best fits you! ---------------------------

# You can comment some of these out and see if it makes a difference
summaries_to_use <- c(
  'R', 
  'A', 
  'TPF', 
  'D', 
  'S', 
  'RNG',
  'gzip'
)



fit <- abc(
  target= obs_summaries[obs_summaries$pseudonym == "lks", summaries_to_use],
  param = sim_summaries$models,
  sumstat = sim_summaries[, summaries_to_use],
  tol = , # fill this field (and see how different values make a difference)
  method = "rejection"
)

table(fit$unadj.values) / sum(table(fit$unadj.values))



# Practice: Do this for all participants and combine the likelihoods --------
# Loop iterates over all participants...
df <- data.frame()
for (i in unique(obs_summaries$id_number)){
  fit <- abc(
    target= obs_summaries[obs_summaries$id_number == i, 3:9],
    param = sim_summaries$models,
    sumstat = sim_summaries[, 7:13],
    tol = .05, method = "rejection"
  )
  result <- cbind(
    data.frame(
      id_number=i,
      pseudonym=obs_summaries$pseudonym[obs_summaries$id_number == i]
    ),
    data.frame(table(fit$unadj.values) / sum(table(fit$unadj.values)))
  )
  df <- rbind(df, result)
}

## We can plot the results... 
df %>%
  ggplot(aes(pseudonym, Freq, fill=Var1)) +
  geom_col(position = "stack", color="black") +
  scale_fill_brewer(type = "qual") +
  guides(fill=guide_legend(title="Model")) +
  theme_minimal(21) +
  xlab("Pseudonym") +
  ylab("p") +
  theme(legend.position = "bottom") +
  theme(
    panel.background = element_rect(fill="#eeeeee", colour="#eeeeee"),
    plot.background = element_rect(fill="#eeeeee", colour="#eeeeee"),
    axis.text.x = element_text(angle=90, hjust = 1)
  )


## Combine these probabilities into evidence for models at the sample level

# Your code here 
# 



# Parameter Estimation ----------------------------------------------------
# change this to you:
pseudonym <- "lks"
winning_model <- "MCREC"

# then run this:
params <- c()
if (winning_model %in% c("MH", "MC3", "RECMH", "RECMC3")) {
  params <- c(params, "sigma")
} else {
  params <- c(params, "L")
}

if (winning_model %in% c("RECMH", "RECMC3", "RECHMC", "MCREC")){
  params <- c(params, "alpha")
}

# then this:
fit <- abc(
  target= obs_summaries[
    obs_summaries$pseudonym == pseudonym, 3:9 
  ],
  param = sim_summaries[
    sim_summaries$models == winning_model, 
    params 
  ],
  sumstat = sim_summaries[
    sim_summaries$models == winning_model, 7:13 
  ],
  tol = .1,
  method = "rejection"
)

# plot
if (length(params) == 1){
  as_tibble(fit$unadj.values) %>% 
    magrittr::set_colnames(params) %>% 
    pivot_longer(all_of(params)) %>% 
    ggplot(aes(value)) + 
    geom_histogram() + 
    xlab(params)
} else {
  as_tibble(fit$unadj.values) %>% 
    magrittr::set_colnames(c("A", "B")) %>% 
    ggplot(aes(A, B)) + 
    geom_point(shape=21, fill="grey", alpha=.7) + 
    xlab(params[1]) + ylab(params[2])
}
