#########################
###      ANALYSES    ####
#########################

path <- "C:/datafiles/"
path <- "C:/"

AESdata <- source(paste(path,"AES.dmp",sep=""))$value

Y <- AESdata[,2:25]
v <- AESdata[,1]
X <- rbind(as.matrix(AESdata[v==1,28:51]),as.matrix(AESdata[v==2,52:75]))
u <- c(F,F,F,F,F,T,T,T,T,T,T,F,T,T,T,T,F,T,T,T,T,F,T,T)

#########################
###      tau1         ###
#########################

T1 <- compute.T1(X,Y)
ET1 <- approximate.ET1(X,v,u,m)
MT1 <- compute.MT1(X,v,u)
tau1 <- (T1 - ET1)/(MT1 - ET1)

#########################
###      tau2         ###
#########################

T2 <- compute.T2(X,Y)
ET2 <- approximate.ET2(X,v,u,m)
MT2 <- compute.MT2(X,v,u)
tau2 <- (T2 - ET2)/(MT2 - ET2)

#########################
### tau1* and tau2*   ###
#########################

G <- compute.G(X,Y,v)
MG <- compute.MG(X,v,u)
tau1.star <- tau1 + 0.5 * G/MG
tau2.star <- tau2 + 0.5 * G/MG

#########################
###  R functions      ###
#########################

compute.T1 <- function(X,Y){
  N <- nrow(X)
  J <- ncol(X)
  M <- array(0,c(N,N,J))
  for (j in 1:J){M[,,j] <- outer(Y[,j],Y[,j],"==")}
  XX <- array(0,c(N,N,J))
  for (j in 1:J){XX[,,j] <- outer(X[,j],X[,j])}
  T1 <- apply(M * (1 - XX),c(1,2),sum)
  diag(T1) <- 0
  return(T1)
}

exp.T1.function <- function(X.v,X.w,J,U,m){
  P.v <- X.v/J
  P.w <- X.w/J
  Q.v <- (1-P.v)/(m-1)
  Q.w <- (1-P.w)/(m-1)
  E.U <- U * ((m-2) * Q.v * Q.w + P.v * Q.w + Q.v * P.w)
  E.C <- (J-U) * (m-1) * Q.v * Q.w
  return(E.C + E.U)
}

approximate.ET1 <- function(X,v=one.version,u=all.common,m=4){
 J <- ncol(X)
 N <- nrow(X)
 one.version <- rep(1,N)
 all.common <- rep(F,J)
 V <- outer(v,v,"==")
 X. <- X %*% matrix(1,nrow=J)
 length.u <- length(u[u==F])

 ET1.same <- matrix(0,J+1,J+1)
 for(i in 0:J) for(j in 0:J) ET1.same[i+1,j+1] <- exp.T1.function(i,j,J,0,m)

 ET1.diff <- matrix(0,J+1,J+1)
 for(i in 0:J) for(j in 0:J) ET1.diff[i+1,j+1] <- exp.T1.function(i,j,J,length.u,m)

 ET1 <- matrix(0,N,N)
 for (i in 1:N) for(j in 1:N){
   ET1[i,j] <- ifelse(V[i,j], ET1.same[X.[i]+1,X.[j]+1],ET1.diff[X.[i]+1,X.[j]+1])
 }
 diag(ET1) <- 0
 return(ET1)
}
max.T1.function <- function(L,K,J,U){
   if (K > L){G <- L; L <- K; K <- G}
   A <- min(U,L,J-K)
   B <- min(U-A,K,J-L)
   C <- min(J-A-B,J-L-B,J-K-A)
   return(A + B + C)
}

