---
title: "Statistical Properties of the N-Back Sequences"
output:
  html_notebook: default
  pdf_document: default
editor_options: 
  chunk_output_type: inline
---

# Problems

Statistical properties of n-back sequences bias behaviors. These bias, under specified structure, allows multiple cognitive strategies, producing heterogeneous behavior in a "gold standard" cognitive task.

# Gaps

- Unclear how to parameterize interesting variations for sequence generation
- How do we model these multiple strategies (which requires identifying which sequence variations matter) 
    - local vs. global properties, which one matters the most?
    - Local:  lumpiness, short sequence patterns -> could be exploited by “reactive”/automaticity 
    - Global:  No lures, large vocabulary -> pattern repeats implies a target


## Formulating Generating the N-Back Sequences as a CSP instance

$P=\langle V,D,C,W\rangle$

$V=\{x_N,x_{T},x_{T,local},x_L,x_{L,local},x_V,x_U,x_S,x_{S,local},x_G\}$

$D=\{\}$


Constraints:

$$
\\

x_n = N, W_n = 1 - |10 \times dnorm(x_n-N,sd=4)|

\\\\

x_t = T \times trials, W_t = 1 - |10\times dnorm(T\times trials-x_t,sd=4)|

\\\\

x_{tl} = {T \times w \over trials}, W_{tl} = 1 - |10\times dnorm(x_{tl} - {T \over trials} \times w,sd=4)|

\\\\

x_{l} = L \times trials
\\\\

x_{ll} = L \times w
\\\\

x_{v} = |V|
\\

x_{ul} = w
\\\\

x_{s} = {trials \over |V|}
\\\\

x_{sl} = max(1, {w \over |V|})
\\\\

x_{g} = {trials \over w}

\\\\

x_{vl} = min(|V|, w)
$$

```{r libraries, message=FALSE, include=FALSE, paged.print=FALSE}
library(ggplot2)
library(tidyverse)
library(stringi)
library(plsRglm)
library(plsdof)
library(caret)
```

```{r params}
load('./data/CL2015.RData')

window_size <- 8
```



```{r history}

with_lures <- function(stimulus, stimulus_type, n) {
  sapply(1:length(stimulus), function(i) {
    lures <- c(as.character(stimulus[i-n-1]), as.character(stimulus[i-n+1]))
    are_valid_trials <- i>n && all(!is.na(c(lures,stimulus[i])))
    ifelse(are_valid_trials && stimulus[i] %in% lures,
           "lure", 
           as.character(stimulus_type[i]))
  })
}

with_history <- function(stimuli, length=16, fixed=F) {
  seq <- paste(stimuli, collapse = '')
  
  sapply(1:length(stimuli), function(i) {
    stri_reverse(str_sub(seq, max(1,i-length+1), i))
  })
  #ifelse(fixed, h[str_length(h)==size], h)
}

# $x_{s,local}$
with_skewness <- function(history) {
  sapply(history, function(h) {
    freqs <- table(unlist(str_split(h,"")))
    sum(sort(freqs, decreasing = T)[1:2]) - 1
  })
}

# $x_{u,local}$
with_lumpiness <- function(history) {
  sapply(history, function(h) {
    freqs <- table(unlist(str_split(h,"")))
    max(freqs) - 1
  })
}


with_targets_ratio <- function(stimulus_type, length=16) {
  sapply(1:length(stimulus_type), function(i) {
    trials <- stimulus_type[max(1,i-length):i]
    length(trials[trials=="target"]) / length(trials)
  })
}

with_lures_ratio <- function(stimulus_type, length=16) {
  sapply(1:length(stimulus_type), function(i) {
    trials <- stimulus_type[max(1,i-length):i]
    length(trials[trials=="lure"]) / length(trials)
  })
}

#TODO change to list column workflow with broom for model fitting and evaluating the fits
# duh! we are using list nested insided a tibble, so put all new columns in a new list column
# instead of adding a new column for each.
NB2 <- NB %>%
  group_by(participant, condition, block) %>%
  nest() %>% unnest(data) %>%
  mutate(n = ifelse(condition=='2-back',2,3)) %>%
  mutate(stimulus_type = with_lures(stimulus, stimulus_type, n)) %>%
  mutate(history = with_history(stimulus, window_size)) %>%
  mutate(x_sl = with_skewness(history)) %>%
  mutate(x_ul = with_lumpiness(history)) %>%
  mutate(x_t = with_targets_ratio(stimulus_type, window_size)) %>%
  mutate(x_l = with_lures_ratio(stimulus_type, window_size)) %>%
  ungroup()
  
pca <- prcomp(~x_sl+x_ul+x_t+x_l, NB2, center = TRUE,scale. = TRUE, na.action=na.exclude)
NB2 <- NB2 %>% mutate(pc1=pca$x[,'PC1'], pc2=pca$x[,'PC2'])

# caret
library(caret)
# Compile cross-validation settings


any(is.na(NB2))
NB2 <- na.omit(NB2)

# set.seed(100)
# trainingfold <- createMultiFolds(NB2@correct, k = 5, times = 10)
# 
# # PLS
# mod1 <- train(correct ~ ., data = NB2[,c("correct","x_sl","x_ul","x_t","x_l")],
#  method = "pls",
#  metric = "Accuracy",
#  tuneLength = 20,
#  trControl = trainControl("repeatedcv", index = trainingfold, selectionFunction = "oneSE"),
#  preProc = c("zv","center","scale"))
#  
# # Check CV
# plot(mod1)


plsResult <- plsR(rt ~ ., data=NB2[,c("rt","x_sl","x_ul","x_t","x_l")],3)
plsResult <- plsR(correct ~ ., data=NB2[,c("correct","x_sl","x_ul","x_t","x_l")],3)


plsResult
```

