######################################################################
## R functions implementing the eCADS normalization method for two- ##
## channel microarrays, as presented in:                            ##
##                                                                  ##
##   Dabney A.R. and Storey J.D. (2007) Normalization of two-       ##
##   channel microarrays accounting for experimental design and     ##
##   intensity-dependent relationships.  Genome Biology, 8:R44.     ##
##                                                                  ##
## See the manual at:                                               ##
##                                                                  ##
##   http://www.stat.tamu.edu/~adabney/ecads                        ##
##                                                                  ##
## for instructions on using the software, and for examples.        ##
##                                                                  ##
## Software written by Alan Dabney.                                 ##
######################################################################

ecadsArgs = function(Y, Z, varNames, varDims) {
  m = nrow(Y)
  n = ncol(Y) / 2
  p = varDims[1]
  q = length(varNames)

  ZZ = vector("list", q)
  ZZ[[1]] = Z[, 1:p]
  if(q > 1) 
    for(i in 2:q)
      ZZ[[i]] = Z[, 1:varDims[i] + sum(varDims[1:(i - 1)])]

  XX = rep(1, 2 * n)
  for(i in 1:q)
    for(j in 1:(varDims[i] - 1))
      XX = cbind(XX, ZZ[[i]][, j] - ZZ[[i]][, varDims[i]])

  bb = t(solve(t(XX) %*% XX) %*% t(XX) %*% t(Y))

  out = cbind(bb[, 1] + bb[, 2:p], bb[, 1] - drop(bb[, 2:p, drop = F] %*% rep(1, p - 1)))
  return(out)
}

ecadsReshapeDesignMatrix = function(Z, varDims) {
  p = varDims[1]
  q = length(varDims)

  Z0 = Z[, 1:p]
  for(i in 1:(q - 1)) {
    Z0 = cbind(Z0, Z[, 1:(varDims[i + 1] - 1) + sum(varDims[1:i])] - Z[, varDims[i + 1] + sum(varDims[1:i])])
  }

  return(Z0)
}

ecadsFit = function(Y, Z, Args, varNames, varDims, df = 5) {
  require(splines)
  m = nrow(Y)
  n = ncol(Y) / 2
  p = varDims[1]
  q = length(varNames)
  d = df * (sum(varDims[-1]) - (q - 1))

  ## reshape the design matrix
  Z0 = ecadsReshapeDesignMatrix(Z, varDims)

  ## basis matrices
  B = ns(as.numeric(Args), df = df, intercept = T)
  BB = vector("list", p)
  for(i in 1:p)
    BB[[i]] = predict(B, Args[, i])

  Gx = 0
  for(i in 1:p)
    Gx = Gx + kronecker(Z[, i], Args[, i])

  ## compute (XtX)^(-1)
  XX = matrix(0, nrow = d, ncol = d)
  for(i in 1:(2 * n)) {
    XB = 0
    for(j in 1:p)
      XB = XB + Z0[i, j] * BB[[j]]
    XB = kronecker(Z0[i, -(1:p), drop = F], XB)
    XX = XX + t(XB) %*% XB
  }
  XX = solve(XX)

  ## estimated parameters
  theta = rep(0, d)
  ii = 1:m
  for(i in 1:(2 * n)) {
    XB = 0
    for(j in 1:p)
      XB = XB + Z0[i, j] * BB[[j]]
    XB = kronecker(Z0[i, -(1:p), drop = F], XB)

    theta = theta + drop(XX %*% t(XB) %*% (as.numeric(Y) - Gx)[ii])
    ii = ii + m
  }

  beta = vector("list", q - 1)
  beta[[1]] = matrix(theta[1:(df * (varDims[2] - 1))], ncol = varDims[2] - 1)
  beta[[1]] = cbind(beta[[1]], -drop(beta[[1]] %*% rep(1, varDims[2] - 1)))
  for(i in 2:(q - 1)) {
    beta[[i]] = matrix(theta[1:(df * (varDims[i + 1] - 1)) + df * (sum(varDims[2:i]) - (i - 1))], 
      ncol = varDims[i + 1] - 1)
    beta[[i]] = cbind(beta[[i]], -drop(beta[[i]] %*% rep(1, varDims[i + 1] - 1)))
  }
  names(beta) = varNames[-1]

  ## normalized data 
  yTilde = Y
  for(i in 1:(2 * n)) {
    XB = 0
    for(j in 1:p)
      XB = XB + Z0[i, j] * BB[[j]]
    XB = kronecker(Z0[i, -(1:p), drop = F], XB)    

    yTilde[, i] = yTilde[, i] - drop(XB %*% theta)
  }

  ## output 
  out = list("yTilde" = drop(yTilde), "B" = B, "beta" = beta) 
 
  return(out) 
}