.packageName <- "EZ2"
`Data2EZ` <-
function (Pc,VRT,MRT,s=0.1)
# input:
#				Pc - Proportion correct
#				VRT - sample variance of the RT's
#				MRT - sample mean of the RT's
#				s - diffusion standard deviation
# returns:	Object with properties v, a, and Ter, containing EZ-estimates of these parameters
{
	#var L, x, v, a, y, MDT, Ter, s2, arr = [];
	p = Pc; 
      pow = function(a,b) a^b;
      logit = function(p) log(p/(1-p))
	if (p == 0)
	{
		stop("Oops, only errors!");
	}
	if (p == 0.5)
	{
		stop("Oops, chance performance!");
	}
	if (p == 1)
	{
		stop("Oops, only correct responses!");
	}
	s2     = s*s;
	L      = logit(p);
	x      = L*(L*p*p - L*p + p - 0.5) / VRT;
	v      = sign(p-0.5)*s*pow(x,(1/4));
	a      = s2*logit(p)/v;
	y      = -v*a/s2;
	MDT    = (a/(2*v))*(1-exp(y))/(1+exp(y));
	Ter    = ifelse(!missing(MRT),MRT - MDT,NA); # compute Ter only if MRT was provided
	return(list(v= v, a= a, Ter= Ter));
}

`EZ2` <-
function (pstart, ObsValPair, ..., method = "Nelder-Mead", control = list(), 
    hessian = FALSE) 
{
    tab = c(ObsValPair, ...)
    if (is.null(names(pstart)) | any(names(pstart) == "")) 
        stop("All elements of pstart should have names.")
    optim(pstart, EZ2.objf, method = method, control = control, 
        hessian = hessian, nms = names(pstart), ObsValPairs = tab, 
        grad = FALSE)
}