<!--chapter:end:ccn2019.rev2.Rmd-->

---
title: "Evaluating N-Back Sequences"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(tidyverse)
library(ggplot2)
library(stringi)
library(GA)
library(dbscan)
library(inspectdf)

load('./data/CL2015.RData')
```

### Variables
- $T$ number of targets
- $L$ number of lures
- $S$ Skewness score
- $U$ Uniformity (!repetition)
- $RT_{mean}$
- $Accuracy_{mean}$
- $dprime$
- $criterion$

## Constraints

- fixed number of targets
- fixed number of lures (a.k.a, foils)
- uniform distribution of choices
- controlled local lumpiness


Each constraint is an up side down quadratic function to be minimized.

```{r, eval=F}
targets_cost <- function(x) 1.0 - 10*dnorm(x,mean=0,sd=4)
lures_cost <- function(x) 1.0 - 10*dnorm(x,mean=0,sd=4)
skewness_cost <- function(x, choices) {
  uniform_ratio <- length(x) / length(choices)
  deviation_from_uniform <- setNames(vector('numeric', length(choices)), choices)
  for (c in choices) {
    deviation_from_uniform[c] = abs(length(x[x==c]) - uniform_ratio)
  }
  #TODO convert to gaussian loss
  max(deviation_from_uniform)
}

lumpiness_cost <- function(x, choices) {
  #trials = len(seq)
  #freqs = [float(seq.count(c)) for c in choices]
  #ralph_skewed = sum(heapq.nlargest(int(len(choices) / 2), freqs)) > (trials * 2 / 3)
  #return ralph_skewed
  NA
}

#merged_cost <- function(x) targets_fitness(x) + lures_fitness(x) + skewness_fitness(x)
#GA <- ga(type = "real-valued", fitness = merged_cost, lower = -10, upper = 10)
#plot(GA)
```


```{r}

with_lures <- function(condition, stim, stim_type, history = NA) {
  sapply(1:length(stim), 
    function(i) {
      switch(as.character(condition[i]),
             "2-back" = {
               ifelse(
                 stim[i]==stri_sub(history[i],-2,-2) ||  stim[i]==stri_sub(history[i],-4,-4), 
                 'lure', 
                 as.character(stim_type[i])
               )},
             "3-back" = {
               ifelse(
                 stim[i]==stri_sub(history[i],-3,-3) ||  stim[i]==stri_sub(history[i],-5,-5), 
                 'lure', 
                 as.character(stim_type[i])
               )}
      )

    })
}

with_targets_ratio <- function(stimulus_type, history) {
  sapply(1:length(history), function(i) {
    trials <- stimulus_type[(i-stri_length(history[i])):i]
    trials <- unlist(trials, use.names=FALSE)
    length(trials[trials=="target"])
  })
}

with_lures_ratio <- function(stimulus_type, history) {
  sapply(1:length(history), function(i) {
    trials <- stimulus_type[(i-str_length(history[i])):i]
    trials <- unlist(trials, use.names=FALSE)
    length(trials[trials=="lure"])
  })
}

