transitive_inference.Rmd

---
title: "Transitive Inference"
author: "Tamas Foldes"
output: html_document
---


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

context1 <- context2 <- rep(0, n_features)

# setup stimuli
a[1:20] <- b [21:40] <- c[41:60] <- d[61:80] <- e[81:100] <- f[101:120] <- 1
context[121:160] <- sample(c(0,1), 40, replace = TRUE)
outcome[161:n_features] <- 0

# number of replications
n_replications <- 10

# learning rates
p_encode <- c(0.2, 0.4, 0.6, 0.8, 1)

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

```{r}
# specify events
training_a_event <- a + context
training_b_event <- b + context
training_c_event <- c + context
training_d_event <- d + context
training_e_event <- e + context
training_f_event <- f + context

training_ab_event <- a+b + context
training_bc_event <- b+c + context
training_cd_event <- c+d + context
training_de_event <- d+e + context
training_ef_event <- e+f + context

test_bc_event <- b + c + context 
recall_bc_event <- b + c + context 

test_bd_event <- b + d + context 
recall_bcd_event <- b + c + d + context 

test_be_event <- b + e + context 
recall_bcde_event <- b + c + d + e + context 


# specify probe
probe_a <- a + context
probe_b <- b + context
probe_c <- c + context
probe_d <- d + context
probe_e <- e + context
probe_f <- f + context

stim_trials <- 50
# setup trial number
n_trials <- 5*stim_trials*3

#training array
training_array <- rep(list(training_ab_event, training_bc_event, training_cd_event, training_de_event, training_ef_event), trials)
library('plot.matrix')
par(mar=c(5.1, 4.1, 4.1, 4.1)) # adapt margins
#plot(matrix(unlist(training_array), ncol=240, byrow = TRUE))
#plot(matrix(unlist(all_pair_testing_array_context2), ncol=240, byrow = TRUE))
#random shuffle
training_array <- sample(training_array)
#plot(matrix(unlist(training_array), ncol=160))

# prepare results object as 3D array[p_encode, n_trials, n_replications]
sim_results_0_pretrials <- array(NA, c(length(p_encode), n_trials, n_replications))
sim_results_10_pretrials <- array(NA, c(length(p_encode), n_trials, n_replications))
sim_results_30_pretrials <- array(NA, c(length(p_encode), n_trials, n_replications))
sim_results_100_pretrials <- array(NA, c(length(p_encode), n_trials, n_replications))

sim_results_0_adjacent <- array(NA, c(length(p_encode), n_trials, n_replications))
sim_results_0_close <- array(NA, c(length(p_encode), n_trials, n_replications))
sim_results_0_distant <- 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_0_pretrials <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))
normalized_echos_10_pretrials <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))
normalized_echos_30_pretrials <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))
normalized_echos_100_pretrials <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))

normalized_echos_0_adjacent <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))
normalized_echos_0_close <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))
normalized_echos_0_distant <- array(list(rep(NA, n_features)), dim = c(length(p_encode), n_trials, n_replications))


# all pair testing
# implementing the larger than rule
outcome[161:n_features] <- 1
context1[121:140] <- 1
context1[141:160] <- 0
context2[121:140] <- 0
context2[141:160] <- 1

items <- list(a,b,c,d,e,f)
all_pair_testing_array_context1 <- list()

for (i in 1:length(items)) {
  for (j in (i):length(items)) {
    if(i != j){
          #print(c(i,j))
          all_pair_testing_array_context1 <- append(all_pair_testing_array_context1,list(items[[i]]+items[[j]]+context1+outcome))
    }
  }
}

all_pair_testing_array_context2 <- list()

for (i in 1:length(items)) {
  for (j in (i):length(items)) {
    if(i != j){
          #print(c(i,j))
          all_pair_testing_array_context2 <- append(all_pair_testing_array_context2,list(items[[i]]+items[[j]]+context2))
    }
  }
}

all_pair_testing_array_context_combined <- append(all_pair_testing_array_context1, all_pair_testing_array_context2)

# create new probe items
test_bc_event_context1 <- b + c +  context + outcome 
recall_bc_event_context1 <- b + c + context1 + outcome

test_bd_event_context1 <- b + d +  context + outcome 
recall_bcd_event_context1 <- b + c + d + context1 + outcome

test_be_event_context1 <- b + e + context + outcome 
recall_bcde_event_context1 <- b + c + d + e + context1 + outcome

```

