### This is a series of files for generating an AMD from a Q-Matrix.
### Version one assumes that there is exactly one observable per task,
### but this is a fairly common case in retrofit designs.

### Currently assume that the layout has the following fields:
## TaskModel, EvidenceModel, TaskName, ObsName, FormCode, CPTType, Difficulty
## and that all of the columns following difficulty represent skills.

## Reads a Q matrix from a CSV file.
## Most of the arguments are passed to read.csv which does most of the work.
## Massages the first couple of columns to be character instead of factors.
"read.Q" <- function (filename, ...) {
  Qmat <- read.csv(filename,...)
  Qmat$TaskModel <- as.character(Qmat$TaskModel)
  Qmat$EvidenceModel <- as.character(Qmat$EvidenceModel)
  Qmat$TaskName <- as.character(Qmat$TaskName)
  Qmat$ObsName <- as.character(Qmat$ObsName)
  Qmat$CPTType <- as.character(Qmat$CPTType)
  if (length(grep("Distribution", Qmat$CPTType)) == 0) {
    Qmat$CPTType <- paste(Qmat$CPTType, "Distribution")
  }
  ## Should be able to leave Anchor as a factor.
  ##Qmat$Anchor <- as.character(Qmat$Anchor)
  Qmat
}

## Extracts names (and number) of skill columns
"QSkillNames" <- function (Qmat) {
  start <- match("Difficulty",names(Qmat)) +1
  if (is.na(start)) {
    stop ("Difficulty column not found, is argument a Q-Matrix?")
  }
  names(Qmat)[start:length(names(Qmat))]
}

"QSkillNames<-" <-
function (Qmat, value) 
{
    start <- match("Difficulty", names(Qmat)) + 1
    if (is.na(start)) {
        stop("Difficulty column not found, is argument a Q-Matrix?")
    }
    names(Qmat)[start:length(names(Qmat))] <- value
    Qmat
}

## Build a collection of proficiency variables from the names in the Q-Matrix.
"QProficiencyVars" <- function (names, levels) {
  sapply(names,function(name)
         varDescription(name=name,role="SM",levels=levels),
         simplify=FALSE)
}

## Find the first task ID associated with each evidence model ID.
## Used later to select exemplar rows when building code for evidence model.
"QExemplarTasks" <- function (Qmat) {
  ems <- unique(Qmat$EvidenceModel)
  tasks <- sapply(ems,function (em) {
    Qmat$TaskName[Qmat$EvidenceModel==em][1]
    })
  tasks
}

"QbuildEMList" <- function (Qmat, caf, proficiencies, observables,
                            emlist = unique(Qmat$EvidenceModel), 
                            stds = list(diff = 2.0,disc=.5, dinc=.4)) {
  ## Assumes 1 TM per EM
  lapply(emlist, function (emname) {
    tmname <- Qmat$TaskModel[Qmat$EvidenceModel==emname][1]
    QbuildEM(Qmat,emname,tmname,caf,proficiencies, observables,stds=stds)
  })

}

"QbuildEM" <- function(Qmat,emname, tmname, caf, proficiencies, observables,
                       taskname=Qmat$TaskName[Qmat$EvidenceModel==emname
                       & Qmat$TaskModel == tmname][1],
                       stds = list(diff = 2.0,disc=.5, dinc=.4)) {
  linkid <- paste(caf,emname,"",sep="/") # Force trailing /
  tid <- paste(caf,tmname,"",sep="/")
  varnames <- character(0)
  rows <- Qmat[Qmat$EvidenceModel == emname & Qmat$TaskModel == tmname
               & Qmat$TaskName == taskname, ]
  dists <- vector(nrow(rows),mode="list")
  for (irow in 1:nrow(rows)) {
    dist <- buildDistFromRow(rows[irow,],proficiencies,
                                    observables, NULL, stds)
    cond <- dist$conditions
    cons <- dist$consequences
    varnames <- union( union(varnames, cons), cond)

    dists[[irow]] <-dist
  }
  vars <- c(observables,proficiencies)[varnames]
  notfound <- sapply(vars, is.null)
  if (any(notfound)) {
    stop("Couldn't find variable(s) ", varnames[notfound])
  }
  gmModel(linkid,tid,type="EM",variables=vars,distributions=dists)
}

"buildDistFromRow" <- function (row, proficiencies, observables,
                                dinc = NULL,
                                stds = list(diff = 2.0,disc=.5, dinc=.4)) {
  consequences <- row$ObsName
  # This next assumes there is exactly one consequence, currently only
  # case supported by StatShop.  
  obs <- observables[[consequences]]
  type <- as.character(row$CPTType)
  diff <- row$Difficulty
  pnames <- intersect(names(row), names(proficiencies))
  Q <- row[, pnames]
  conditions <- pnames[Q > 0]
  disc <- log(Q[Q>0])
  params <- buildParams(diff,conditions,disc,obs,dinc)
  cov <- buildCov(params,conditions,obs,stds)
  parameterSet <- parameterVector(params,covariance = cov)
  distribution(conditions, consequences, type, parameterSet=parameterSet)
}

buildParams <- function(diff, conditions, disc, obs, dinc=NULL) {
  names(diff) <- "Task Difficulty"
  names(disc) <-paste("Relative Importance(",conditions,")",sep="")
  params <- c(diff,disc)
  incs <- obs$levels
  nlevel <- length(incs)
  if (nlevel > 2) {
    if (is.null(dinc)) {
      dinc <- qnorm( (1:(nlevel-1)-.5)/(nlevel-1))
      dinc <- rev(diff(dinc))
    }
    incs <- incs[1:(length(incs)-2)]
    ## User could have specified dinc as single value or vector, this
    ## forces it to the right length
    dvals <- rep(0,length=length(incs))
    dvals <- dvals + dinc
    names(dvals) <- paste("Level Difficulty Increment(",obs$name,
                          ",",incs,")",sep="")
    params <- c(params,dvals)
  }
  params
}

"buildCov" <-
  function (params, conditions, obs, stds = list(diff = 2.0,disc=.5, dinc=.4)) { 
  cov <- matrix(0,length(params),length(params),
                dimnames=list(names(params),names(params)))
  diag(cov)[1] <- stds$diff*stds$diff
  npar <- length(conditions)
  if (npar > 0)
    diag(cov)[2:(1+npar)] <- stds$disc*stds$disc
  nlevels <- length(obs$levels)
  if (nlevels > 2)
    diag(cov)[(2+npar):(2+npar+nlevels-2-1)] <- stds$dinc*stds$dinc
  cov
}

"QbuildAMD" <-
function (Qmat, sm, observables, stds = list(diff = 2, disc = 0.5, 
    dinc = 0.4)) 
{
    caf <- getECDIDCaf(sm$id)
    smname <- getECDIDModel(sm$id)
    smlist <- list(sm)
    names(smlist) <- smname
    emlist <- QbuildEMList(Qmat, caf, sm$variables, observables, 
        stds = stds)
    amd(smlist, emlist, list(), filename = NULL)
}