`EZ2.cmrt` <-
function (nu, z, a, s = 0.1) 
{
    .expr2 <- s^2
    .expr3 <- 4 * nu/.expr2
    .expr5 <- exp(.expr3 * a)
    .expr6 <- 2 * nu
    .expr7 <- z + a
    .expr10 <- exp(.expr6 * .expr7/.expr2)
    .expr12 <- .expr6/.expr2
    .expr14 <- exp(.expr12 * a)
    .expr18 <- exp(.expr6 * z/.expr2)
    .expr19 <- .expr5 + .expr10 - .expr14 - .expr18
    .expr23 <- 2 * .expr14 - 2 * .expr10
    .expr25 <- .expr19 * z + .expr23 * a
    .expr26 <- .expr25/nu
    .expr27 <- .expr14 - .expr18
    .expr28 <- .expr26/.expr27
    .expr30 <- -1 + .expr14
    .expr37 <- .expr10 * (2 * .expr7/.expr2)
    .expr41 <- .expr14 * (2/.expr2 * a)
    .expr45 <- .expr18 * (2 * z/.expr2)
    .expr60 <- .expr27^2
    .expr65 <- .expr30^2
    .expr68 <- .expr10 * .expr12
    .expr69 <- .expr18 * .expr12
    .expr73 <- 2 * .expr68
    .expr84 <- .expr14 * .expr12
    .value <- .expr28/.expr30
    .grad <- array(0, c(length(.value), 4), list(NULL, c("nu", 
        "z", "a", "Ter")))
    .grad[, "nu"] <- ((((.expr5 * (4/.expr2 * a) + .expr37 - 
        .expr41 - .expr45) * z + (2 * .expr41 - 2 * .expr37) * 
        a)/nu - .expr25/nu^2)/.expr27 - .expr26 * (.expr41 - 
        .expr45)/.expr60)/.expr30 - .expr28 * .expr41/.expr65
    .grad[, "z"] <- (((.expr68 - .expr69) * z + .expr19 - .expr73 * 
        a)/nu/.expr27 + .expr26 * .expr69/.expr60)/.expr30
    .grad[, "a"] <- (((.expr5 * .expr3 + .expr68 - .expr84) * 
        z + ((2 * .expr84 - .expr73) * a + .expr23))/nu/.expr27 - 
        .expr26 * .expr84/.expr60)/.expr30 - .expr28 * .expr84/.expr65
    .grad[, "Ter"] <- 0
    attr(.value, "gradient") <- .grad
    .value
}
`EZ2.cvrt` <-
function (nu, z, a, s = 0.1) 
{
    .expr2 <- -4 * nu
    .expr4 <- s^2
    .expr5 <- 2 * nu/.expr4
    .expr7 <- exp(.expr5 * a)
    .expr8 <- .expr2 * .expr7
    .expr10 <- 2 * z
    .expr13 <- exp(.expr10 * nu/.expr4)
    .expr14 <- -1 + .expr13
    .expr15 <- .expr8 * .expr14
    .expr17 <- 4 * nu/.expr4
    .expr19 <- exp(.expr17 * a)
    .expr20 <- .expr19 - .expr13
    .expr21 <- .expr15 * .expr20
    .expr22 <- a^2
    .expr25 <- 2 * (z + a)
    .expr28 <- exp(.expr25 * nu/.expr4)
    .expr29 <- 4 * .expr28
    .expr30 <- .expr29 * nu
    .expr31 <- .expr7 - 1
    .expr32 <- .expr31^2
    .expr33 <- .expr30 * .expr32
    .expr34 <- z^2
    .expr37 <- 8 * .expr28
    .expr38 <- .expr37 * nu
    .expr39 <- .expr38 * .expr32
    .expr40 <- .expr39 * a
    .expr43 <- 2 * .expr4
    .expr44 <- .expr43 * .expr7
    .expr45 <- .expr44 * .expr14
    .expr46 <- .expr45 * .expr31
    .expr48 <- -.expr7 + .expr13
    .expr49 <- .expr46 * .expr48
    .expr52 <- .expr4 * .expr32
    .expr54 <- 4 * z
    .expr57 <- exp(.expr54 * nu/.expr4)
    .expr58 <- -.expr19 + .expr57
    .expr59 <- .expr52 * .expr58
    .expr61 <- .expr21 * .expr22 - .expr33 * .expr34 + .expr40 * 
        z + .expr49 * a - .expr59 * z
    .expr62 <- .expr61/.expr32
    .expr63 <- nu^3
    .expr64 <- .expr62/.expr63
    .expr65 <- .expr7 - .expr13
    .expr66 <- .expr65^2
    .expr70 <- .expr7 * (2/.expr4 * a)
    .expr76 <- .expr13 * (.expr10/.expr4)
    .expr82 <- .expr19 * (4/.expr4 * a)
    .expr88 <- .expr28 * (.expr25/.expr4)
    .expr94 <- 2 * (.expr70 * .expr31)
    .expr132 <- .expr32^2
    .expr147 <- .expr66^2
    .expr150 <- .expr13 * .expr5
    .expr156 <- .expr28 * .expr5
    .expr159 <- 4 * .expr156 * nu * .expr32
    .expr166 <- 8 * .expr156 * nu * .expr32
    .expr191 <- .expr7 * .expr5
    .expr195 <- .expr19 * .expr17
    .expr203 <- 2 * (.expr191 * .expr31)
    .value <- .expr64/.expr66
    .grad <- array(0, c(length(.value), 4), list(NULL, c("nu", 
        "z", "a", "Ter")))
    .grad[, "nu"] <- ((((((.expr2 * .expr70 - 4 * .expr7) * .expr14 + 
        .expr8 * .expr76) * .expr20 + .expr15 * (.expr82 - .expr76)) * 
        .expr22 - ((4 * .expr88 * nu + .expr29) * .expr32 + .expr30 * 
        .expr94) * .expr34 + ((8 * .expr88 * nu + .expr37) * 
        .expr32 + .expr38 * .expr94) * a * z + (((.expr43 * .expr70 * 
        .expr14 + .expr44 * .expr76) * .expr31 + .expr45 * .expr70) * 
        .expr48 + .expr46 * (.expr76 - .expr70)) * a - (.expr4 * 
        .expr94 * .expr58 + .expr52 * (.expr57 * (.expr54/.expr4) - 
        .expr82)) * z)/.expr32 - .expr61 * .expr94/.expr132)/.expr63 - 
        .expr62 * (3 * nu^2)/.expr63^2)/.expr66 - .expr64 * (2 * 
        ((.expr70 - .expr76) * .expr65))/.expr147
    .grad[, "z"] <- ((.expr8 * .expr150 * .expr20 - .expr15 * 
        .expr150) * .expr22 - (.expr159 * .expr34 + .expr33 * 
        .expr10) + (.expr166 * a * z + .expr40) + (.expr44 * 
        .expr150 * .expr31 * .expr48 + .expr46 * .expr150) * 
        a - (.expr52 * (.expr57 * .expr17) * z + .expr59))/.expr32/.expr63/.expr66 + 
        .expr64 * (2 * (.expr150 * .expr65))/.expr147
    .grad[, "a"] <- (((.expr2 * .expr191 * .expr14 * .expr20 + 
        .expr15 * .expr195) * .expr22 + .expr21 * (2 * a) - (.expr159 + 
        .expr30 * .expr203) * .expr34 + ((.expr166 + .expr38 * 
        .expr203) * a + .expr39) * z + (((.expr43 * .expr191 * 
        .expr14 * .expr31 + .expr45 * .expr191) * .expr48 - .expr46 * 
        .expr191) * a + .expr49) - (.expr4 * .expr203 * .expr58 - 
        .expr52 * .expr195) * z)/.expr32 - .expr61 * .expr203/.expr132)/.expr63/.expr66 - 
        .expr64 * (2 * (.expr191 * .expr65))/.expr147
    .grad[, "Ter"] <- 0
    attr(.value, "gradient") <- .grad
    .value
}
`EZ2.mrt` <-
function (nu, z, a, s = 0.1, Ter) 
{
    .expr3 <- s * s
    .expr4 <- -2 * nu/.expr3
    .expr6 <- exp(.expr4 * a)
    .expr7 <- .expr6 - 1
    .expr8 <- a/.expr7
    .expr9 <- .expr8/nu
    .expr11 <- exp(.expr4 * z)
    .expr13 <- 1/nu
    .expr16 <- .expr13 * a
    .expr20 <- 2/.expr3
    .expr22 <- .expr6 * (.expr20 * a)
    .expr24 <- .expr7^2
    .expr27 <- nu^2
    .expr35 <- 1/.expr27
    .expr48 <- .expr6 * .expr4
    .value <- .expr9 * .expr11 - .expr13 * z - .expr16/.expr7 + 
        Ter
    .grad <- array(0, c(length(.value), 4), list(NULL, c("nu", 
        "z", "a", "Ter")))
    .grad[, "nu"] <- (a * .expr22/.expr24/nu - .expr8/.expr27) * 
        .expr11 - .expr9 * (.expr11 * (.expr20 * z)) + .expr35 * 
        z + (.expr35 * a/.expr7 - .expr16 * .expr22/.expr24)
    .grad[, "z"] <- .expr9 * (.expr11 * .expr4) - .expr13
    .grad[, "a"] <- (1/.expr7 - a * .expr48/.expr24)/nu * .expr11 - 
        (.expr13/.expr7 - .expr16 * .expr48/.expr24)
    .grad[, "Ter"] <- 1
    attr(.value, "gradient") <- .grad
    .value
}
`EZ2.objf` <-
function(p,ObsValPairs,nms=names(p),grad=FALSE)
{
    #
    # the basic idea underlying this function is:
    #   ObsValPair = 0.3 ~ EZ2.vrt(nu=nu,z=a-z,a=a)
    #   with(list(nu=.1,z=0.09,a=0.2), eval(ObsValPair[[3]])
    #
    names(p) = nms
    p = as.list(p)
    f = sapply(ObsValPairs, function(.x) {
        y = with(p, eval(.x[[3]]))
        z = rep(0, length(nms))
        names(z) = nms
        .g = attr(y, "gradient")[1, ]
        z[names(.g)] = .g
        c(pred = y, z)
    })
    e = sapply(ObsValPairs, function(x) eval(x[[2]])) - f[1, ]
    s = 1e+06 * sum(e^2)
    if (grad) {
        df = 1e+06 * f[-1, ]
        .grad = -2 * df %*% e
        attr(s, "gradient") = drop(.grad)
    }
    return(s)
}