```{r}
# main test
#-----------------------Main Condition--------------------------------#

for (m in 1:length(model)){
  for (r in 1:n_replications) {
    training_array <- sample(training_array)

    for(i in 1:length(p_encode)) {
      # Memory is empty on first trial
      normalized_echo <- probe_memory(training_array[[1]], NULL, cue_features, model = model[m])
      expectancy <- expect_event(training_array[[1]], normalized_echo)
      memory <- learn(
        normalized_echo
        , training_array[[1]]
        , p_encode[i]
        , NULL
        , model = model[m]
      )
      
      normalized_echo <- probe_memory(test_bd_event, NULL, cue_features, model = model[m])
      expectancy <- expect_event(recall_bcd_event, normalized_echo)

      sim_results_0_close[[i,1,r]] <- expectancy
      normalized_echos_0_close[[i,1,r]] <- normalized_echo
      
      normalized_echo <- probe_memory(test_be_event, NULL, cue_features, model = model[m])
      expectancy <- expect_event(recall_bcde_event, normalized_echo)

      sim_results_0_distant[[i,1,r]] <- expectancy
      normalized_echos_0_distant[[i,1,r]] <- normalized_echo
      
      normalized_echo <- probe_memory(test_bc_event, NULL, cue_features, model = model[m])
      expectancy <- expect_event(recall_bc_event, normalized_echo)

      sim_results_0_adjacent[[i,1,r]] <- expectancy
      normalized_echos_0_adjacent[[i,1,r]] <- normalized_echo
      
       # Adjacent trials
      for(j in 2:(n_trials/3)) {
          normalized_echo <- probe_memory(training_array[[j]], memory, cue_features, model = model[m])
          expectancy <- expect_event(outcome, normalized_echo)
          memory <- learn(
            normalized_echo
            , training_array[[j]]
            , p_encode[i]
            , memory
            , model = model[m]
          )
          
      normalized_echo <- probe_memory(test_bd_event, memory, cue_features, model = model[m])
      expectancy <- expect_event(recall_bcd_event, normalized_echo)

      sim_results_0_close[[i,j,r]] <- expectancy
      normalized_echos_0_close[[i,j,r]] <- normalized_echo
      
      normalized_echo <- probe_memory(test_be_event, memory, cue_features, model = model[m])
      expectancy <- expect_event(recall_bcde_event, normalized_echo)

      sim_results_0_distant[[i,j,r]] <- expectancy
      normalized_echos_0_distant[[i,j,r]] <- normalized_echo
      
      normalized_echo <- probe_memory(test_bc_event, memory, cue_features, model = model[m])
      expectancy <- expect_event(recall_bc_event, normalized_echo)

      sim_results_0_adjacent[[i,j,r]] <- expectancy
      normalized_echos_0_adjacent[[i,j,r]] <- normalized_echo
      
      }
      
      for(j in ((n_trials/3) + 1):n_trials) {
          item <- unlist(sample(all_pair_testing_array_context1,1))
          normalized_echo <- probe_memory(item, memory, cue_features, model = model[m])
          expectancy <- expect_event(outcome, normalized_echo)
          memory <- learn(
            normalized_echo
            , item
            , p_encode[i]
            , memory
            , model = model[m]
          )
          
      normalized_echo <- probe_memory(test_bd_event_context1, memory, cue_features, model = model[m])
      expectancy <- expect_event(recall_bcd_event_context1, normalized_echo)

      sim_results_0_close[[i,j,r]] <- expectancy
      normalized_echos_0_close[[i,j,r]] <- normalized_echo
      
      normalized_echo <- probe_memory(test_be_event_context1, memory, cue_features, model = model[m])
      expectancy <- expect_event(recall_bcde_event_context1, normalized_echo)

      sim_results_0_distant[[i,j,r]] <- expectancy
      normalized_echos_0_distant[[i,j,r]] <- normalized_echo
      
      normalized_echo <- probe_memory(test_bc_event_context1, memory, cue_features, model = model[m])
      expectancy <- expect_event(recall_bc_event_context1, normalized_echo)

      sim_results_0_adjacent[[i,j,r]] <- expectancy
      normalized_echos_0_adjacent[[i,j,r]] <- normalized_echo
      
      }
      
    }
  }
}

#-----calculate mean and standard error of results objects-----

divide_by_n_rep <- function(x, n_rep = n_replications){x/n_rep}
# Sum of expectancies for each trial over all replications
sum_sim_results_0_adjacent <- apply(sim_results_0_adjacent, c(1,2), FUN = sum)
# Mean of expectancies for each trial over all replications
mean_sim_results_0_adjacent <- apply(sum_sim_results_0_adjacent, c(1,2), FUN = divide_by_n_rep)
# Sum of expectancies for each trial over all replications
sum_sim_results_0_close <- apply(sim_results_0_close, c(1,2), FUN = sum)
# Mean of expectancies for each trial over all replications
mean_sim_results_0_close <- apply(sum_sim_results_0_close, c(1,2), FUN = divide_by_n_rep)
# Sum of expectancies for each trial over all replications
sum_sim_results_0_distant <- apply(sim_results_0_distant, c(1,2), FUN = sum)
# Mean of expectancies for each trial over all replications
mean_sim_results_0_distant <- apply(sum_sim_results_0_distant, c(1,2), FUN = divide_by_n_rep)

# critical trial until x|cue >= 0.95 in each replication for every learning rate
# crit_trials_100_pretrials <- apply(sim_results_100_pretrials, 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_100_pretrials <- apply(crit_trials_100_pretrials, 1, FUN = std_err)
```

