#######################################################
##    Bayesian Compound Poisson GLMM                 ##
## Author: Wayne Zhang, actuary_zhang@hotmail.com    ##
#######################################################

bcpglmm <- function(formula, link = "log", data, inits = NULL,
                   weights, offset, subset, na.action, contrasts = NULL, 
                   n.chains=3, n.iter=2000, n.burnin=floor(n.iter/2),
                   n.thin=max(1, floor(n.chains * (n.iter - n.burnin) / n.sims)),
                   n.sims=1000, n.report=1000, prior.beta.mean=NULL, prior.beta.var=NULL, 
                   bound.phi=100, bound.p=c(1.01,1.99),  prior.Sigma = NULL,
                   tune.iter=4000, n.tune=10, tune.weight=0.25,
                   basisGenerators = c("tp","tpU","bsp","sp2d"),
                   method="dtweedie",...) {
      
  call <- amer:::expand.call()  
  if (missing(data)) 
    data <- environment(formula)   
  link.power <- make.link.power(link)
  # identify smooth terms 
  tf <- terms.formula(formula, specials = eval(call$basisGenerators, 
        parent.frame(2)))
  n.f <- length(unlist(attr(tf, "specials")))
  # create model frame and get factor list
  if (n.f) {
    call2 <- as.list(call)[-1]
    m <- match(c("formula", "data", "weights", "offset", 
                        "contrasts", "basisGenerators"), names(call2),0L)
    call2 <- call2[m]
    #setup <- do.call(amer:::amerSetup, as.list(call2))     
    setup <- do.call(frFL, as.list(call2))
    fr <- setup$m$fr 
    FL <- setup$m$FL
  } else {  
    fr <- lme4:::lmerFrames(call, formula, contrasts)
    FL <- lme4:::lmerFactorList(formula, fr, 0L, 0L)
  }
  
  dm <- lme4:::mkZt(FL, NULL)  
  n.obs <- unname(dm$dd['n'])
  n.beta <- unname(dm$dd['p'])
 
  # check initial values
  if (!is.null(inits))
    check.inits.bcpglmm(inits=inits,n.beta=ncol(fr$X),
                       n.term=as.numeric(dm$dd['nt']), n.chains=n.chains)
  M1 <- length(levels(dm$flist[[1]]))
  if (M1 >= n.obs) {
        msg1 <- "Number of levels of a grouping factor for the random effects\n"
        msg3 <- "n, the number of observations"
        if (dm$dd["useSc"]) 
            stop(msg1, "must be less than ", msg3)
        else if (M1 == n.obs) 
            message(msg1, "is *equal* to ", msg3)
    }
        
  # default prior info
  if (is.null(prior.beta.mean))
    prior.beta.mean <- rep(0, n.beta)
  if (is.null(prior.beta.var))
    prior.beta.var <- rep(10000, n.beta)
  if (is.null(prior.Sigma))
    prior.Sigma <- prior.Sigma.default(dm$ST)
  
  # default weights and offsets if NULL    
 # default offset and prior wts
 if (is.null(fr$wts) || length(fr$wts)==0)  
    wts <- as.double(rep(1,n.obs)) else 
    wts <- fr$wts 
  if (is.null(fr$off) || length(fr$off)==0)
    off <- as.double(rep(0,n.obs)) else 
    off <- fr$off
  
  # dimensions used in simulation
  nc <- unlist(lapply(dm$ST, ncol))
  nlev <- diff(dm$Gp)/nc  
  n.keep <- floor((n.iter-n.burnin) / n.thin)
  n.sims <- n.chains * n.keep    
  dims <- list(n.obs= n.obs,
           n.beta=as.integer(unname(n.beta)),           
           n.pos= as.integer(sum(fr$Y>0)),
           n.term = as.integer(dm$dd['nt']),
           n.u = as.integer(dm$dd['q']),
           n.all = as.integer(dm$dd['p'] + dm$dd['q'] + sum(nc^2)+2),
           n.chains=as.integer(n.chains), 
           n.iter=as.integer(n.iter), 
           n.burnin=as.integer(n.burnin),
           n.thin=as.integer(n.thin), 
           n.keep=as.integer(n.keep),
           n.sims=as.integer(n.sims),
           n.report=as.integer(n.report),
           tune.iter=as.integer(tune.iter),
           n.tune=as.integer(n.tune))
  dims <- unlist(dims)
  
  C <- 2.38 * 2.38
  # proposal covariance matrix              
  ebeta.var <- diag(1, nrow=n.beta) * C /n.beta
  eu.var <- diag(1, nrow=dm$dd['q']) * C/dm$dd['q']
  ephi.var <- ep.var <-  C 
                   
  # generate initial values if necessary
  if (is.null(inits)) {
    pstart <- 1.5
    fit.start <- glm.fit(fr$X, fr$Y, weights = wts, offset = off, 
        family = tweedie(var.power=pstart,link.power=link.power), 
        intercept = attr(attr(fr$mf, "terms"), "intercept") > 0)
    betastart <- as.numeric(fit.start$coefficients)
    ustart <- rnorm(dm$dd[["q"]])
    phistart <- sum((fr$Y-fit.start$fitted.values)^2/fit.start$fitted.values^pstart)/
                    fit.start$df.residual
		Sigmastart <- lapply(dm$ST,function(x) x%*%t(x))
  	Sigmastartv <- unlist(lapply(Sigmastart, as.numeric))
    inits.start <- c(betastart, phistart, pstart, ustart, Sigmastartv)
    inits <- vector("list",n.chains)
	  inits[[1]] <- c(betastart, phistart, pstart, ustart, Sigmastartv)
	  if (n.chains>1){
		    for (i in 2:n.chains)
			    inits[[i]] <- c(as.numeric(betastart + rnorm(n.beta,0,0.5)),
						        runif(1,phistart/2,1.5*phistart),
						        runif(1,(bound.p[1]+pstart)/2,(bound.p[2]+pstart)/2),
                    ustart+rnorm(length(ustart),0,0.5) ,
                    Sigmastartv)
	   }
    #update proposal
    #ebeta.var <- as.matrix(vcov(fit.start)*C/n.beta)
    eu.var <- lapply(1:dm$dd['nt'], function(x) 
              kronecker(Sigmastart[[x]],diag(1,nrow=nlev[x])))
    eu.var <- as.matrix(do.call(bdiag, eu.var))*C/dm$dd['q']
  } else{
    betastart <- inits[[1]]$beta
    phistart <- inits[[1]]$phi
    pstart <- inits[[1]]$p
    ustart <- inits[[1]]$u
    Sigmastart <- inits[[1]]$Sigma
    inits <- lapply(inits, function(x) c(x$beta, x$phi, x$p, x$u,
                                      unlist(lapply(x$Sigma, as.numeric))))
  }

  # run MCMC   
    # input for the C function 	     
    input <- list(X=fr$X,
               y=as.double(fr$Y),
               Zt = dm$Zt, 
               ygt0= as.integer(which(fr$Y>0L)-1),
               offset=as.double(off),
               pWt=as.double(wts),
               mu = double(dims["n.obs"]),
               eta = double(dims["n.obs"]),
               inits = inits,
               beta=as.double(betastart),
               u= as.double(ustart),
               phi=as.double(phistart),
               p=as.double(pstart),
               link.power=as.double(link.power),
               pbeta.mean=as.double(prior.beta.mean),
               pbeta.var=as.double(prior.beta.var),
               pSigma = prior.Sigma,
               bound.p=as.double(bound.p),
               bound.phi=as.double(bound.phi),    
               ebeta.var=as.double(ebeta.var),
               eu.var = as.double(eu.var),    
               ep.var= as.double(ep.var),
               ephi.var = as.double(ephi.var),               
               dims=dims,
               tune.weight=as.double(tune.weight),
               Gp = unname(dm$Gp),
               Sigma = Sigmastart,
               ncol = as.integer(nc), 
               nlev = as.integer(nlev),
               lambda = as.double(2),
               k = as.integer(1),
               simT = as.integer(rep(1,dims["n.pos"])))

  if (method=="dtweedie")
    sims.list<- .Call("bcpglmm_gibbs_tw",input) 
  if (method=="latent")
    sims.list<- .Call("bcpglmm_gibbs_lat",input)
  xnames <- names(fr$fixef)
  unames <- paste("u", 1:dm$dd['q'],sep="")
  snames <- sapply(1:length(dm$ST), function(x){
              tm <- apply(expand.grid(1:nc[x],1:nc[x]),1,
                          paste, collapse=",")
              paste("Sigma",x,"[",tm,"]",sep="")
  })
  snames <- unlist(snames)
  sims.list <- lapply(sims.list, function(x){ 
                  dimnames(x) <- list(NULL, 
                      c(xnames,"phi","p",unames,snames))
                  return(x)}) 
  # coerce to mcmc object                  
  sims <- lapply(sims.list, as.mcmc)
  sims <- as.mcmc.list(sims)
 
  # coerce sims.list to mcmc.list from coda
  ans <- new("bcpglmm", 
             n.chains=as.integer(n.chains), 
             n.iter=as.integer(n.iter), 
             n.burnin=as.integer(n.burnin),
             n.thin=as.integer(n.thin), 
             n.sims=as.integer(dims["n.sims"]), 
             sims.list=sims,
             link.power=link.power,
             call=call,
             formula=formula,
             model.frame = fr$mf,
             contrasts=contrasts,
             inits = inits,
             prop.var = list(beta.var = matrix(input$ebeta.var,n.beta,n.beta), 
                             u.var = matrix(input$eu.var,dm$dd['q'],dm$dd['q']),
                             phi.var = input$ephi.var,
                             p.var = input$ep.var))  
  return(ans)
                  
}