`EZ2.objf.old` <-
function(p,nms,ObsValPairs,grad=FALSE)
{
    #
    # the basic idea underlying this function is:
    #   ObsValPair = c(func="EZ2.vrt", obs=0.3,v= ~v, z=~a-z, a=~a)
    #   args = lapply(ObsValPair[-(1:2)],function(exp) eval(exp[[length(exp)]]))
    #   do.call("EZ2.vrt",args)
    #
###   eval(parse(text=paste(nms,'=',p,collapse=';'))) # define the parameters ---> changed to using 'with'
    names(p)=nms # nlm removes names 
    p=as.list(p)
    f = sapply(ObsValPairs, #1, 
		function(x) do.call(x[[1]],lapply(x[-(1:2)],function(expr) with(p,eval(expr[[length(expr)]]))))
        )
    e = unlist(sapply(ObsValPairs,function(x)x[[2]]))-unlist(f) #unlist(ObsValPairs[,2])-unlist(f)
    s = sum(e^2)
    if(grad & FALSE){ # this is currently not used...
	.grad = -2*t(df)%*%e
	attr(s,"gradient") = .grad
    }
    return(s)
}

`EZ2.pe` <-
function (nu, z, a, s = 0.1) 
{
    .expr3 <- s * s
    .expr4 <- -2 * nu/.expr3
    .expr6 <- exp(.expr4 * a)
    .expr8 <- exp(.expr4 * z)
    .expr9 <- .expr6 - .expr8
    .expr10 <- .expr6 - 1
    .expr12 <- 2/.expr3
    .expr14 <- .expr6 * (.expr12 * a)
    .expr20 <- .expr10^2
    .expr27 <- .expr6 * .expr4
    .value <- .expr9/.expr10
    .grad <- array(0, c(length(.value), 4), list(NULL, c("nu", 
        "z", "a", "Ter")))
    .grad[, "nu"] <- -((.expr14 - .expr8 * (.expr12 * z))/.expr10 - 
        .expr9 * .expr14/.expr20)
    .grad[, "z"] <- -(.expr8 * .expr4/.expr10)
    .grad[, "a"] <- .expr27/.expr10 - .expr9 * .expr27/.expr20
    .grad[, "Ter"] <- 0
    attr(.value, "gradient") <- .grad
    .value
}
`EZ2.vrt` <-
function (nu, z, a, s = 0.1) 
{
    .expr4 <- s^2
    .expr5 <- -2 * nu/.expr4
    .expr7 <- exp(.expr5 * z)
    .expr8 <- .expr7 - 1
    .expr9 <- .expr8^2
    .expr10 <- -nu * .expr9
    .expr11 <- a^2
    .expr13 <- 4 * nu
    .expr14 <- .expr13 * .expr8
    .expr16 <- .expr10 * .expr11 - .expr14 * .expr11
    .expr18 <- exp(.expr5 * a)
    .expr19 <- .expr18 - 1
    .expr20 <- .expr19^2
    .expr23 <- -3 * nu
    .expr25 <- .expr13 * z
    .expr26 <- .expr25 * a
    .expr29 <- .expr23 * .expr11 + .expr26 + .expr4 * a
    .expr31 <- .expr29 * .expr8 + .expr26
    .expr35 <- .expr16/.expr20 + .expr31/.expr19 - .expr4 * z
    .expr36 <- nu^3
    .expr38 <- 2/.expr4
    .expr40 <- .expr7 * (.expr38 * z)
    .expr53 <- .expr18 * (.expr38 * a)
    .expr57 <- .expr20^2
    .expr61 <- 4 * z * a
    .expr80 <- .expr13 * a
    .expr82 <- .expr7 * .expr5
    .expr98 <- 2 * a
    .expr103 <- .expr18 * .expr5
    .value <- .expr35/.expr36
    .grad <- array(0, c(length(.value), 4), list(NULL, c("nu", 
        "z", "a", "Ter")))
    .grad[, "nu"] <- (((nu * (2 * (.expr40 * .expr8)) - .expr9) * 
        .expr11 - (4 * .expr8 - .expr13 * .expr40) * .expr11)/.expr20 + 
        .expr16 * (2 * (.expr53 * .expr19))/.expr57 + (((.expr61 - 
        3 * .expr11) * .expr8 - .expr29 * .expr40 + .expr61)/.expr19 + 
        .expr31 * .expr53/.expr20))/.expr36 - .expr35 * (3 * 
        nu^2)/.expr36^2
    .grad[, "z"] <- ((.expr80 * .expr8 + .expr29 * .expr82 + 
        .expr80)/.expr19 - (nu * (2 * (.expr82 * .expr8)) * .expr11 + 
        .expr13 * .expr82 * .expr11)/.expr20 - .expr4)/.expr36
    .grad[, "a"] <- ((.expr10 * .expr98 - .expr14 * .expr98)/.expr20 - 
        .expr16 * (2 * (.expr103 * .expr19))/.expr57 + (((.expr23 * 
        .expr98 + .expr25 + .expr4) * .expr8 + .expr25)/.expr19 - 
        .expr31 * .expr103/.expr20))/.expr36
    .grad[, "Ter"] <- 0
    attr(.value, "gradient") <- .grad
    .value
}
`EZ2batch` <-
function (pstart, ObsValPair, ..., data, nrestart = 1, method = "Nelder-Mead", 
    control = list(), hessian = FALSE) 
{
    mdl = c(ObsValPair, ...)
    t(apply(data, 1, function(x) {
        attach(as.list(x), warn = FALSE)
        fit <- EZ2(pstart, mdl, method = method, control = control, 
            hessian = hessian)
        for (i in 1:nrestart) fit <- EZ2(fit$par, mdl, method = method, 
            control = control, hessian = hessian)
        unlist(fit)
    }))
}

