retardation_test_conditioned_inhibition.Rmd

---
title: "Retardation test for conditioned inhibition"
author: "Deniz Tuzsus"
date: "12 November 2019"
output: html_document
---

```{r}
# setup constant variables over all simulations
n_features <- 120
cue_features <- 1:100
a <- b  <- c <- d <- context <- outcome <- rep(0, n_features)

# setup stimuli
a[1:20] <- b [21:40] <- c[41:60] <- d[61:80] <- 1
context[81:100] <- 1
outcome[101:120] <- 1

# number of replications
n_replications <- 100

# learning rates
p_encode <- c(0.33, 0.67, 1)

# specify models
model = "Minerva AL"
```

```{r}
# specify events
training_ax_event <- a + context + outcome
training_ab_event <- a + b + context
test_bcx_event <- b + c + context + outcome
test_cx_event <- c + context + outcome

probe_a <- a + context
probe_ab <- a + b + context
probe_bc <- b + c + context
probe_c <- c + context

# setup trial number
n_trials <- 200

# prepare results object as 3D array[p_encode, n_trials, n_replications]
sim_results_exp <- array(NA, c(length(p_encode), n_trials, n_replications))
sim_results_control <- array(NA, c(length(p_encode), n_trials, n_replications))
# specify object for normalized echo for every trial as 3D array[p_encode, n_trials, n_replications]
# with cells containing normalized echo vector 
normalized_echos_exp <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))
normalized_echos_control <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))
```

```{r}
#------------------------Experimental Condition--------------------------------#

for (m in 1:length(model)){
  for (r in 1:n_replications) {
    for(i in 1:length(p_encode)) {
      # Memory is empty on first trial
      normalized_echo <- probe_memory(probe_ab, NULL, cue_features, model = model[m])
      memory <- learn(
        normalized_echo
        , training_ab_event
        , p_encode[i]
        , NULL
        , model = model[m]
      )
          normalized_echos_exp[[i,1,r]] <- normalized_echo
      
       # A -> X trials on every even number
      for(j in 2:100) {
        if (j %% 2 == 0){
          normalized_echo <- probe_memory(probe_a, memory, cue_features, model = model[m])
          memory <- learn(
            normalized_echo
            , training_ax_event
            , p_encode[i]
            , memory
            , model = model[m]
          )
          normalized_echos_exp[[i,j,r]] <- normalized_echo
        }
        else {
          # AB trials on every odd number
          normalized_echo <- probe_memory(probe_ab, memory, cue_features, model = model[m])
          memory <- learn(
            normalized_echo
            , training_ab_event
            , p_encode[i]
            , memory
            , model = model[m]
          )
          normalized_echos_exp[[i,j,r]] <- normalized_echo
        }
      }
      # BC -> X trials
      for(j in 101:200) {
        normalized_echo <- probe_memory(probe_bc, memory, cue_features, model = model[m])
        expectancy <- expect_event(outcome, normalized_echo)
        memory <- learn(
          normalized_echo
          , test_bcx_event 
          , p_encode[i]
          , memory
          , model = model[m]
        )
        sim_results_exp[[i,j,r]] <- expectancy
        normalized_echos_exp[[i,j,r]] <- normalized_echo
      }
    }
  }
}

#-----calculate mean and standard error of results objects-----#
sim_results_exp <- sim_results_exp[,101:200,]
divide_by_n_rep <- function(x, n_rep = n_replications){x/n_rep}
# Sum of expectancies for each trial over all replications
mean_sim_results_exp <- apply(sim_results_exp, c(1,2), FUN = sum)
# Mean of expectancies for each trial over all replications
mean_sim_results_exp <- apply(mean_sim_results_exp, c(1,2), FUN = divide_by_n_rep)
# critical trial until x|cue >= 0.95 in each replication for every learning rate
crit_trials <- apply(sim_results_exp, c(1,3), FUN = find_crit_trial)
# se of critical trials until x|cue >= 0.95 in each replication for every learning rate
std_error <- apply(crit_trials, 1, FUN = std_err)
```

```{r}
res_mat <- matrix(NA, nrow = 2, ncol = 3)

#compute means over replications
x <- rowMeans(crit_trials)
#compute standard errors over replications
y <- std_error

for (i in 1:length(p_encode)){
  res_mat[1,i] <- table_fun(x[i], y[i])
}

```


