#####
## Generic optimizers

#######################################################################
## EM parameter optimization.  The default uses optim, which computes analytic derivatives
## and hence might be slow.

### The generic methods use unconstratined optimization, but the
### problem is that we generally need CONSTRAINED optimization.  The
### solution is to have the lprior and llike methods return -Inf if the
### conditions are not met.  For example, if the variance is not
### positve or the covariance matrix is not positive definite, it
### should return -Inf.  As long as optim starts at a legal value,
### this should be fine.

setGeneric("optimalPMParams",
           function (pm,latent,weights=1,
                    background=NA,param=parameters(pm),
                    control=list(),maxit=100,nowarn=FALSE,
                    Bayes=TRUE)
             standardGeneric("optimalPMParams"))

setMethod("optimalPMParams", "ProficiencyModel",
          function (pm,latent,weights=1,
                    background=NA,param=parameters(pm),
                    control=list(),maxit=100,nowarn=FALSE,
                    Bayes=TRUE) {
            stop("No parameter optimizer for class ",class(pm))
          })

genericPMoptimizer <-
  function (pm,latent,weights=1,
            background=NA,param=parameters(pm),
            control=list(),maxit=100,nowarn=FALSE,
            Bayes=TRUE) {
    if (is(param,"NullParameter")) {
      ## No initialization, randomly initialize.
      param <- drawPMParam(pm)
    }

    control$maxit <- maxit
    pv0 <- pvec(param)
    if (Bayes) {
      llike <- function (pv) {
        pvec(param) <- pv
        ll <- lpriorLatent(pm,latent,background,param)
        -2*(sum(weights*ll)+lpriorPMParam(pm,param,background))
      }
    } else {
      llike <- function (pv) {
        pvec(param) <- pv
        ll <- lpriorLatent(pm,latent,background,param)
        -2*sum(weights*ll)
      }
    }
    optp <- optim(pv0,llike,control)
    if (!nowarn && optp$convergence !=0) {
      warning("Optimization did not coverge.")
      warning(optp$message)
    }
    pvec(param) <- optp$par
    list(param=param,deviance=optp$value,
         convergence=(optp$convergence==0),optout=optp)
  }



setMethod("optimalEMParams", "EvidenceModel",
          function (em,obs,latent,weights=1,
                    background=NA,QQ=Q(em),param=parameters(em),
                    control=list(),maxit=100,nowarn=FALSE,
                    Bayes=TRUE) {
            stop("No parameter optimizer for class ",class(em))
          })

genericEMoptimizer <-
  function (em,obs,latent,weights=1,
            background=NA,QQ=Q(em),param=parameters(em),
            control=list(),maxit=100,nowarn=FALSE,
            Bayes=TRUE) {

    control$maxit <- maxit
    pv0 <- pvec(param)
    if (Bayes) {
      llike <- function (pv) {
        pvec(param) <- pv
        ll <- llikeObs(em,obs,latent,background,QQ,param)
        -2*(sum(weights*ll)+lpriorEMParam(em,param,background,QQ))
      }
    } else {
      llike <- function (pv) {
        pvec(param) <- pv
        ll <- llikeObs(em,obs,latent,background,QQ,param)
        -2*sum(weights*ll)
      }
    }
    optp <- optim(pv0,llike,control)
    if (!nowarn && optp$convergence !=0) {
      warning("Optimization did not coverge.")
      warning(optp$message)
    }
    pvec(param) <- optp$par
    list(param=param,deviance=optp$value,
         convergence=(optp$convergence==0),optout=optp)
  }