with_lumpiness_score <- function(stimulus, history) {
  sapply(1:length(history), function(i) {
    trials <- stimulus[(i-str_length(history[i])):i]
    trials <- unlist(trials, use.names=FALSE)
    max(table(trials)) - 1
  })
}

with_lag <- function(stimulus, history) {
  # find last occurance the of stimulus
}

with_skewness_score <- function(stimulus, history) {
  sapply(1:length(history), function(i) {
    trials <- stimulus[(i-str_length(history[i])):i]
    trials <- unlist(trials, use.names=FALSE)
    sum(sort(table(trials), decreasing = T)[1:2]) - 1
  })
}

with_history <- function(stims, size=16) {
  res <- c('')
  for (i in 2:length(stims)) {
    res[i] <- stri_sub(paste(res[i-1], stims[i], sep=''),from=-size,length=size)
  }
  #res <- ifelse(stri_length(res)==size, res, NA)
  res
}

normalize_scores <- function(targets_ratio, lures_ratio, skewness, lumpiness) {
  #TODO
  sapply(1:length(targets_ratio), function(i) 0)
}

window_size <- 8

NB_modified <- NB %>%
  group_by(participant, condition, block) %>%
  mutate(history = with_history(stimulus, size=window_size)) %>%
  #mutate(stimulus_type = map_chr(.x=stimulus, stim_type=stimulus_type, history=history,.f=with_lures))
  mutate(stimulus_type_2 = with_lures(condition, stimulus, stimulus_type, history)) %>%
  mutate(targets = with_targets_ratio(stimulus_type_2, history)) %>%
  mutate(lures = with_lures_ratio(stimulus_type_2, history)) %>%
  mutate(skewness = with_skewness_score(stimulus, history)) %>%
  mutate(lumpiness = with_lumpiness_score(stimulus, history)) %>%
  filter(stri_length(history)==window_size) %>%
  mutate(correct = ifelse(stimulus_type=='burn-in',NA,correct)) %>%
    #normalize_scores(targets_ratio, lures_ratio, skewness, lumpiness) %>%
  ungroup()

pca <- prcomp(NB_modified[,c('targets','lures','skewness','lumpiness')], center = TRUE,scale. = TRUE)

NB_modified <- NB_modified %>% mutate(pc1=pca$x[,'PC1'], pc2=pca$x[,'PC2'])


## participant-level averaged NB, a single row represent an observation for a single subject
## in a single condition
NB_avg <- NB_modified %>%
  group_by(participant, condition) %>%
  mutate(correct = ifelse(stimulus_type=='burn-in',NA,correct)) %>%
  summarise(
    targets=sum(targets), 
    lures=sum(lures), 
    skewness=sum(skewness), 
    lumpiness=sum(lumpiness), 
    rt = mean(rt, na.rm=T), 
    correct=sum(correct,na.rm=T)/90) %>%
  ungroup()

# print
# NB_modified %>%
#   filter(participant=='P1') %>%
#   View()
#   
  

fit <- lm(correct ~ t * s * u * l * d, NB_modified)

```


```{r}

# DBSCAN Clustering (RT+ACCURACY against skewness)
NB_avg <- NB_avg %>%
  mutate(cluster = dbscan(cbind(correct,rt), eps = 0.3, minPts = 3)$cluster)

NB_avg %>%
  ggplot(aes(skewness, correct, color=factor(cluster))) +
  ggtitle(" clusters (window = 16 trials)", "NOTE: each point is a single participant") +
  geom_point(alpha=0.3) +
  #geom_smooth(method='lm', se = F) +
  facet_wrap(~condition)

```


```{r}
## single-subject figures
NB_modified %>%
  ggplot(aes(t,s,color=correct)) +
  geom_jitter() +
  geom_point() +
  stat_summary(fun.y="mean")

NB_modified %>%
  inspect_cor(show_plot = T)

NB_avg %>%
  inspect_cor(show_plot = T)

NB_modified %>%
  ggplot(aes(rt,correct,color=u)) +
  geom_jitter() +
  geom_point() +
  stat_summary(fun.y="mean")

NB_modified %>%
  filter(!is.na(correct)) %>%
  ggplot(aes(jitter(s),jitter(t),color=correct)) +
    geom_jitter() +
    geom_point(alpha=0.1)

NB_modified %>%
  filter(!is.na(correct)) %>%
  ggplot(aes(jitter(s),jitter(u),color=correct)) +
    geom_jitter() +
    geom_point() +
    facet_wrap(~condition)