```{r, cache=TRUE}
# -------------------- Diagnostic Plots --------------------------------------#

# par(mfrow = c(2,2))
# 
# # memory encoding of the outcome over the trials
# plot(
#   1:200
#   , memory[, 101]
#   , type = "l"
#   , col = scales::alpha("black", 0.3)
#   , ylim = c(-2, 2)
#   , xlab = "Trial"
#   , ylab = "Feature encoding"
#   , main = "Features of outcome X"
#   , las = 1
# )
# for(i in 102:120) {
#   lines(
#     1:200
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
# 
# # memory encoding of cue A over the trials
# plot(
#   1:200
#   , memory[, 1]
#   , type = "l"
#   , col = scales::alpha("black", 0.3)
#   , ylim = c(-2, 2)
#   , xlab = "Trial"
#   , ylab = "Feature encoding"
#   , main = "Features of Cue A"
#   , las = 1
# )
# for(i in 2:20) {
#   lines(
#     1:200
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
# 
# # memory encoding of cue B over the trials
# plot(
#   1:n_trials
#   , memory[, 21]
#   , type = "l"
#   , col = scales::alpha("black", 0.3)
#   , ylim = c(-2, 2)
#   , xlab = "Trial"
#   , ylab = "Feature encoding"
#   , main = "Features of Cue B"
#   , las = 1
# )
# for(i in 22:40) {
#   lines(
#     1:n_trials
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
# 
# # memory encoding of cue C over the trials
# plot(
#   1:n_trials
#   , memory[, 41]
#   , type = "l"
#   , col = scales::alpha("black", 0.3)
#   , ylim = c(-2, 2)
#   , xlab = "Trial"
#   , ylab = "Feature encoding"
#   , main = "Features of Cue C"
#   , las = 1
# )
# for(i in 42:60) {
#   lines(
#     1:n_trials
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
```

```{r, cache=TRUE}
# normalized echos over the trials

# echos <- vector()
# for (i in 1:n_trials){
#   echos <- c(echos,normalized_echos_exp[[3,i,1]])
# }
# 
# data <- data.frame(Stimuli = c(1:n_features), Echo = echos, frame = rep(1:n_trials, each = n_features))
# 
# 
# # Make an animated boxplot
# p <- ggplot(data, aes(x=Stimuli, y=Echo, fill=Stimuli)) + 
#   geom_bar(stat='identity') +
#   theme_bw() +
#   # gganimate specific bits:
#   transition_states(
#     states = frame,
#     transition_length = 1,
#     state_length = 1, 
#     wrap = TRUE
#   ) 
# 
# p <- p + ggtitle('Now showing Trial:{closest_state}')
# animate(p, nframes = 2*n_trials)
```

```{r}
#------------------------Control Condition--------------------------------#

for (m in 1:length(model)){
  for (r in 1:n_replications) {
    for(i in 1:length(p_encode)) {
      # Memory is empty on first trial
      normalized_echo <- probe_memory(probe_ab, NULL, cue_features, model = model[m])
      memory <- learn(
        normalized_echo
        , training_ab_event
        , p_encode[i]
        , NULL
        , model = model[m]
      )
          normalized_echos_control[[i,1,r]] <- normalized_echo
      
       # A -> X trials on every even number
      for(j in 2:100) {
        if (j %% 2 == 0){
          normalized_echo <- probe_memory(probe_a, memory, cue_features, model = model[m])
          memory <- learn(
            normalized_echo
            , training_ax_event
            , p_encode[i]
            , memory
            , model = model[m]
          )
          normalized_echos_control[[i,j,r]] <- normalized_echo
        }
        else {
          # AB trials on every odd number
          normalized_echo <- probe_memory(probe_ab, memory, cue_features, model = model[m])
          memory <- learn(
            normalized_echo
            , training_ab_event
            , p_encode[i]
            , memory
            , model = model[m]
          )
          normalized_echos_control[[i,j,r]] <- normalized_echo
        }
      }
      # C -> X trials
      for(j in 101:200) {
        normalized_echo <- probe_memory(probe_c, memory, cue_features, model = model[m])
        expectancy <- expect_event(outcome, normalized_echo)
        memory <- learn(
          normalized_echo
          , test_cx_event 
          , p_encode[i]
          , memory
          , model = model[m]
        )
        sim_results_control[[i,j,r]] <- expectancy
        normalized_echos_control[[i,j,r]] <- normalized_echo
      }
    }
  }
}

#-----calculate mean and standard error of results objects-----#
sim_results_control <- sim_results_control[,101:200,]
divide_by_n_rep <- function(x, n_rep = n_replications){x/n_rep}
# Sum of expectancies for each trial over all replications
mean_sim_results_control <- apply(sim_results_control, c(1,2), FUN = sum)
# Mean of expectancies for each trial over all replications
mean_sim_results_control <- apply(mean_sim_results_control, c(1,2), FUN = divide_by_n_rep)
# critical trial until x|cue >= 0.95 in each replication for every learning rate
crit_trials <- apply(sim_results_control, c(1,3), FUN = find_crit_trial)
# se of critical trials until x|cue >= 0.95 in each replication for every learning rate
std_error <- apply(crit_trials, 1, FUN = std_err)
```