compute.MT1 <- function(X,v=one.version,u=all.common){
  J <- ncol(X)
  N <- nrow(X)
  one.version <- rep(1,N)
  all.common <- rep(F,J)
  V <- outer(v,v,"==")
  X. <- X %*% matrix(1,nrow=J)
  length.u <- length(u[u==F])

  MT1.diff <- matrix(0,J+1,J+1)
  for(i in 0:J) for(j in 0:J) MT1.diff[i+1,j+1] <- max.T1.function(i,j,J,length.u)

  MT1.same <- matrix(0,J+1,J+1)
  for(i in 0:J) for(j in 0:J) MT1.same[i+1,j+1] <- max.T1.function(i,j,J,0)

  MT1 <- matrix(0,N,N)
  for (i in 1:N) for(j in 1:N){
    MT1[i,j] <- ifelse(V[i,j], MT1.same[X.[i]+1,X.[j]+1],MT1.diff[X.[i]+1,X.[j]+1])
  }
  return(MT1)
}
compute.T2 <- function(X,Y){
  N <- nrow(X)
  J <- ncol(X)
  M <- array(0,c(N,N,J))
  for (j in 1:J){M[,,j] <- outer(Y[,j],Y[,j],"==")}
  XX. <- array(0,c(N,N,J))
  for (j in 1:J){XX.[,,j] <- outer(X[,j],1-X[,j])}
  X.X. <- array(0,c(N,N,J))
  for (j in 1:J){X.X.[,,j] <- outer(1-X[,j],1-X[,j])}
  T2 <- t(apply(M * (2 * XX. + X.X.),c(1,2),sum))
  diag(T2) <- 0
  return(T2)
}

exp.T2.function <- function(X.v,X.w,J,U,m){
  P.v <- X.v/J
  P.w <- X.w/J
  Q.v <- (1-P.v)/(m-1)
  Q.w <- (1-P.w)/(m-1)
  E.U <- U * ((m-2) * 0.5 * Q.v * Q.w + Q.v * P.w)
  E.C <- (J-U) * 0.5 * (m-1) * Q.v * Q.w
  return(2 * (E.C + E.U))
}

approximate.ET2 <- function(X,v=one.version,u=all.common,m=4){
 J <- ncol(X)
 N <- nrow(X)
 one.version <- rep(1,N)
 all.common <- rep(F,J)
 V <- outer(v,v,"==")
 X. <- X %*% matrix(1,nrow=J)
 length.u <- length(u[u==F])

 ET2.same <- matrix(0,J+1,J+1)
 for(i in 0:J) for(j in 0:J) ET2.same[i+1,j+1] <- exp.T2.function(i,j,J,0,m)

 ET2.diff <- matrix(0,J+1,J+1)
 for(i in 0:J) for(j in 0:J) ET2.diff[i+1,j+1] <- exp.T2.function(i,j,J,J-length.u,m)

 ET2 <- matrix(0,N,N)
 for (i in 1:N) for(j in 1:N){
   ET2[i,j] <- ifelse(V[i,j], ET2.same[X.[i]+1,X.[j]+1],ET2.diff[X.[i]+1,X.[j]+1])
 }
 diag(ET2) <- 0
 return(ET2)
}

max.T2.function <- function(L,K,J,U){
  A <- min(U, K, J-L)
  B <- min(J-K,J-A-L)
  return(2 * A + B)
}

compute.MT2 <- function(X,v=one.version,u=all.common){
  J <- ncol(X)
  N <- nrow(X)
  one.version <- rep(1,N)
  all.common <- rep(F,J)
  V <- outer(v,v,"==")
  X. <- X %*% matrix(1,nrow=J)
  length.u <- length(u[u==F])

  MT2.diff <- matrix(0,J+1,J+1)
  for(i in 0:J) for(j in 0:J) MT2.diff[i+1,j+1] <- max.T2.function(i,j,J,J-length.u)

  MT2.same <- matrix(0,J+1,J+1)
  for(i in 0:J) for(j in 0:J) MT2.same[i+1,j+1] <- max.T2.function(i,j,J,0)

  MT2 <- matrix(0,N,N)
  for (i in 1:N) for(j in 1:N){
    MT2[i,j] <- ifelse(V[i,j], MT2.same[X.[i]+1,X.[j]+1],MT2.diff[X.[i]+1,X.[j]+1])
  }
  return(MT2)
}