```{r}
#-----Performance plot

library("tidyverse")
adjacent_df <- as.data.frame(mean_sim_results_0_adjacent, row.names = p_encode)
adjacent_df_long <- adjacent_df %>% rownames_to_column("Learning_rate") %>% gather(Trial, MeanExp, -Learning_rate) %>% add_column(pair="adjacent")

close_df <- as.data.frame(mean_sim_results_0_close, row.names = p_encode)
close_df_long <- close_df %>% rownames_to_column("Learning_rate") %>% gather(Trial, MeanExp, -Learning_rate) %>% add_column(pair="close")

distant_df <- as.data.frame(mean_sim_results_0_distant, row.names = p_encode)
distant_df_long <- distant_df %>% rownames_to_column("Learning_rate") %>% gather(Trial, MeanExp, -Learning_rate) %>% add_column(pair="distant")

summary_df <- rbind(adjacent_df_long, close_df_long, distant_df_long)
summary_df <- summary_df %>% transform(., Trial = parse_number(Trial))

ggplot(data=summary_df, aes(x=Trial, y=MeanExp, color=Learning_rate)) +
  geom_line()+
  geom_point() + facet_grid("pair")

```


```{r}
#-----------------------0 Pretrials 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_a, NULL, cue_features, model = model[m])
      expectancy <- expect_event(outcome, normalized_echo)
      memory <- learn(
        normalized_echo
        , test_ax_event
        , p_encode[i]
        , NULL
        , model = model[m]
      )
      sim_results_0_pretrials[[i,1,r]] <- expectancy
      normalized_echos_0_pretrials[[i,1,r]] <- normalized_echo
      
       # A -> X trials
      for(j in 2:n_trials) {
          normalized_echo <- probe_memory(probe_a, memory, cue_features, model = model[m])
          expectancy <- expect_event(outcome, normalized_echo)
          memory <- learn(
            normalized_echo
            , test_ax_event
            , p_encode[i]
            , memory
            , model = model[m]
          )
          sim_results_0_pretrials[[i,j,r]] <- expectancy 
          normalized_echos_0_pretrials[[i,j,r]] <- normalized_echo
      }
    }
  }
}

#-----calculate mean and standard error of results objects-----#
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_0_pretrials <- apply(sim_results_0_pretrials, c(1,2), FUN = sum)
# Mean of expectancies for each trial over all replications
mean_sim_results_0_pretrials <- apply(mean_sim_results_0_pretrials, c(1,2), FUN = divide_by_n_rep)
# critical trial until x|cue >= 0.95 in each replication for every learning rate
crit_trials_0_pretrials <- apply(sim_results_0_pretrials, 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_0_pretrials <- apply(crit_trials_0_pretrials, 1, FUN = std_err)
```


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

#compute means over replications
x <- rowMeans(crit_trials_0_pretrials)
#compute standard errors over replications
y <- std_error_0_pretrials

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(1,2))
# 
# # memory encoding of the outcome over the trials
# plot(
#   1:n_trials
#   , 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:n_trials
#     , memory[, i]
#     , col = scales::alpha("black", 0.3)
#   )
# }
# 
# # memory encoding of cue A over the trials
# plot(
#   1:n_trials
#   , 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: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_0_pretrials[[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)
```