```{r}

#compute means over replications
x <- rowMeans(crit_trials)
#compute standard errors over replications
y <- std_error

for (i in 1:length(p_encode)){
  res_mat[2,i] <- table_fun(x[i], y[i])
}

knitr::kable(cbind(Condition = c("Experimental", "Control"), res_mat), col.names = c("Condition", round(p_encode, 2)))
```


```{r, cache=TRUE}
# -------------------- Diagnostic Plots --------------------------------------#

# par(mfrow = c(2,2))
# 
# # memory encoding of the outcome over the trials
# plot(
#   1:200
#   , memory[, 101]
#   , type = "l"
#   , col = scales::alpha("black", 0.3)
#   , ylim = c(-2, 2)
#   , xlab = "Trial"
#   , ylab = "Feature encoding"
#   , main = "Features of outcome X"
#   , las = 1
# )
# for(i in 102:120) {
#   lines(
#     1:200
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
# 
# # memory encoding of cue A over the trials
# plot(
#   1:200
#   , memory[, 1]
#   , type = "l"
#   , col = scales::alpha("black", 0.3)
#   , ylim = c(-2, 2)
#   , xlab = "Trial"
#   , ylab = "Feature encoding"
#   , main = "Features of Cue A"
#   , las = 1
# )
# for(i in 2:20) {
#   lines(
#     1:200
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
# 
# # memory encoding of cue B over the trials
# plot(
#   1:n_trials
#   , memory[, 21]
#   , type = "l"
#   , col = scales::alpha("black", 0.3)
#   , ylim = c(-2, 2)
#   , xlab = "Trial"
#   , ylab = "Feature encoding"
#   , main = "Features of Cue B"
#   , las = 1
# )
# for(i in 22:40) {
#   lines(
#     1:n_trials
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
# 
# # memory encoding of cue C over the trials
# plot(
#   1:n_trials
#   , memory[, 41]
#   , type = "l"
#   , col = scales::alpha("black", 0.3)
#   , ylim = c(-2, 2)
#   , xlab = "Trial"
#   , ylab = "Feature encoding"
#   , main = "Features of Cue C"
#   , las = 1
# )
# for(i in 42:60) {
#   lines(
#     1:n_trials
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
```

```{r, cache=TRUE}
# normalized echos over the trials

# echos <- vector()
# for (i in 1:n_trials){
#   echos <- c(echos,normalized_echos_control[[3,i,1]])
# }
# 
# data <- data.frame(Stimuli = c(1:n_features), Echo = echos, frame = rep(1:n_trials, each = n_features))
# 
# 
# # Make an animated boxplot
# p <- ggplot(data, aes(x=Stimuli, y=Echo, fill=Stimuli)) + 
#   geom_bar(stat='identity') +
#   theme_bw() +
#   # gganimate specific bits:
#   transition_states(
#     states = frame,
#     transition_length = 1,
#     state_length = 1, 
#     wrap = TRUE
#  ) 

# p <- p + ggtitle('Now showing Trial:{closest_state}')
# animate(p, nframes = 2*n_trials)
```