# 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
ix <- which( d$TYP == "no MS" | d$mt <= 0)
d <- d[-ix, ]
d[] <- lapply(d,function(x) x[drop=TRUE])
d$TYP <- relevel(d$TYP, ref="single MS")

exp <- c("VV", "VA", "AV", "AA", "VVb")
for (i in 1:5) {
d.exp <- d[d$EXP==exp[i] & d$TYP !="last MS" & d$mt > 250 & d$mt <= 400, ]
#d.exp <- d[d$EXP==exp[i] & d$TYP !="last MS", ]
d.exp[] <- lapply(d.exp,function(x) x[drop=TRUE])

d.exp.rs <- melt(d.exp, id=c("ID", "CV", "MTC", "TYP"), measure="rt")
table <- cast(d.exp.rs, CV+MTC ~ variable, function(x) c(M=round(mean(x)), N=length(x)) )
table$prob1 <- table$rt_N/sum(table$rt_N)
table

x <- (table$CV=="valid")*(.64-table$prob1)/.48
x <- x[x>0]
print(paste(exp[i],": x=", round(x[1], 2), "N=", sum(table$rt_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))
}


# Compute prob directly for observed probs.(only old model!)
x <- (table$CV[c(1,2,5,6)]=="valid")*(.64-table$prob2[c(1,2,5,6)])/.48
x <- x[x>0]
x