# rt/accuracy and lures
NB_modified %>%
  filter(!is.na(correct)) %>%
  ggplot(aes(jitter(l),rt,color=correct,alpha=0.01)) +
    geom_jitter() +
    geom_point(shape=16) +
    geom_smooth(method="lm",se = F) + 
    facet_wrap(~condition, scales="free")


NB_modified %>%
  filter(!is.na(correct)) %>%
  ggplot(aes(pc1,pc2,color=correct)) +
    geom_point() +
    geom_smooth(method="lm",se = F) + 
    facet_wrap(~condition, scales="free")

NB_modified %>%
  filter(!is.na(correct)) %>%
  ggplot(aes(pc1,rt,color=correct)) +
    geom_point(alpha=0.3) +
    geom_smooth(method="lm",se = F) + 
    facet_wrap(~condition, scales="free")

NB_modified %>%
  filter(!is.na(correct)) %>%
  ggplot(aes(pc1,correct,color=correct)) +
    geom_point() +
    geom_smooth(method="lm",se = F) + 
    facet_wrap(~condition, scales="free")


```

## TODO
  - data %>% mutate(constraint1=fitness1(history), constrain2=fitness2(history),
  - constraint3=fitness(history))
  - kmeans(NB)
  - ggplot(kmeans_clusters$accuracy) 
  - ggplot(kmeans_clusters$rt)
  - 




```{python}
a=2
p
```

<!--chapter:end:ccn2019.Rmd-->

---
title: "PLS Training"
output: html_notebook
---

PLS:


```{r}
#detach("package:MASS","plsdof") # to avoid conflict with dplyr::select
library(tidyverse)
library(pls)

## 1. load sample data
#data <- read.csv("http://wiki.q-researchsoftware.com/images/d/db/Stacked_colas.csv")

rm(NB)
load("./data/CL2015.RData")
data <- NB
str(data)

## 2. clean data (remove brand and URLID)
data <- data %>% 
  mutate(n=ifelse(condition=='2-back', 2, 3)) %>%
  select(-condition,
         -stimulus,
         -block,
         -trial)
# %>%
#   rename(
#          ev.participant=participant,
#          ev.n=n,
#          ev.block=block,
#          ev.stimulus_type=stimulus_type,
#          rv.choice=choice,
#          rv.rt=rt,
#          rv.correct=correct
#          )

## 3. use cross validatation to find the optimal number of dimensions
pls.model = plsr(rt ~ ., data = data, validation = "CV")

## 3.1. find the model with lowest cv error
cv <- RMSEP(pls.model)
best_dims <- which.min(cv$val[estimate = "adjCV", , ]) - 1

## 4. rebuild the model
pls.model <- plsr(rt ~ ., data = data, ncomp = best_dims)

## 5. Sort, and visualize top coefficients
coefs <- coef(pls.model)

barplot(sort(coefs[,1,1], decreasing = T)[1:4])
```


```{r simulate}
X <- matrix(rnorm(1100), 100, 11)
Y <- matrix(rnorm(400), 100, 4)

pls.model <- plsr(Y ~ X, validation = "CV")

cv <- RMSEP(pls.model)
best_dims <- which.min(cv$val[estimate = "adjCV", , ]) - 1
pls.model <- plsr(Y ~ X, ncomp = best_dims)
coefs <- sort(coef(pls.model)[,1,1], decreasing = T)

barplot(coefs)

```


```{r cca-simulate}
X <- matrix(rnorm(1100), 100, 11)
Y <- matrix(rnorm(400), 100, 4)

M <- cor(cbind(X,Y))
corrplot(M, method="ellipse", order="hclust", addrect=2, addCoef.col="black")
cc <- cancor(X, Y)

#NB: cc <- cancor(cbind(rt,correct, accuracy) ~ xt + xl + xtl, data = data)

```


```{r plsrglm}
rm(list = ls())
library(plsRglm)

data(Cornell)
df <- Cornell
x <- subset(df, select = -c(Y))
y <- df$Y
## K is the number of folds in CV, and nt is the maximum number of components, 
#cv.modpls<-cv.plsRglm(dataY=y,dataX=x ,nt=10,modele="pls-glm-logistic",K=8)

modpls <- plsRglm(dataY = y,dataX = x, nt = 10, modele = "pls-glm-logistic", sparse=TRUE,sparseStop=TRUE)
res.cv.modpls<-cvtable(summary(cv.modpls))

res6<-plsR(Y~.,data=Cornell, nt=6, typeVC="missing", pvals.expli=TRUE)

```

<!--chapter:end:pls_playground.Rmd-->