is.one <- function(x,y){ifelse(x==1 & y==1,T,F)}

compute.G <- function(X,Y,v=one.version){
 N <- nrow(X)
 J <- ncol(X)
 one.version <- rep(1,N)
 # TVs = number of different test versions
 TVs <- length(unique(v))
 item.order <- list()
 Nv <- list()
 G <- matrix(0,N,N)
 for (i in 1:TVs){
   Nv[[i]] <- length(v[v==i])
   item.order[[i]] <- rev(order(rep(1/Nv[[i]],Nv[[i]]) %*% X[v==i,] + rnorm(J,0,1e-10)))
   X <- X[,item.order[[i]]]
   Y <- Y[,item.order[[i]]]
   # M = matching response options
   M <- array(0,c(N,N,J))
   for (j in 1:J){M[,,j] <- outer(Y[,j],Y[,j],"==")}
   # A = both item scores equal to 1
   A <- array(0,c(N,N,J))
   for (j in 1:J){A[,,j] <- outer(X[,j],X[,j],is.one)}
   # M = matching item scores equal to 1
   M <- M*A
   # Diagonal is tot number of Guttman errors. Off-diagonal elements are the
   # number of Guttman errors of row examinee minus the number of Guttman errors
   # of row examinee's obtained on the nonmatching items with column examinee.
   for (i in which(v==i)) for (j in 1:N){
     G[i,j] <- sum(X[i,]* M[i,j,] * cumsum(abs(X[i,]-1)))
   }
 }
 diag(G) <- 0
 return(G)
}

max.G.function <- function(X,v=one.version,u=all.common){
  N <- nrow(X)
  J <- ncol(X)
  all.common <- rep(F,J)
  one.version <- rep(1,N)
  TVs <- length(unique(v))
  item.order <- list()
  Nv <- list()
  max.G <- list()
  for (copier in 1:TVs){
     Nv[[copier]] <- length(v[v==copier])
     item.order[[copier]] <- rev(order(rep(1/Nv[[copier]],Nv[[copier]])
        %*% X[v==copier,] + rnorm(J,0,1e-10)))
     max.G[[copier]] <- list()
     for (source in 1:TVs){
        max.G[[copier]][[source]] <- rep(0,J+1)
        if (copier==source){
           u.o <- all.common
        }else{
           u.o <- u[item.order[[copier]]]
        }
        t. <- 1
        common.item <- length(u.o[u.o==F])
        unique.item <- length(u.o[u.o==T])
        X.t <- rep(0,J)
        repeat{
           t. <- t. + 1
           X.t.1 <- X.t.2 <- X.t
           dc <- ifelse(common.item > 0,which(!u.o)[common.item],NA)
           du <- ifelse(unique.item > 0,which(u.o)[unique.item],NA)
           X.t.1[dc] <- 1
           X.t.2[du] <- 1
           G.1 <- ifelse(!is.na(dc),sum((X.t.1 * !u.o) * cumsum(abs(X.t.1-1))),0)
           G.2 <- ifelse(!is.na(du),sum((X.t.2 * !u.o) * cumsum(abs(X.t.2-1))),0)
           if (G.1 >= G.2){
              max.G[[copier]][[source]][t.] <- G.1
              X.t <- X.t.1
              common.item <- common.item - 1
           } else{
              max.G[[copier]][[source]][t.] <- G.2
              X.t <- X.t.2
              unique.item <- unique.item - 1
           }
           if(t.==J+1)break
        }
     }
  }
  return(max.G)
}

compute.MG <- function(X,v=one.version,u=all.common){
  J <- ncol(X)
  N <- nrow(X)
  one.version <- rep(1,N)
  all.common <- rep(F,J)
  max.G <- max.G.function(X,v,u)
  sum <- X %*% rep(1,J)
  MG <- matrix(0,N,N)
  for(i in 1:N) for(j in 1:N){
     MG[i,j] <- max.G[[v[i]]][[v[j]]][sum[i]+1]
  }
  return(MG)
}