# 15 February, 2008
# based on P2/p4.r

library(ggplot2)
library(reshape)
library(lme4)

rm(list=ls())
load("MM.rda")

# Select trials with MS in CTI; only Exp 1 to 4
ix <- which( d$TYP == "no MS" | d$mt <= 0 | d$EXP == "VVb")
d <- d[-ix, ]
d[] <- lapply(d,function(x) x[drop=TRUE])
d$TYP <- relevel(d$TYP, ref="single MS")
d$BINS <- cut(d$mt, breaks=c(0, 160, 320, 480, 640, 800, 960))
d$CUE <- ifelse(d$EXP =="VV" | d$EXP == "VA", "visual", "auditory")

### CHUNK 1: Figure 4 of LKRE
d2 <- d[d$TYP !="last MS" & is.finite(d$BINS), ]
d2.rs <- melt(d2, id=c("ID", "EXP", "CV", "MTC", "BINS"), measure="mt")

table <- cast(d2.rs, EXP + BINS ~  CV + MTC + variable, function(x) c(N=length(x)) )
table$N <- rowSums(table[,3:6])
table$p <- table[ ,6]/table$N
table$x <- (table$p - 0.16)/0.48   # x is prob that MS direction corresponds to direction of covert attention shift
table <- data.frame(table)
table$bins <- as.numeric(table$BINS)
m <- qplot(bins, x, data=table, geom=c("point", "line"), facets= . ~ EXP)
m <- m + scale_x_continuous("MS Time of Occurrence [ms]", breaks=c(1:6), labels=c("80","240","400","560","720", "880"))
m <- m + scale_y_continuous("Probability(MS Direction == Attention)" )
m

#########################################

# CHUNK2: Table 2 of LKRE
#d3 <- d[d$TYP !="last MS" & d$mt > 160 & d$mt <= 800, ]
d3 <- d[d$TYP !="last MS" & d$mt > 0 & d$mt <= 800, ]
d3.rs <- melt(d3, id=c("ID", "CUE", "CV", "MTC"), measure="mt")

table <- cast(d3.rs, CUE ~  CV + MTC + variable, function(x) c(N=length(x)) )
table$N <- rowSums(table[,2:5])
table[ ,2:5]/table$N



# MODEL COMPUTATIONS
# parameters
a <- 1.0       # prob of cue; note a=1 for model w/o neutral cues
b <- 0.80      # cue validity 
k <- b         # probability matching

p.ma <- c(1.0, 0.75, 0.5)  # probability that microsaccade goes with attention

# prob of model states
state <- vector()

for (i in 1:3) {
x <- p.ma[i]

state[1] <- a*b*k*x
state[2] <- a*b*k*(1-x)
state[3] <- a*b*(1-k)*x
state[4] <- a*b*(1-k)*(1-x)

state[5] <- a*(1-b)*k*x
state[6] <- a*(1-b)*k*(1-x)
state[7] <- a*(1-b)*(1-k)*x
state[8] <- a*(1-b)*(1-k)*(1-x)

state[9]  <- (1-a)*0.50*x
state[10] <- (1-a)*0.50*(1-x)
state[11] <- (1-a)*0.50*x
state[12] <- (1-a)*0.50*(1-x)

p.v.con <- state[1]+state[4]
p.v.inc <- state[2]+state[3]

p.i.con <- state[6]+state[7]
p.i.inc <- state[5]+state[8]

p.n.con <- state[9]+state[12]
p.n.inc <- state[10]+state[11]

print(p <- c(p.v.con, p.v.inc, p.i.con, p.i.inc, p.n.con, p.n.con))
}