## Normal--normal model


## y[,,t] <- x[,,t](H[,,t]*Q[,,t]) + v[,,t]
## v[,,t] ~ N(z[,,t],R[,,t])

## This model is invariant wrt background variables Z.

FixedQNormalParam <-
setClass("FixedQNormalParam",
         slots = c(Hvec="vector",
                   R="matrix", z="vector",
                   is.Rpd ="logical",Rinv="matrix"),
         contains=("ParameterObject"))
FixedQNormalParam <- function (Hvec=numeric(),z=numeric(),R=diag(length(z)),
                               is.Rpd = is.pdm(R), Rinv=solve(R)) {
  if (missing(R) && !missing(Rinv)) {
    is.Rpd <- is.dpm(Rinv)
    if (is.Rpd) {
      R <- solve(Rinv)
    } else {
      R <- matrix(NA)
    }
  }
  if (missing(Rinv) && !is.Rpd) {
    Rinv <- matrix(NA)
  }
  new("FixedQNormalParam",Hvec,R,z,is.Rpd,Rinv)
}


setGeneric("getHvec",function(param) standardGeneric("getHvec"))
setGeneric("getH",function(param,Q) standardGeneric("getH"))
setGeneric("getH<-",function(param,Q,value) standardGeneric("getH<-"))
setGeneric("getR",function(param) standardGeneric("getR"))
setGeneric("getR<-",function(param,value) standardGeneric("getR<-"))
setGeneric("getz",function(param) standardGeneric("getz"))

setMethod("pvec","FixedQNormalParam",
          function(param) c(as.vector(param@Hvec),param@z,
                            diag(param@R),param@R[lower.tri(param@R)]))
setMethod("pvec<-","FixedQNormalParam",
          function(param,value) {
            h <- length(param@Hvec)
            p <- length(param@z)
            if (length(value) != h+p+p + p*(p-1)/2)
              stop("Parameter vector is incorrect length.")
            param@Hvec <- value[1:h]
            param@z <- value[h+1:p]
            R <- diag(value[h+p+(1:p)])
            R[lower.tri(R)] <- value[(h+2*p+1):length(value)]
            R <- t(R)
            R[lower.tri(R)] <- value[(h+2*p+1):length(value)]
            ## Calculates Rinv and pd for us.
            getR(param) <- R
            param
          })


setMethod("getHvec","FixedQNormalParam",
          function(param) param@Hvec)

setMethod("getH","FixedQNormalParam",
          function(param,Q) {
            result <- matrix(0,nrow(Q),ncol(Q))
            result[abs(Q)>0] <- param@Hvec
            result
          })
setMethod("getH<-","FixedQNormalParam",
          function(param,Q,value) {
            param@Hvec <- value[abs(Q)>0]
            param
          })



setMethod("getR","FixedQNormalParam",
          function(param) param@R)
setMethod("getR<-","FixedQNormalParam",
          function(param,value) {
            param@R <- value
            param@is.Rpd <- is.pdm(value)
            if (param@is.Rpd) {
              param@Rinv <- solve(value)
            } else {
              param@Rinv <- matrix(NA)
            }
            param
          })


setMethod("getz","FixedQNormalParam",
          function(param) param@z)


## Priors
## H  ~ N(hMean,hStd)
## R ~ Wishart(Rmean,Rdf)
## z ~ N(zMean,zStd)

FixedQNormalEM <-
  setClass("FixedQNormalEM",
           slots=c(hMean="vector",hStd="vector",
                   zMean="vector",zVar="matrix",ziVar="matrix",
                   Rscale="matrix",Rdf="numeric",
                   Riscale="matrix"),
           contains="EvidenceModel")

FixedQNormalEM <- function (Q=matrix(),param=NullParameter(),
                            hMean=rep(1,sum(Q)),hStd=rep(1,length(hMean)),
                            zMean=rep(0,nrow(Q)),zVar=diag(length(zMean)),
                            ziVar=sovle(zVar),
                            Rscale=diag(length(zMean)),
                            Riscale = solve(Rscale),
                            Rdf=length(zMean)+1) {
  Q <- array(as.logical(Q),dim(Q))
  d <- nrow(Q)                          #Residual dimensions
  hh <- sum(Q)                          #Number of non-zero parameters
  if (length(hMean) != hh || length(hStd) !=hh) {
    stop("Length of H mean or SD does not match number of non-zero elements in Q-matrix.")
  }
  if (length(zMean) != d) {
    stop("Mean vector does not match number of rows in Q-matrix.")
  }
  if (missing(zVar) && !missing(ziVar)) {
    zVar <- solve(ziVar)
  }
  if (nrow(zVar) !=d) {
    stop("Dimensions of variance of z does not match length of z.")
  }
  if (!is.pdm(zVar)) {
    stop("zVar is not positive definite.")
  }
  if (missing(Rscale) && !missing(Riscale)) {
    Rscale <- solve(Riscale)
  }
  if (nrow(Rscale) !=d) {
    stop("Dimensions of Rscale does not match length of z.")
  }
  if (!is.pdm(Rscale)) {
    stop("Rscale is not positive definite.")
  }
  if (Rdf <= d) {
    stop("Improper prior distribution for R, df less than matrix rank.")
  }
  if (!is(param,"FixedQNormalParam") && !is(param,"NullParameter")) {
      stop("Expected a FixedQNormal Parameter.")
  }

  new("FixedQNormalEM", param = param,
      paramtype="FixedQNormalParam",
      Q=Q,hMean=hmean(),hStd=hStd,
      zMean=zMean,zVar=zVar,ziVar=ziVar,
      Rscale=Rscale,Riscale=Riscale,Rdf=Rdf)
}


setMethod("lpriorEMParam", "FixedQNormalEM",
          function(em,param=parameters(em),background=NA,QQ=Q(em)) {
            sum(dnorm(getHvec(param),em@hMean,em@hStd,log=TRUE)) +
            dmvnorm(getz(param),em@zMean,em@zVar,log=TRUE) +
              log(dwish(getR(param),em@Rdf,em@Rscale))
          })


setMethod("drawEMParam","FixedQNormalEM",
function(em,background=NA,QQ=Q(em)) {
  Hvec <- rnorm(sum(abs(QQ)>0),em@hMean,em@hStd)
  z <- rmvnorm(1,em@zMean,em@zVar)
  Rmat <- riwish(em@Rdf,em@Rscale)
  FixedQNormalParam(Hvec=Hvec,R=Rmat,z=z)
})


setMethod("llikeObs", "FixedQNormalEM",
          function(em,obs,latent,background=NA,QQ=Q(em),param=parameters(em)) {
            Y <- obs - latent%*%t(getH(param,QQ)*scaleQ(QQ))
            dmvnorm(Y, getz(param),getR(param),log=TRUE)
 })

setMethod("drawObs", "FixedQNormalEM",
          function(em,latent,background=NA,QQ=Q(em),param=parameters(em)) {
            ncases <- nrow(latent)
            v <- rmvnorm(ncases,getz(param),getR(param))
            latent%*%t(getH(param,QQ)*scaleQ(QQ)) + v
          })


setMethod("optimalEMParams", "FixedQNormalEM",
          function (em,obs,latent,weights=1,
                    background=NA,QQ=Q(em),param=parameters(em),
                    control=list(),maxit=100,nowarn=FALSE,
                    Bayes=TRUE) {



            })