## This model has a normal initial status and normal growth with a
## time dependent variance, hence Brownian motion.

## x[,,t] <- x[,,t-1]F[,,t] + u[,,t]B[,,t]+w[,,t]
## w[,,t] ~ N(0,Sigma[,,t])

## x==truth/est, u==act

## x[,,0] ~ N(mu,A)

## This model is invariant across the background data background.

TimeInvariantNormalParam <-
setClass("TimeInvariantNormalParam",
         representation(Fmat="array",Bmat="array",
                        Sigma="array",mu="vector",Amat="array"))

## Priors
## F  ~ N(Fmean,Fstd)
## B  ~ N(Bmean,Bstd)
## mu ~ N(muMean,muStd)
## Sigma ~ Wishart(Smean,Sdf)
## A ~ Wishart(Amean,Adf)

TimeInvariantNormalPGM <-
  setClass("TimeInvariantNormalPGM",
  representation(fMean="matrix",fStd="matrix",
                 bMean="matrix",bStd="matrix",
                 Smean="matrix",Sdf="numeric",
                 muMean="vector",muStd="matrix",
                 Amean="matrix",Adf="numeric"
                 ),
  prototype(fMean=matrix(c(.6,.3,.3,.6),2,2), fStd=matrix(.25,2,2),
            bMean=matrix(c(.1,.1),1,2), bStd=matrix(.05,1,2),
            Smean=matrix(c(.7,.3,.3,.7),2,2), Sdf=3,
            muMean=c(Mechanics=2,Fluency=2),  muStd=diag(c(Mechanics=1,Fluency=1)),
            Amean=matrix(c(.6,.4,.4,.6),2,2), Adf=3,
            usesSuccess=FALSE),
  contains="ProficiencyGrowthModel")

setMethod("drawPMParam","TimeInvariantNormalPGM",
function(pm,background) {
  Fmat <- matrix(rnorm(length(pm@fMean)),dim(pm@fMean))
  Fmat <- Fmat * pm@fStd + pm@fMean

  Bmat <- matrix(rnorm(length(pm@bMean)),dim(pm@bMean))
  Bmat <- Bmat*pm@bStd + pm@bMean

  mu <- as.vector(rmvnorm(1,pm@muMean,pm@muStd))

  Sigma <- riwish(pm@Sdf,pm@Smean)
  Amat <- riwish(pm@Adf,pm@Amean)
  new("TimeInvariantNormalParam",
      Fmat=Fmat, Bmat=Bmat, Sigma=Sigma, mu=mu, Amat=Amat )
})

setMethod("nlatent","TimeInvariantNormalPGM",
          function(pm) {
            length(pm@muMean)
          })

setMethod("naction","TimeInvariantNormalPGM",
          function(pm) {
            nrow(pm@bMean)
          })


setMethod("drawInitialLatent", "TimeInvariantNormalPGM",
          function(pm,ncases=1,background=NA,param=parameters(pm)) {
            latent <- rmvnorm(ncases,param@mu,param@Amat)
          })


setMethod("lpriorLatent", "TimeInvariantNormalPGM",
          function(pm,latent,background=NA,param=parameters(pm)) {
            sum(dmvnorm(latent,param@mu,param@Amat,log=TRUE))
          })



setMethod("advanceLatent",
          c("TimeInvariantNormalPGM","TimeInvariantNormalParam"),
function (pm,param,latent,success=NA,action=NA,background=NA,time=1) {
  ncases <- dim(latent)[1]
  nvar <- nlatent(pm)
  if (is.null(dim(action))) {
    action <- matrix(action,ncases,length(action),byrow=TRUE)
  }
  if (time != round(time)) {
    stop("Only integral time values supported.")
  } else {
    Fmat <- param@Fmat
    i <- 1
    while (i<time) {
      Fmat <- Fmat %*% param@Fmat
      i <- i+1
    }
  }
  w <- rmvnorm(ncases,rep(0,nvar),param@Sigma)
  latent%*%Fmat + action%*%param@Bmat + w
})

setMethod("llikeAdvanceLatent",
          c("TimeInvariantNormalPGM","TimeInvariantNormalParam"),
function (pm,param,oldlatent,newlatent,success=NA,action=NA,background=NA,time=1) {
  ncases <- dim(newlatent)[1]
  nvar <- nlatent(pm)
  if (is.null(dim(action))) {
    action <- matrix(action,ncases,length(action),byrow=TRUE)
  }
  if (time != round(time)) {
    stop("Only integral time values supported.")
  } else {
    Fmat <- param@Fmat
    i <- 1
    while (i<time) {
      Fmat <- Fmat %*% param@Fmat
      i <- i+1
    }
  }
  dmvnorm(newlatent - oldlatent%*%Fmat + action%*%param@Bmat,
              rep(0,nvar),
              param@Sigma, log=TRUE)
})