.packageName <- "bayesSurv"
###########################################
#### AUTHOR:     Arnost Komarek        ####
####             (2004)                ####
####                                   ####
#### FILE:       bayesDensity.R        ####
####                                   ####
#### FUNCTIONS:  bayesDensity          ####
####             print.bayesDensity    ####
####             plot.bayesDensity     ####
###########################################

### ======================================
### bayesDensity
### ======================================
bayesDensity <- function(dir = getwd(),
                         stgrid,
                         grid,
                         n.grid = 100,
                         skip = 0,                         
                         standard = TRUE,
                         unstandard = TRUE)
{
  thispackage = "bayesSurv"
  
  ## Check whether needed files are available
  ## * further,  whether at least first row has correct number of elements
  ## * and determine the MC sample size
  ## ==============================================================
  filesindir <- dir(dir)    ## character vector with available files   
  if (sum(!is.na(match(filesindir, "mweight.sim")))){
    mix <- read.table(paste(dir, "/mweight.sim", sep = ""), nrows = 1)
    k.max <- length(mix)
  }
  else
    stop("File with simulated values of mixture weights not found.")

  if (sum(!is.na(match(filesindir, "mmean.sim")))){
    mix <- read.table(paste(dir, "/mmean.sim", sep = ""), nrows = 1)
    kmax2 <- length(mix)
    if (k.max != kmax2) stop("Different kmax indicated by files mweight.sim and mmean.sim.")
  }
  else
    stop("File with simulated values of mixture means not found.")
   
  if (sum(!is.na(match(filesindir, "mvariance.sim")))){
    mix <- read.table(paste(dir, "/mvariance.sim", sep = ""), nrows = 1)     
    kmax2 <- length(mix)
    if (k.max != kmax2) stop("Different kmax indicated by files mweight.sim and mvariance.sim.")
  }
  else
    stop("File with simulated values of mixture variances not found.")

  if (sum(!is.na(match(filesindir, "mixmoment.sim")))){
    mix <- read.table(paste(dir, "/mixmoment.sim", sep = ""), header = TRUE)
    M <- dim(mix)[1]    
  }
  else
    stop("File mixmoment.sim not found.")     

  k.cond <- 1:k.max
    
  if (missing(skip)) skip <- 0
  else{
    if (skip > M) stop("You ask to skip more rows from the file than available.")
    if (skip < 0) skip <- 0
  }


  ## Try to guess the grid (from first at most 20 mixtures) if not given by the user
  if (!unstandard) grid <- 0
  if (!standard) stgrid <- 0
  if (missing(grid)){
    mus <- scan(paste(dir, "/mmean.sim", sep = ""), nlines = 20, skip = 1)
    sigma2s <- scan(paste(dir, "/mvariance.sim", sep = ""), nlines = 20, skip = 1)    
    mu.min <- min(mus)
    mu.max <- max(mus)
    sd.min <- sqrt(min(sigma2s))
    sd.max <- sqrt(max(sigma2s))
    sd.mean <- sqrt(mean(sigma2s))
    grid <- seq(mu.min - 2*sd.mean, mu.max + 2*sd.mean, length = n.grid)
  }
  if (missing(stgrid)){
    stgrid <- seq(-3, 3, length = n.grid)
  }
  ngrid <- length(grid)
  nstgrid <- length(stgrid)
  type <- 0*(standard & unstandard) + 1*(unstandard & (!standard) + 2*(!unstandard & standard))

  mcmc <- .C("bayesDensity", aver = double(ngrid * (1 + k.max)),  staver = double(nstgrid * (1 + k.max)),
                             intercept = double(M),               scale = double(M),
                             M.k = integer(1 + k.max),            as.character(dir),
                             as.double(grid),                     as.double(stgrid),
                             as.integer(k.max),                   as.integer(M),
                             as.integer(skip),
                             as.integer(ngrid),                   as.integer(nstgrid),
                             as.integer(type),                    err = integer(1),
             PACKAGE = thispackage)

  if (mcmc$err) stop("No results produced, something is wrong.")

  if (unstandard){
    mcmc$aver <- matrix(mcmc$aver, nrow = ngrid)
    mcmc$aver <- cbind(grid, mcmc$aver)
    mcmc$aver <- as.data.frame(mcmc$aver)
    rownames(mcmc$aver) <- paste(1:ngrid)
    colnames(mcmc$aver) <- c("grid", "unconditional", paste("k = ", 1:k.max, sep = ""))    
  }
  if (standard){
    mcmc$staver <- matrix(mcmc$staver, nrow = nstgrid)
    mcmc$staver <- cbind(stgrid, mcmc$staver)
    mcmc$staver <- as.data.frame(mcmc$staver)
    rownames(mcmc$staver) <- paste(1:nstgrid)
    colnames(mcmc$staver) <- c("grid", "unconditional", paste("k = ", 1:k.max, sep = ""))        
  }
  names(mcmc$M.k) <- c("unconditional", paste(1:k.max))
  names(mcmc$intercept) <- paste(1:M)
  names(mcmc$scale) <- paste(1:M)    
    
  charact <- data.frame(intercept = mcmc$intercept, scale = mcmc$scale)
  rownames(charact) <- paste(1:M)          

  density <- list()
  if (standard) density$standard <- mcmc$staver
  else          density$standard <- "Not asked."

  if (unstandard) density$unstandard <- mcmc$aver
  else            density$unstandard <- "Not asked."

  attr(density, "sample.size") <- mcmc$M.k
  attr(density, "moments") <- charact
  attr(density, "k") <- data.frame(k = mix[,1])
  
  class(density) <- "bayesDensity"
  return(density)   
}  


### ======================================
### print.bayesDensity
### ======================================
print.bayesDensity <- function(x, ...)
{
  cat("\nStandardized McMC average: \n")
  print(x$standard)  

  cat("\nUnstandardized McMC average: \n")
  print(x$unstandard)  
  
  return(invisible(x))
}  


### ======================================
### plot.bayesDensity
### ======================================
plot.bayesDensity <- function(x,
                              k.cond,
                              dim.plot = TRUE,
                              over = TRUE,
                              alegend = TRUE,
                              standard = TRUE,
                              type = "l",
                              bty = "n",
                              xlab = expression(epsilon),
                              ylab = expression(f(epsilon)),
                              lty,
                              xlim,
                              ylim,
                              xleg,
                              yleg,
                              main,
                              ...
                              )
{
  MM.k <- attr(x, "sample.size")
  kmax <- length(MM.k) - 1
  kmax.insample <- max((0:kmax)[MM.k > 0])
  
  if (missing(k.cond)) k.cond <- 0:kmax.insample
  not.sampled <- k.cond > kmax.insample
  k.cond <- k.cond[!not.sampled]

  num.plots <- length(k.cond)
  dims <- switch(num.plots, c(1, 1),
                            c(1, 2),
                            c(2, 2), c(2, 2),
                            c(2, 3), c(2, 3),
                            c(3, 3), c(3, 3), c(3, 3),
                            c(3, 4), c(3, 4), c(3, 4),
                            c(4, 4), c(4, 4), c(4, 4), c(4, 4),
                            c(5, 4), c(5, 4), c(5, 4), c(5, 4))
  if (dim.plot)        par(mfrow = dims)
  if (dim.plot & over) par(mfrow = c(1, 1))
  if (missing(lty)) lty <- 1:num.plots
  

  empty.plot <- function(mess = "Not in the sample")
  {
    plot(0:100, 0:100, type = "n", xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "")
    text(50, 90, labels = mess, adj = 1)
  }    
  

  M <- MM.k[1]
  M.k <- MM.k[-1]
  nulls <- M.k == 0
  
  if (over & k.cond[1] & !sum(M.k[k.cond])){
    empty.plot(mess = "Conditional densities you asked are not in the sample.")
    return(invisible(x))
  }    

  if (standard) mcmc <- x$standard
  else          mcmc <- x$unstandard

  if (is.character(mcmc)) stop("You have to first compute McMC averages.")
  grid <- mcmc$grid
  
  remove <- c(TRUE, FALSE, nulls)
  mcmc0 <- mcmc[!remove]
  if (missing(xlim)) xlim <- range(grid)
  if (missing(ylim)) ylim <- c(0, max(sapply(mcmc0, max, na.rm = TRUE), na.rm = TRUE))

  legg <- character(0)
  relatM <- round((M.k/M)*100, 2)
  if (over){
    i <- 1
    if (!k.cond[1]){
      dens <- mcmc$unconditional
      legg <- c(legg, paste("Uncond.,  ", "M = ", M, sep = ""))
    }
    else{
      while (!M.k[k.cond[i]]){
        legg <- c(legg, paste("k = ", k.cond[i], "   (", relatM[k.cond[i]], " %)",  sep = ""))
        i <- i + 1 
      }  
      dens <- mcmc[[paste("k = ", k.cond[i], sep = "")]]
      legg <- c(legg, paste("k = ", k.cond[i], "   (", relatM[k.cond[i]], " %)",  sep = ""))      
    }
    plot(grid, dens, type = type, bty = bty, xlab = xlab, ylab = ylab, lty = lty[i], xlim = xlim, ylim = ylim)
    if (i < length(k.cond)){
      for (j in (i+1):length(k.cond)){
        if (!M.k[k.cond[j]]){
          legg <- c(legg, paste("k = ", k.cond[j], "   (", relatM[k.cond[j]], " %)", sep = ""))
          next
        }        
        dens <- mcmc[[paste("k = ", k.cond[j], sep = "")]]
        legg <- c(legg, paste("k = ", k.cond[j], "   (", relatM[k.cond[j]], " %)", sep = ""))      
        lines(grid, dens, lty = lty[j])  
      }
    }      
    if (missing(xleg)) xleg <- min(grid)
    if (missing(yleg)) yleg <- ylim[2] - 0.1*(ylim[2] - ylim[1]) 
    if (alegend) legend(xleg, yleg, legend = legg, lty = lty, xjust = 0, yjust = 1, bty = "n")
    if (missing(main)) title(main = "McMC averages of the density")
    else               title(main = main)
  }    
  else{    ## not over
    for (k in k.cond){
      if (k == 0){
        dens <- mcmc$unconditional
        titul <- "unconditional"
        subb <- paste("M = ", M, sep = "")
        if (missing(xlim)) xlim <- range(grid)
        if (missing(ylim)) ylim <- range(dens)
        plot(grid, dens, type = type, bty = bty, xlab = xlab, ylab = ylab, lty = lty[1], xlim = xlim, ylim = ylim)
      }      
      else{
        dens <- mcmc[[paste("k = ", k, sep = "")]]
        titul <- paste("k = ", k, sep = "")
        subb <- paste("M = ", M.k[k], "   (", relatM, " %)", sep = "")
        if (!M.k[k]){
          empty.plot()
          subb <- ""
        }        
        else{
          if (missing(xlim)) xlim <- range(grid)
          if (missing(ylim)) ylim <- range(dens)          
          plot(grid, dens, type = type, bty = bty, xlab = xlab, ylab = ylab, lty = lty[1], xlim = xlim, ylim = ylim)
        }  
      }  
      title(main = titul, sub = subb)
    }
  }    

  return(invisible(x))  
}  




## Subfunction common for all 'bayessurvreg' functions
##  - extract a design information

## 18/03/2004

bayessurvreg.design <- function(m, formula, random, data, transform, dtransform)
{

   tempF <- c("", "formula", "data", "subset", "na.action")
   mF <- m[match(tempF, names(m), nomatch=0)]
   mF[[1]] <- as.name("model.frame")
   special <- c("cluster")
   TermsF <- if(missing(data)) terms(formula, special)
             else              terms(formula, special, data = data)
   mF$formula <- TermsF
   mF <- eval(mF, parent.frame())
 
   ## Response matrix
   Y <- model.extract(mF, "response")
   Yinit <- Y
   if (!inherits(Y, "Surv"))
      stop("Response must be a survival object. ")
   type <- attr(Y, "type")
   if (type == 'counting') stop ("Invalid survival type ('counting' is not implemented). ")
   n <- nrow(Y)
   nY <- ncol(Y)

   ## Cluster indicators   
   cluster <- attr(TermsF, "specials")$cluster
   dropx <- NULL
   if (length(cluster)) {
     tempc <- untangle.specials(TermsF, "cluster", 1:10)
     ord <- attr(TermsF, "order")[tempc$terms]
     if (any(ord > 1))
         stop("Cluster can not be used in an interaction")
     cluster <- strata(mF[, tempc$vars], shortlabel = TRUE)
     dropx <- tempc$terms
   }
   if (length(dropx)){
      newTermsF <- TermsF[-dropx]
      attr(newTermsF, "intercept") <- attr(TermsF, "intercept")  ## I do not why but the command on the previous row
                                                                 ## sets attr(newTermsF, "intercept") always to 1,
                                                                 ## irrespective of what attr(TermsF, "intercept") was...
   }
   else
      newTermsF <- TermsF

   ## Design matrix for both fixed and random effects X
   ## (finally, always without the intercept)
     ## Temporarily, include always at least the intercept
     ##  (to get nice rownames and initial estimates)
   attr(newTermsF, "intercept") <- 1
   Xinit <- model.matrix(newTermsF, mF)
   rnamesX <- row.names(Xinit)
   cnamesX <- colnames(Xinit)
   
     ## Finally, intercept will be always removed
   nX <- ncol(Xinit) - 1
   cnamesX <- cnamesX[-1]
   if (nX){
     X <- Xinit
     X <- X[,-1]                                      ## removal of the intercept
     attr(X, "assign") <- attr(Xinit, "assign")[-1]   ## removal of the intercept
     attr(X, "contrasts") <- attr(Xinit, "contrasts")
   }
   else{
     X <- NULL
   }     
   indb <- if (nX) rep(-1, nX)                          ## initially, all effects are only fixed
           else    0

   
   ## Design matrix for random effects 
   randomInt <- FALSE
   nrandom <- 0
   if (!missing(random)){
     if (!length(cluster)) stop ("You have to indicate clusters. ")
     tempR <- c("", "random", "data", "subset", "na.action")
     mR <- m[match(tempR, names(m), nomatch=0)]
     mR[[1]] <- as.name("model.frame")
     names(mR)[2] <- "formula"
     TermsR <- if(missing(data)) terms(random)
               else              terms(random, data = data)
     lTR <- length(attr(TermsR, "variables"))
     if (lTR == 1 & !attr(TermsR, "intercept")){        ## do nothing, in reality no random terms in the model
       names.random <- character(0)
     }
     else{
       if (lTR == 1 & attr(TermsR, "intercept")){        ## the only random term is the intercept
         randomInt <- TRUE
         names.random <- character(0)
       }
       else{
         mR$formula <- TermsR
         mR <- eval(mR, parent.frame())
         if (attr(TermsR, "intercept")){
           randomInt <- TRUE
           attr(TermsR, "intercept") <- 0     ## remove it from the design matrix of random effects
         }         
         names.random <- colnames(model.matrix(TermsR, mR))
       }
       nrandom <- 1*randomInt + length(names.random)
       if (sum(names.random %in% cnamesX) != nrandom - 1*randomInt) stop("Each random effect has to have also its fixed counterpart.")
       find.indeces <- function(all.eff){     ## Find indeces of random effects in a design matrix to be passed to C++
         where <- names.random %in% all.eff
         if (!sum(where)) return (-1)
         if (sum(where) > 1) stop("Error, contact the author.")
         index <- (1:length(names.random))[where]
         if (!randomInt) index <- index - 1
         return(index)
       }
       if (nX) indb <- as.numeric(apply(matrix(cnamesX, ncol = 1), 1, find.indeces))
     }       
   }
   if (randomInt) names.random <- c("(Intercept)", names.random)
   nfixed <- nX - (nrandom - 1*randomInt)

   ## Give indeces of factors in the design matrix, it was used to define MH blocks in the earlier version of this program
   n.factors <- 0
   n.in.factors <- NULL
   factors <- NULL
   if (nX){
     temp <- attr(X, "assign")
     if (length(temp) == 1) factors <- 0
     else{
       factors <- numeric(length(temp))
       n.in.factors <- numeric(0)
       temp  <- temp - c(0, temp[1:(length(temp)-1)])
       i <- length(temp)
       while (i >= 1){
         if (temp[i] == 0){
           n.factors <- n.factors + 1
           factors[i] <- n.factors
           n.in.factor <- 1
           while (temp[i-1] == 0){
             i <- i - 1
             factors[i] <- n.factors
             n.in.factor <- n.in.factor + 1
           }
           i <- i - 1
           factors[i] <- n.factors
           n.in.factors <- c(n.in.factor + 1, n.in.factors)
         }
         else
           n.in.factors <- c(1, n.in.factors)         
         i <- i - 1
       }         
     }
     if (length(temp) != nX) stop("Something is wrong, contact the author.")
   }

   ## Sort everything according to the cluster indicator
   ##   and find the numbers of observations per cluster
   if (length(cluster)){
     ordering <- order(cluster)
     Y <- Y[ordering, ]
     cluster <- cluster[ordering]
     rnamesX <- rnamesX[ordering]
     if (nX){
       namesX <- cnamesX
       if (nX == 1) X <- matrix(X[ordering], ncol = 1)
       else         X <- as.matrix(X[ordering, ])
       colnames(X) <- namesX
     }
     ncluster <- length(attr(cluster, "levels"))
     helpf <- function(cl){return(sum(cluster %in% attr(cluster, "levels")[cl]))}
     nwithin <- apply(matrix(1:ncluster, ncol = 1), 1, "helpf")
   }
   else{
     if (nX) X <-  as.matrix(X[,])
     cluster <- 1:n
     ncluster <- n
     nwithin <- rep(1, n)
   }     
   
   ## Transform the response
   if (type == 'interval') {
      if (any(Y[,3]==3)) Y <- cbind(eval(call("transform", Y[,1:2])), Y[,3])
      else               Y <- cbind(eval(call("transform", Y[,1])), Y[,3])
   }
   else if (type == 'left'){
          Y <- cbind(eval(call("transform", Y[,1])), 2-Y[,2])   ## change 0 indicator into 2 indicating left censoring
        }
        else  ## type = 'right' or 'interval2'
           Y <- cbind(eval(call("transform", Y[,1])), Y[,2])

   if (!all(is.finite(Y))) stop("Invalid survival times for this distribution (infinity on log-scale not allowed). ")

   design <- list(n = n, ncluster = ncluster, nwithin = nwithin, nY = nY, nX = nX,
                  nfixed = nfixed, nrandom = nrandom, randomInt = randomInt,
                  Y = Y, X = X, Yinit = Yinit, Xinit = Xinit, cluster = cluster, indb = indb,
                  rnamesX = rnamesX, names.random = names.random, factors = factors, n.factors = n.factors, n.in.factors = n.in.factors)

   return(design)
}  
## Subfunction for 'bayessurvreg' functions
##  -> just to make it more readable
##
## Manipulation with the prior and proposal specification
##  for the regression parameters and means of random effects
##
##
## factors ..... a vector of length nX with 0 on places of columns which are not factors
##               and with 1, 2, ... on places of factors
##               * all columns of X corresponding to one factor have same index
##               * the first factor has teh highest number
## n.factors ... number of factor variables

bayessurvreg.priorBeta <- function(prior.beta, nX, indb, factors, n.factors, n.in.factors)
{
  if (!nX){
    priori <- c(0, 0, 0, 0, 0, 0, 0)
    priord <- c(0, 0, 0, 0, 0, 0, 0, 0)
    
    pdi.beta <- list(double = priord, integer = priori)
    attr(pdi.beta, "prior.beta") <- list()    
    return(pdi.beta)
  }    
  
  if (is.null(prior.beta$mean.prior)) stop("Prior means for betas must be specified.")
  if (is.null(prior.beta$var.prior)) stop("Prior vars for betas must be specified.")
  if (length(prior.beta$mean.prior) != nX) stop("Incorrect length of a vector of prior means for betas.")
  if (length(prior.beta$var.prior) != nX) stop("Incorrect length of a vector of prior vars for betas.")  
  if (sum(is.na(prior.beta$mean.prior))) stop("Prior means for betas must not be missing.")
  if (sum(is.na(prior.beta$var.prior))) stop("Prior vars for betas must not be missing.")  
  if (sum(prior.beta$var.prior <= 0)) stop("Prior vars for betas must be all positive.")    

  if (is.null(prior.beta$blocks) | is.null(prior.beta$blocks$ind.block)){

    ## One block with fixed effects, one block with means of random effects
    prior.beta$blocks <- list()
    fixed <- (1:nX)[indb == -1]
    random <- (1:nX)[indb > 0]

    prior.beta$blocks$ind.block <- list()
    prior.beta$blocks$cov.prop <- list()
    nBlocks <- 0
    if (length(fixed) > 0){
      nBlocks <- nBlocks + 1
      prior.beta$blocks$ind.block[[nBlocks]] <- fixed
      prior.beta$blocks$cov.prop[[nBlocks]] <- numeric((length(fixed) * (1 + length(fixed)))/2)
    }
    if (length(random) > 0){
      nBlocks <- nBlocks + 1
      prior.beta$blocks$ind.block[[nBlocks]] <- random
      prior.beta$blocks$cov.prop[[nBlocks]] <- numeric((length(random) * (1 + length(random)))/2)
    }      
    prior.beta$type.upd <- rep("gibbs", nBlocks)
    
    ## ===== OLDER VERSION =====
    ## Each parameter in one block, except factors
    ## cov.prop = diag(prior variance) for each block
##    prior.beta$blocks <- list()
##    nBlocks <- length(n.in.factors)
##    nInBlock <- n.in.factors
##    cumn <- c(0, cumsum(nInBlock))

##    prior.beta$blocks$ind.block <- list()
##    prior.beta$blocks$cov.prop <- list()
##    for (i in 1:nBlocks){
##      prior.beta$blocks$ind.block[[i]] <- (cumn[i]+1):cumn[i+1]
##      covmat <- diag(prior.beta$var.prior[(cumn[i]+1):cumn[i+1]], nrow = nInBlock[i])      
##      prior.beta$blocks$cov.prop[[i]] <- covmat[lower.tri(covmat, diag = TRUE)]
##    }    

  }

  if (is.null(prior.beta$blocks$cov.prop)){
    prior.beta$blocks$cov.prop <- list()
  }     
  
  nBlocks <- length(prior.beta$blocks$ind.block)
  if (is.null(prior.beta$type.upd)) prior.beta$type.upd <- rep("gibbs", nBlocks)
  if (length(prior.beta$type.upd) != nBlocks) stop("Incorrect prior.beta$type.upd parameter.")

  typeUpd <- pmatch(prior.beta$type.upd, table = c("random.walk.metropolis", "adaptive.metropolis", "gibbs"), nomatch = 0, duplicates.ok = TRUE)
  typeUpd[typeUpd == 0] <- 3        ## no matching ==> update using gibbs
  typeUpd <- typeUpd - 1            ## now: 0 = random walk, 1 = adaptive, 2 = gibbs
  
  for (i in 1:nBlocks){
    if (typeUpd[i] == 2){
      lb <- length(prior.beta$blocks$ind.block[[i]])
      prior.beta$blocks$cov.prop[[i]] <- numeric((lb*(1+lb))/2)
    }      
  }    
  
  if (nBlocks != length(prior.beta$blocks$cov.prop)) stop("Not all beta blocks have defined a covariance matrix for the proposal")

  nInBlock <- sapply(prior.beta$blocks$ind.block, length)
  if (sum(nInBlock) != nX) stop("Some beta parameters are not assigned to any block or to more blocks.")

  indBlockLV <- unlist(prior.beta$blocks$ind.block)
  if (length(indBlockLV) != nX) stop("Some beta parameters are assigned to more blocks")
  if (sum(indBlockLV %in% 1:nX) != nX) stop("Incorrect prior.beta$ind.block parameter.")

  covparLV <- unlist(prior.beta$blocks$cov.prop)
  lcovparLV <- sum(0.5*nInBlock*(nInBlock+1))
  if (sum(is.na(covparLV))) stop("Incorrect prior.beta$cov.prop parameter.")
  if (length(covparLV) != lcovparLV) stop("Incorrect prior.beta$cov.prop parameter.")

  if (is.null(prior.beta$mean.sampled)) prior.beta$mean.sampled <- rep(0, nX)
  if (length(prior.beta$mean.sampled) != nX) stop("Incorrect prior.beta$mean.sampled parameter.")
  if (sum(is.na(prior.beta$mean.sampled))) stop("Incorrect prior.beta$mean.sampled parameter.")
  meanSampled <- prior.beta$mean.sampled
  
  if (is.null(prior.beta$eps.AM)) prior.beta$eps.AM <- rep(0.05, nBlocks)
  eps <- prior.beta$eps.AM
  if (length(eps) != nBlocks) stop("Incorrect prior.beta$eps.AM parameter.")
  if (sum(is.na(eps))) stop("Incorrect prior.beta$eps.AM parameter.")
  if (sum(eps < 0)) stop("Incorrect prior.beta$eps.AM parameter.")  
  
  if (is.null(prior.beta$sd.AM)) prior.beta$sd.AM <- (2.4*2.4)/(1:max(nInBlock))
  sdNum <- prior.beta$sd.AM
  if (length(sdNum) < max(nInBlock)) stop("Too short prior.beta$sd.AM parameter.")
  sdNum <- sdNum[1:max(nInBlock)]
  if (sum(is.na(sdNum))) stop("Incorrect prior.beta$sd.AM parameter.")
  if (sum(sdNum < 0)) stop("Incorrect prior.beta$sd.AM parameter.")    

  if (is.null(prior.beta$weight.unif)) prior.beta$weight.unif <- rep(0.5, nBlocks)
  if (sum(is.na(prior.beta$weight.unif))) stop("Incorrect prior.beta$weight.unif parameter.")
  if (length(prior.beta$weight.unif) != nBlocks) stop("Incorrect prior.beta$weight.unif parameter.")  
  prior.beta$weight.unif[prior.beta$weight.unif < 0] <- 0
  prior.beta$weight.unif[prior.beta$weight.unif > 1] <- 1  
  weightUnif <- prior.beta$weight.unif

  if (is.null(prior.beta$half.range.unif)) prior.beta$half.range.unif <- 0.5*sqrt(12*prior.beta$var.prior)
  if (sum(is.na(prior.beta$half.range.unif))) stop("Incorrect prior.beta$half.range.unif parameter.")
  if (length(prior.beta$half.range.unif) != nX) stop("Incorrect prior.beta$half.range.unif parameter.")  
  halfRangeUnif <- prior.beta$half.range.unif
  
  ## Put everything to long vectors
  ##  for indBlockLV, do R --> C++ transformation
  priord <- c(prior.beta$mean.prior, prior.beta$var.prior, meanSampled, halfRangeUnif, covparLV,  weightUnif, eps, sdNum)
  priori <- c(nBlocks, nX, nInBlock, max(nInBlock), lcovparLV, indBlockLV - 1, typeUpd)

  pdi.beta <- list(double = priord, integer = priori)
  attr(pdi.beta, "prior.beta") <- prior.beta

  return(pdi.beta)  
}  
## Subfunction for 'bayessurvreg' functions
##  -> just to make it more readable
##
## Manipulation with proposal and prior specification
##  for random effects
##
bayessurvreg.priorb <- function(prior.b, nrandom, ncluster, toler.chol)
{
  thispackage = "bayesSurv"

  if (!nrandom){
    priori <- c(0, 0, 0, 0, 0, 0, 0, 0)
    priord <- c(0, 0, 0, 0, 0)
    
    pdi.b <- list(double = priord, integer = priori)
    attr(pdi.b, "prior.b") <- list()    
    return(pdi.b)
  }

  ## Type of prior for D
  ## =====================
  if (is.null(prior.b$prior.D)) prior.b$prior.D <- "inv.wishart"
  priorD <-  pmatch(prior.b$prior.D, table = c("inv.wishart", "sduniform"), nomatch = 0)
  if (priorD == 0){
    prior.b$prior.D <- "inv.wishart"
    warning("Non-matching prior.b$prior.D changed to inv.wishart")
    priorD <- 1    
  }
  if (priorD != 1 & nrandom > 1){
    prior.b$prior.D <- "inv.wishart"
    warning("Dimension of the random effect > 1: prior.b$prior.D changed to inv.wishart")
  }    
  priorD <- priorD - 1     ## now: 0 = inv.wishart, 1 = sduniform
  
  ## Type of proposal etc.
  ## =====================
  if (is.null(prior.b$type.upd)) prior.b$type.upd <- "gibbs"
  if (length(prior.b$type.upd) != 1) stop("Incorrect prior.b$type.upd parameter.")
  typeUpd <- pmatch(prior.b$type.upd, table = c("random.walk.metropolis", "adaptive.metropolis", "gibbs"), nomatch = 0)
  if (typeUpd == 2 | typeUpd == 0){
    warning("Update type for random effects changed to 'gibbs'.")
    typeUpd <- 3
  }    
  typeUpd <- typeUpd - 1            ## now: 0 = random walk, 1 = adaptive, 2 = gibbs

  if (typeUpd == 2){
    prior.b$blocks <- NULL
    prior.b$weight.unif <- NULL
    prior.b$half.range.unif <- NULL

    nBlocks <- 1
    nInBlock <- nrandom
    lcovparLV <- 0.5*nInBlock*(nInBlock+1)
    indBlockLV <- 1:nrandom
    covparLV <- numeric(lcovparLV)
    halfRangeUnif <- numeric(nrandom)
    weightUnif <- numeric(nBlocks)    
  }
  else{
    if (is.null(prior.b$blocks)) stop("Blocks for random effects must be specified.")
    if (is.null(prior.b$blocks$ind.block)) stop("Blocks for random effects must be specified.")
    if (is.null(prior.b$blocks$cov.prop)) stop("Covariance matrices for proposal for random effects must be specified.")

    nBlocks <- length(prior.b$blocks$ind.block)
    if (nBlocks != length(prior.b$blocks$cov.prop)) stop("Not all random effects blocks have defined a covariance matrix for the proposal")

    nInBlock <- sapply(prior.b$blocks$ind.block, length)
    if (sum(nInBlock) != nrandom) stop("Some random effects are not assigned to any block or to more blocks.")

    indBlockLV <- unlist(prior.b$blocks$ind.block)
    if (length(indBlockLV) != nrandom) stop("Some b parameters are assigned to more blocks")
    if (sum(indBlockLV %in% 1:nrandom) != nrandom) stop("Incorrect prior.b$ind.block parameter.")

    covparLV <- unlist(prior.b$blocks$cov.prop)
    lcovparLV <- sum(0.5*nInBlock*(nInBlock+1))
    if (sum(is.na(covparLV))) stop("Incorrect prior.b$cov.prop parameter.")
    if (length(covparLV) != lcovparLV) stop("Incorrect prior.b$cov.prop parameter.")
    
    if (is.null(prior.b$weight.unif)) prior.b$weight.unif <- rep(0.5, nBlocks)
    if (sum(is.na(prior.b$weight.unif))) stop("Incorrect prior.b$weight.unif parameter.")
    if (length(prior.b$weight.unif) != nBlocks) stop("Incorrect prior.b$weight.unif parameter.")  
    prior.b$weight.unif[prior.b$weight.unif < 0] <- 0
    prior.b$weight.unif[prior.b$weight.unif > 1] <- 1  
    weightUnif <- prior.b$weight.unif

    if (is.null(prior.b$half.range.unif)) stop("half.range.unif for random effects must be specified.")
    if (sum(is.na(prior.b$half.range.unif))) stop("Incorrect prior.b$half.range.unif parameter.")
    if (length(prior.b$half.range.unif) != nrandom) stop("Incorrect prior.b$half.range.unif parameter.")  
    halfRangeUnif <- prior.b$half.range.unif
  }

  ## Parameters of prior for covariance matrix D
  ## ===========================================
  if (is.null(prior.b$df.D)) prior.b$df.D <- nrandom + 2
  if (is.null(prior.b$scale.D)){
    if (priorD == 1) stop("Upper limit of the uniform prior for std. dev. of the random effect must be given.")
    prior.b$scale.D <- 0.002*diag(nrandom)
    prior.b$scale.D <- prior.b$scale.D[lower.tri(prior.b$scale.D, diag = TRUE)]
  }
  if (!is.null(prior.b$df.D) && priorD == 0)
    if (prior.b$df.D <= nrandom - 1) stop("Too low prior degrees of freedom for D matrix.")

  if (priorD == 0){
    indD <- 0:(nrandom - 1)
    diagI <- (indD * (2*nrandom - indD + 1)) / 2
    lD <- (nrandom * (nrandom + 1))/2
    if (length(prior.b$scale.D) != lD) stop("Incorrect prior scale matrix for D matrix.")
    cholD <- .C("cholesky", as.double(prior.b$scale.D), rank = integer(1), as.integer(nrandom), as.integer(diagI), as.double(toler.chol),
                PACKAGE = thispackage)
    if (cholD$rank < nrandom) stop("Prior scale matrix for D matrix is not positive definite.")
  }
  if (priorD == 1){
    prior.b$scale.D <- prior.b$scale.D[1]
    if (prior.b$scale.D[1] <= 0) stop("Upper limit of the uniform prior for std. dev. of the random effect must be positive.")
    prior.b$df.D <- 1/(prior.b$scale.D*prior.b$scale.D)
  }    

  priordD <- c(prior.b$df.D, prior.b$scale.D)
  if (priorD == 0) names(priordD) <- c("df.D", paste("scale.D", 1:lD))
  if (priorD == 1) names(priordD) <- c("1/B.sq", "B")  

  ## Put everything to long vectors
  ##  for indBlockLV, do R --> C++ transformation  
  ## =============================================
  priori <- c(nrandom, ncluster, priorD, typeUpd, nBlocks, nInBlock, lcovparLV, indBlockLV - 1)
  priord <- c(priordD, covparLV, halfRangeUnif, weightUnif)

  pdi.b <- list(double = priord, integer = priori)
  attr(pdi.b, "prior.b") <- prior.b

  return(pdi.b)
}  
bayessurvreg1 <- function(
     formula,
     random,
     data = parent.frame(),
     subset,
     na.action = na.fail,
     x = FALSE,
     y = FALSE,                          
     onlyX = FALSE,
     nsimul = list(niter = 10, nthin = 1, nburn = 0, nnoadapt = 0, nwrite = 10),
     prior = list(kmax = 5, k.prior = "poisson", poisson.k = 3,
                  dirichlet.w = 1,
                  mean.mu = NULL, var.mu = NULL,
                  shape.invsig2 = 1.5, shape.hyper.invsig2 = 0.8, rate.hyper.invsig2 = NULL,
                  pi.split = NULL, pi.birth = NULL,
                  Eb0.depend.mix = FALSE),
     prior.beta,
     prior.b,
     prop.revjump,
     init = list(iter = 0, mixture = NULL, beta = NULL, b = NULL, D = NULL,
                 y = NULL, r = NULL, otherp = NULL, u = NULL),
     store = list(y = TRUE, r = TRUE, b = TRUE, u = TRUE, MHb = FALSE, regresres = FALSE),
     dir = getwd(),
     toler.chol = 1e-10,
     toler.qr = 1e-10,
     ...)
{
   thispackage = "bayesSurv"
  
   transform = function(t){log(t)}
   dtransform = function(t){1/t}
  
   sim.to.R <- FALSE                                  ## this is here for a compatibility with an older code
   store <- bayessurvreg1.checkStore(store)
   
   ## Give a function call to be recorded in a resulting object.
   call <- match.call(expand.dots = TRUE)

   ## Extract all the design information from the function call
   m <- match.call(expand.dots = FALSE)
   des <- bayessurvreg.design(m, formula, random, data, transform, dtransform)
   if (onlyX) return (des$X)   

   
   ## =========================================================
   ## Manipulate with initial values and the prior information
   ## =========================================================
   priordi <- bayessurvreg1.priorInit(prior, init, des$Yinit, des$Xinit, des$n, des$nX, des$nrandom, des$ncluster, des$indb, des$randomInt, toler.chol)
   prior <- attr(priordi, "prior")
   init <- attr(priordi, "init")
   
   if (missing(prop.revjump)) prop.revjump <- list()
   revjumpdi <- bayessurvreg1.revjump(prop.revjump)
   prop.revjump <- attr(revjumpdi, "prop.revjump")
   
   if (missing(prior.beta)) prior.beta <- list()
   betadi <- bayessurvreg.priorBeta(prior.beta, des$nX, des$indb, des$factors, des$n.factors, des$n.in.factors)
   prior.beta <- attr(betadi, "prior.beta")   
   prior.beta.noadapt <- attr(betadi, "prior.beta")
   if (des$nX){     
     adapts <- pmatch(prior.beta$type.upd, table = c("adaptive.metropolis"), nomatch = 0, duplicates.ok = TRUE)       
     n.adapt <- sum(adapts)
     prior.beta.noadapt$type.upd[adapts == 1] <- "random.walk.metropolis"
   }
   else{
     n.adapt <- 0
   }       
   betadi.noadapt <- bayessurvreg.priorBeta(prior.beta.noadapt, des$nX, des$indb, des$factors, des$n.factors, des$n.in.factors) 
   
   if (missing(prior.b)) prior.b <- list()
   bdi <- bayessurvreg.priorb(prior.b, des$nrandom, des$ncluster, toler.chol)
   prior.b <- attr(bdi, "prior.b")
   

   ## ===================================================================
   ## Compute quantities to determine the space needed to be allocated
   ##   and numbers of iterations in different phases
   ## ===================================================================
   if (nsimul$nburn >= nsimul$niter) nsimul$nburn <- nsimul$niter - 1
   if (nsimul$nburn < 0) nsimul$nburn <- 0

   if (n.adapt == 0) nsimul$nnoadapt <- nsimul$nburn
   if (nsimul$nnoadapt > nsimul$nburn) nsimul$nnoadapt <- nsimul$nburn
   if (nsimul$nnoadapt < 0) nsimul$nnoadapt <- 0
   
   if (nsimul$nburn == 0) nruns <- 1
   else                   if (nsimul$nnoadapt == nsimul$nburn | nsimul$nnoadapt == 0) nruns <- 2
                          else                                                        nruns <- 3

   nrun <- numeric(3)
   nrun[3] <- nsimul$niter - nsimul$nburn
   nrun[2] <- nsimul$nburn - nsimul$nnoadapt
   nrun[1] <- nsimul$nnoadapt

   nwrite.run <- nrun
   nwrite.run[nsimul$nwrite <= nrun] <- nsimul$nwrite   
   max.nwrite <- max(nwrite.run)

   if (!des$nrandom){ store$b <- FALSE;  store$MHb <- FALSE}
   row.need <- ifelse(sim.to.R, max(nsimul$nburn, nafterburn), max.nwrite)
   
   ## =====================================================================================
   ## Write headers to files with stored values
   ## =====================================================================================
   bayessurvreg1.writeHeaders(dir, prior, store, des$nX, des$X, des$names.random, des$ncluster, des$nrandom, des$rnamesX,
                              unique(des$cluster), betadi$integer[1], bdi$integer[4])

   ## =========================================================
   ## Combine similar parameters into one vector
   ## =========================================================
   dims <- c(des$n, des$ncluster, des$nwithin, des$nY, des$nX, des$nfixed, des$nrandom, 1*des$randomInt, row.need)
   storeV <- c(store$y, store$r, store$b, store$u, store$MHb, store$regresres)
   nsimul.run1 <- c(nrun[1], nsimul$nthin, nwrite.run[1])
   nsimul.run2 <- c(nrun[2], nsimul$nthin, nwrite.run[2])
   nsimul.run3 <- c(nrun[3], nsimul$nthin, nwrite.run[3])   
   tolers <- c(toler.chol, toler.qr)

   ## =====================================
   ## Keep some parameters to be returned
   ## =====================================
   keep.init <- init
   
   cat("Simulation started on                       ", date(), "\n", sep = "")      
     ## Run without adaptation of a proposal covariance matrices
     ## Either the whole burn up or first part of burn up
   if (nruns == 3 | (nruns == 2 & nsimul$nnoadapt == nsimul$nburn)){
     fit <- .C("bayessurvreg1", as.character(dir),
                                dims = as.integer(dims),               
                                Y = as.double(des$Y),
                                X = as.double(des$X),
                                indb = as.integer(des$indb),               
                                iter = as.integer(init$iter),
                                loglik = as.double(c(0, 0)),
                                mixture = as.double(init$mixture),
                                mixmoment = as.double(c(0, 0)),
                                beta = as.double(init$beta),
                                b = as.double(init$b),
                                D = as.double(init$D),
                                r = as.integer(init$r),
                                Ys = as.double(init$y),       
                                otherp = as.double(init$otherp),
                                u = as.double(init$u),
                                prior.pari = as.integer(priordi$integer),
                                prior.pard = as.double(priordi$double),
                                revJump.pari = as.integer(revjumpdi$integer),
                                revJump.pard = as.double(revjumpdi$double),
                                prior.betai = as.integer(betadi.noadapt$integer),
                                prior.betad = as.double(betadi.noadapt$double),
                                prior.bi = as.integer(bdi$integer),
                                prior.bd = as.double(bdi$double),
                                nsimul = as.integer(nsimul.run1),
                                store = as.integer(storeV),
                                tolers = as.double(tolers),
                                err = integer(1),
               PACKAGE = thispackage)     
     if (fit$err != 0) stop ("Something went wrong during the simulation.")
     cat("Simulation without adaptation finished on   ", date(), "   (iteration ", fit$iter, ")", "\n", sep = "")
   
       ## Recalculate intial proposal covariance matrices
       ##  * do not change it if sample covariance matrix is not positive definite
     if (des$nX){
       betasam <- read.table(paste(dir, "/beta.sim", sep = ""), header = TRUE)
       betaaver <- apply(betasam, 2, mean)
       prior.beta$mean.sampled <- betaaver
       to.adapt <- (1:betadi$integer[1])[adapts == 1]
       if (length(to.adapt) > 0){
         for (i in 1:length(to.adapt)){
           colsX <- prior.beta$blocks$ind.block[[to.adapt[i]]]
           sample <- as.data.frame(betasam[, colsX])
           colnames(sample) <- colnames(betasam)[colsX]
           covmat <- sampleCovMat(sample)
           covmat <- prior.beta$sd.AM[length(colsX)] * (covmat + prior.beta$eps.AM[to.adapt[i]]*diag(length(colsX)))
           vind <- 0:(dim(covmat)[1] - 1)
           diagI <- (vind * (2*dim(covmat)[1] - vind + 1)) / 2
           covmatt <- covmat[lower.tri(covmat, diag = TRUE)]
           chol <- .C("cholesky", A = as.double(covmatt), rank = integer(1), as.integer(dim(covmat)[1]),
                                    as.integer(diagI), as.double(toler.chol),
                      PACKAGE = thispackage)
           if (chol$rank == dim(covmat)[1]){
             prior.beta$blocks$cov.prop[[to.adapt[i]]] <- as.numeric(covmatt)
           } 
           else{
             warning("Sample covariance matrix after no adapt period was not positive definite.")
           }
         }
       }       
       betadi <- bayessurvreg.priorBeta(prior.beta, des$nX, des$indb, des$factors, des$n.factors, des$n.in.factors)
     }       
     
       ## Give new initials
     init$iter <- fit$iter;   init$mixture <- fit$mixture;   init$beta <- fit$beta;
     init$b <- fit$b;         init$D <- fit$D;               init$r <- fit$r;
     init$y <- fit$Ys;        init$otherp <- fit$otherp;     init$u <- fit$u;

       ## Rewrite sampled values by new files
     bayessurvreg1.writeHeaders(dir, prior, store, des$nX, des$X, des$names.random, des$ncluster, des$nrandom, des$rnamesX,
                                unique(des$cluster), betadi$integer[1], bdi$integer[4])
   }     

     ## Burn up with adaptation
     ## Either the whole burn up or the second part of burn up
   if (nruns == 3 | (nruns == 2 & nsimul$nnoadapt == 0)){
     fit <- .C("bayessurvreg1", as.character(dir),
                                dims = as.integer(dims),               
                                Y = as.double(des$Y),
                                X = as.double(des$X),
                                indb = as.integer(des$indb),               
                                iter = as.integer(init$iter),
                                loglik = as.double(c(0, 0)),
                                mixture = as.double(init$mixture),
                                mixmoment = as.double(c(0, 0)),               
                                beta = as.double(init$beta),
                                b = as.double(init$b),
                                D = as.double(init$D),
                                r = as.integer(init$r),
                                Ys = as.double(init$y),       
                                otherp = as.double(init$otherp),
                                u = as.double(init$u),
                                prior.pari = as.integer(priordi$integer),
                                prior.pard = as.double(priordi$double),
                                revJump.pari = as.integer(revjumpdi$integer),
                                revJump.pard = as.double(revjumpdi$double),
                                prior.betai = as.integer(betadi$integer),
                                prior.betad = as.double(betadi$double),
                                prior.bi = as.integer(bdi$integer),
                                prior.bd = as.double(bdi$double),
                                nsimul = as.integer(nsimul.run2),
                                store = as.integer(storeV),
                                tolers = as.double(tolers),
                                err = integer(1),
               PACKAGE = thispackage)
     if (fit$err != 0) stop ("Something went wrong during the simulation.")
     cat("Burn-up finished on                         ", date(), "   (iteration ", fit$iter, ")", "\n", sep = "")
     
       ## Give new initials
     init$iter <- fit$iter;   init$mixture <- fit$mixture;   init$beta <- fit$beta;
     init$b <- fit$b;         init$D <- fit$D;               init$r <- fit$r;
     init$y <- fit$Ys;        init$otherp <- fit$otherp;     init$u <- fit$u;     

       ## Rewrite sampled values by new files
     bayessurvreg1.writeHeaders(dir, prior, store, des$nX, des$X, des$names.random, des$ncluster, des$nrandom, des$rnamesX,
                                unique(des$cluster), betadi$integer[1], bdi$integer[4])
   }     
   
     ## Main simulation
   fit <- .C("bayessurvreg1", as.character(dir),
                              dims = as.integer(dims),               
                              Y = as.double(des$Y),
                              X = as.double(des$X),
                              indb = as.integer(des$indb),               
                              iter = as.integer(init$iter),
                              loglik = as.double(c(0, 0)),
                              mixture = as.double(init$mixture),
                              mixmoment = as.double(c(0, 0)),             
                              beta = as.double(init$beta),
                              b = as.double(init$b),
                              D = as.double(init$D),
                              r = as.integer(init$r),
                              Ys = as.double(init$y),       
                              otherp = as.double(init$otherp),
                              u = as.double(init$u),
                              prior.pari = as.integer(priordi$integer),
                              prior.pard = as.double(priordi$double),
                              revJump.pari = as.integer(revjumpdi$integer),
                              revJump.pard = as.double(revjumpdi$double),
                              prior.betai = as.integer(betadi$integer),
                              prior.betad = as.double(betadi$double),
                              prior.bi = as.integer(bdi$integer),
                              prior.bd = as.double(bdi$double),
                              nsimul = as.integer(nsimul.run3),
                              store = as.integer(storeV),
                              tolers = as.double(tolers),
                              err = integer(1),
             PACKAGE = thispackage)
   if (fit$err != 0) stop ("Something went wrong during the simulation.")
   cat("Simulation finished on                      ", date(), "   (iteration ", fit$iter, ")", "\n", sep = "")   

   toreturn <- fit$iter
   attr(toreturn, "call") <- call
   attr(toreturn, "prior") <- attr(priordi, "prior")
   attr(toreturn, "init") <- attr(priordi, "init")
   attr(toreturn, "prop.revjump") <- attr(revjumpdi, "prop.revjump")
   attr(toreturn, "prior.beta") <- attr(betadi, "prior.beta")
   attr(toreturn, "prior.b") <- attr(bdi, "prior.b")
   if (x) attr(toreturn, "x") <- des$X
   if (y) attr(toreturn, "y") <- des$Y
   class(toreturn) <- "bayessurvreg"
   
   return(toreturn)
}


bayessurvreg1.checkStore <- function(store)
{
  if (is.null(store$y)) store$y <- FALSE
  if (is.null(store$r)) store$r <- FALSE
  if (is.null(store$b)) store$b <- FALSE
  if (is.null(store$u)) store$u <- FALSE
  if (is.null(store$MHb)) store$MHb <- FALSE
  if (is.null(store$regresres)) store$regresres <- FALSE

  return(store)  
}  
#########################################################
#### AUTHOR:     Arnost Komarek                      ####
####             (2004)                              ####
####                                                 ####
#### FILE:       bayessurvreg1.files2init.R          ####
####                                                 ####
#### FUNCTIONS:  bayessurvreg1.files2init            ####
#########################################################

### ======================================
### bayessurvreg1.files2init
### ======================================
bayessurvreg1.files2init <- function(dir = getwd(), 
                                     row, 
                                     kmax)
{
  files <- dir(dir)    ## character vector with available files (except D and mixture files)
  fnames <- paste(c("iteration", "beta", "b", "Y", "r", "otherp", "u"), ".sim", sep = "")  
  compnames <- c("iter", "beta", "b", "y", "r", "otherp", "u")

  if (missing(row)){
    if (sum(!is.na(match(files, "iteration.sim")))){
      iter <- read.table(file = paste(dir, "/", "iteration.sim", sep = ""), header = TRUE, skip = 0)
      skip <- nrow(iter) - 1
    }
    else
      skip <- 0
  }
  else
    skip <- row - 1
  
  init <- list()

    ## Read everything except mixture and D matrix
    ##  (some files (r, Y, u, b) can be without header)
    ## * if there is no header, read the first line as initial values
  for (i in 1:length(fnames)){
    init[[compnames[i]]] <- NULL
    skip.now <- skip
    if (sum(!is.na(match(files, fnames[i])))){
      help <- read.table(file = paste(dir, "/", fnames[i], sep = ""), header = FALSE, nrows = 1, as.is = TRUE)
      if (is.character(help[1, 1]))
        header <- TRUE
      else{
        header <- FALSE
        skip.now <- 0
      }        
      help <- scan(file = paste(dir, "/", fnames[i], sep = ""), skip = skip.now + 1*header, nlines = 1)
      if (!length(help)) stop("Incorrect 'row' or 'header' parameter.")
      init[[compnames[i]]] <- as.numeric(help)
    }
  }

    ## Read D matrix
  init$D <- NULL
  if (sum(!is.na(match(files, "D.sim")))){
    header <- TRUE
    help <- scan(file = paste(dir, "/", "D.sim", sep = ""), skip = skip + 1*header, nlines = 1)
    if (!length(help)) stop("Incorrect 'row' or 'header' parameter.")
    init$D <- as.numeric(help[-1])     ## Remove determinant
  }
  
    ## Read mixture
  init$mixture <- NULL
  if (sum(!is.na(match(files, "mixmoment.sim"))) + sum(!is.na(match(files, "mweight.sim"))) +
      sum(!is.na(match(files, "mmean.sim"))) + sum(!is.na(match(files, "mvariance.sim"))) == 4){

    header <- TRUE
    
    if (header){
      mix <- read.table(file = paste(dir, "/mweight.sim", sep = ""), nrows = 1)
      kmax <- length(mix)
      mix <- read.table(paste(dir, "/mmean.sim", sep = ""), nrows = 1)
      kmax2 <- length(mix)
      if (kmax != kmax2) stop("Different kmax indicated by files mweight.sim and mmean.sim.")
      mix <- read.table(paste(dir, "/mvariance.sim", sep = ""), nrows = 1)     
      kmax2 <- length(mix)
      if (kmax != kmax2) stop("Different kmax indicated by files mweight.sim and mvariance.sim.")         
    }
    if (!header & missing(kmax)){
      stop("kmax must be given")
    }      
    
    kk <- scan(file = paste(dir, "/", "mixmoment.sim", sep = ""), skip = skip + 1*header, nlines = 1)[1]
    mweight <- scan(file = paste(dir, "/", "mweight.sim", sep = ""), skip = skip + 1*header, nlines = 1)
    mmean <- scan(file = paste(dir, "/", "mmean.sim", sep = ""), skip = skip + 1*header, nlines = 1)    
    mvariance <- scan(file = paste(dir, "/", "mvariance.sim", sep = ""), skip = skip + 1*header, nlines = 1)
    k.now <- length(mweight)
    if (k.now == 0) stop("Invalid mixture weights read.")
    k.now2 <- length(mmean)
    if (k.now != k.now2) stop("Different k indicated by files mweight.sim and mmean.sim.")
    k.now2 <- length(mvariance)
    if (k.now != k.now2) stop("Different k indicated by files mweight.sim and mmean.sim.")    

    init$mixture <- c(kk, mweight, rep(0, kmax-k.now), mmean, rep(0, kmax-k.now), mvariance, rep(0, kmax-k.now))
  }      

  return(init)
}

## Subfunction for bayessurvreg1.R
##  -> just to make it more readable
##
## Manipulation with initial values and prior specifications

bayessurvreg1.priorInit <- function(prior, init, Yinit, Xinit, n, nX, nrandom, ncluster, indb, randomInt, toler.chol){

   ## ============================================================================================
   ## Prior parameters (calculate these that were not given by the user and change the notation)
   ## part 1
   ## ============================================================================================
   prior.pari <- numeric(3)
   names(prior.pari) <- c("kmax", "k.prior", "Eb0.depend.mix")
   if (is.null(prior$kmax)) prior$kmax <- 5
   if (is.null(prior$k.prior)) prior$k.prior <- "poisson"
   if (is.null(prior$Eb0.depend.mix)) prior$Eb0.depend.mix <- FALSE
   if (is.null(prior$poisson.k)) prior$poisson.k <- 3
   prior.pari["k.prior"] <- pmatch(prior$k.prior, c("poisson", "uniform", "fixed"), nomatch = -1) - 1
               ## 0 = Poisson, 1 = Uniform, 2 = Fixed
   if (prior.pari["k.prior"] < 0) stop("Prior for k (number of mixture components) must be either poisson, uniform or fixed.")
   prior.pari["kmax"] <- prior$kmax
   prior.pari["Eb0.depend.mix"] <- 1*(prior$Eb0.depend.mix)

   prior.pard <- numeric(2*prior.pari["kmax"] + 7)
   names(prior.pard) <- c(paste("pi.split", 1:prior.pari["kmax"], sep = ""), paste("pi.birth", 1:prior.pari["kmax"], sep = ""),
                          "lambda", "delta", "xi", "kappa", "zeta", "g", "h")
   prior.pard["lambda"] <- prior$poisson.k


   ## =========================================================
   ## Get initial estimates
   ## =========================================================
   init.error <- "Something is wrong with your initials."
   fit.init <- survreg(Yinit ~ Xinit - 1, dist = "lognormal")          ## intercept is already included in Xinit

     ## index of the first iteration           
   if (is.null(init$iter) | is.na(init$iter)) init$iter <- 0                                 
   else                                       init$iter <- init$iter[1]
   
     ## initial mixture   
   if (is.null(init$mixture)){
     if (prior.pari["k.prior"] == 2) stop("init$mixture must be given when prior$k.prior is 'fixed'.")
     init$mixture <- numeric(1 + 3*prior$kmax)
     init$mixture[1] <- 1                                        ## initial k
     init$mixture[2] <- 1.0                                      ## initial weight of the first mixture component   
     init$mixture[2 + prior$kmax] <- fit.init$coefficients[1]    ## initial mean of the first mixture component
     init$mixture[2 + 2*prior$kmax] <- fit.init$scale^2          ## initial variance of the first mixture component
   }
   else{
     if (length(init$mixture) != 1 + 3*prior$kmax) stop("Incorrect init$mixture parameter supplied.")
     if (is.na(init$mixture[1])) stop("Incorrect init$mixture parameter supplied.")
     wi <- init$mixture[2:(1 + init$mixture[1])]
     mui <-init$mixture[(2 + prior$kmax):(1 + prior$kmax + init$mixture[1])]       
     sig2i <- init$mixture[(2 + 2*prior$kmax):(1 + 2*prior$kmax + init$mixture[1])]   
     if (sum(is.na(wi)) | sum(is.na(mui)) | sum(is.na(sig2i))) stop("Incorrect init$mixture parameter supplied.")
     if (sum(wi < 0)) stop("Incorrect init$mixture parameter supplied.")
     if (sum(sig2i <= 0)) stop("Incorrect init$mixture parameter supplied.")
     wi <- wi/sum(wi)                              ## to make sure that the sum is 1
     ordermu <- order(mui)
     k.temp <- init$mixture[1]
     init$mixture <- numeric(1 + 3*prior$kmax)     
     wi <- wi[ordermu]
     mui <- mui[ordermu]
     sig2i <- sig2i[ordermu]
     init$mixture[1] <- k.temp
     init$mixture[2:(1 + init$mixture[1])] <- wi
     init$mixture[(2 + prior$kmax):(1 + prior$kmax + init$mixture[1])] <- mui
     init$mixture[(2 + 2*prior$kmax):(1 + 2*prior$kmax + init$mixture[1])] <- sig2i     
   }

     ## initial beta parameter
   if (!nX) init$beta <- 0
   else{
     if (is.null(init$beta)){
       init$beta <- fit.init$coefficients[-1]                  ## remove the intercept
     }
     else{
       if (length(init$beta) < nX) stop("Incorrect init$beta parameter supplied.")
       init$beta <- init$beta[1:nX]
     }
     if (sum(is.na(init$beta))) stop("Incorrect init$beta parameter supplied.")
   }     
   
     ## initial values of the random effects
   if (!nrandom) init$b <- 0
   else{
     if (is.null(init$b)){
       bb <- fit.init$coefficients[-1][indb > 0]
       if (randomInt) bb <- c(0, bb)
       init$b <- rep(bb, ncluster)
     } 
     else{
       if (length(init$b) < nrandom*ncluster) stop("Incorrect init$b parameter supplied.")
       init$b <- init$b[1:(nrandom*ncluster)]
     }
     if (sum(is.na(init$b))) stop("Incorrect init$b parameter supplied.")
   }     
   
     ## initial values of the (transformed) latent response
   if (is.null(init$y)){
     init$y <- as.numeric(log(Yinit[,1]))          
   } 
   else{
     if (length(init$y) < n) stop("Incorrect init$y parameter supplied.")    
     init$y <- init$y[1:n]
   }
   if (sum(is.na(init$y))) stop("Incorrect init$y parameter supplied.")
   
     ## initial values of the component pertinences
   if (is.null(init$r)){
     init$r <- numeric(n) + 1         ## initially, everyone belongs to the first mixture component    
   }
   else{
     if (length(init$r) < n) stop("Incorrect init$r parameter supplied.")
     init$r <- init$r[1:n]
   } 
   if (sum(is.na(init$r)) | sum(init$r <= 0) | sum(init$r > init$mixture[1]))
     stop("Incorrect init$r parameter supplied.")
   
     ## initial values of the matrix D (covariance matrix of the random effects)
     ## (lower triangle of the matrix D in column major order)
   if (!nrandom) init$D <- 0   
   else{
     if (is.null(init$D)){
       init$D <- diag(nrandom)[lower.tri(diag(nrandom), diag = TRUE)]   ## identity matrix as initial D
     }
     else{
       if (length(init$D) < 0.5*nrandom*(1+nrandom)) stop("Incorrect init$D parameter supplied.")
       init$D <- init$D[1:(0.5*nrandom*(1+nrandom))]
     }
     if (sum(is.na(init$D))) stop("Incorrect init$D parameter supplied.")
   }     
   
     ## initial values of the remaining parameters (only eta at this moment)
   if (is.null(init$otherp)){
     init$otherp <- 1              ## initial of eta
   }
   else{
     init$otherp <- init$otherp[1]
   }
   if (sum(is.na(init$otherp))) stop("Incorrect init$otherp parameter supplied.")
   
     ## initial proposal vector for a split-combine move
   if (is.null(init$u)){
     init$u <- c(runif(1), 0, 0, runif(3*(prior$kmax - 1)))
   }
   else{
     if (length(init$u) < 3*prior$kmax) stop("Incorrect init$u parameter supplied.")  
     init$u <- init$u[1:(3*prior$kmax)]
   }     
   if (sum(is.na(init$u))) stop("Incorrect init$u parameter supplied.")  
   if (sum(init$u < 0 | init$u > 1)) stop("Incorrect init$u parameter supplied.")  

   
   ## =========================================================================================
   ## Prior parameters (calculate these that were not given by the user and change the notation)
   ## part 2
   ## =========================================================================================
   if (is.null(prior$dirichlet.w)) prior$dirichlet.w <- 1
   if (prior$dirichlet.w < 1) stop ("prior$dirichlet.w must be at least 1.")
   prior.pard["delta"] <- prior$dirichlet.w

   if (is.null(prior$mean.mu)) prior$mean.mu <- init$mixture[2 + prior.pari["kmax"]]    ## mean of the first comp.
   if (is.null(prior$var.mu)) prior$var.mu <- 2*init$mixture[2 + 2*prior.pari["kmax"]]  ## 2*variance of the first comp.
   if (prior$var.mu <= 0) stop("prior$var.mu must be positive.")
   prior.pard["xi"] <- prior$mean.mu
   prior.pard["kappa"] <- prior$var.mu

   if (is.null(prior$shape.invsig2)) prior$shape.invsig2 <- 1.5
   if (is.null(prior$shape.hyper.invsig2)) prior$shape.hyper.invsig2 <- 0.8
   if (is.null(prior$rate.hyper.invsig2)) prior$rate.hyper.invsig2 <- prior$var.mu
   if (prior$shape.invsig2 <= 0) stop("prior$shape.invsig2 must be positive.")
   if (prior$shape.hyper.invsig2 <= 0) stop("prior$shape.hyper.invsig2 must be positive.")
   if (prior$rate.hyper.invsig2 <= 0) stop("prior$rate.hyper.invsig2 must be positive.")      
   prior.pard["zeta"] <- prior$shape.invsig2
   prior.pard["g"] <- prior$shape.hyper.invsig2
   prior.pard["h"] <- prior$rate.hyper.invsig2

   if (is.null(prior$pi.split)) prior$pi.split <- c(1, rep(0.5, prior.pari["kmax"]-2), 0)
   if (is.null(prior$pi.birth)) prior$pi.birth <- c(1, rep(0.5, prior.pari["kmax"]-2), 0)
   if (prior$pi.split[1] != 1) stop("prior$pi.split[1] must be equal to 1.")
   if (prior$pi.birth[1] != 1) stop("prior$pi.birth[1] must be equal to 1.")
   if (prior$pi.split[prior.pari["kmax"]] != 0) stop("prior$pi.split[kmax] must be equal to 0.")
   if (prior$pi.birth[prior.pari["kmax"]] != 0) stop("prior$pi.birth[kmax] must be equal to 0.")      
   if (length(prior$pi.split) != prior.pari["kmax"]) stop("Incorrect length of a vector prior$pi.split.")
   if (length(prior$pi.birth) != prior.pari["kmax"]) stop("Incorrect length of a vector prior$pi.birth.")      
   prior.pard[paste("pi.split", 1:prior.pari["kmax"], sep = "")] <- prior$pi.split
   prior.pard[paste("pi.birth", 1:prior.pari["kmax"], sep = "")] <- prior$pi.birth
   
   priordi <- list(integer = prior.pari, double = prior.pard)
   attr(priordi, "init") <- init
   attr(priordi, "prior") <- prior
   
   return(priordi)
   
 }
## Subfunction for bayessurvreg1.R
##  -> just to make it more readable
##
## Manipulation with the specification of proposal jumps algorithms

bayessurvreg1.revjump <- function(prop.revjump)
{
  if (is.null(prop.revjump$algorithm)) prop.revjump$algorithm <- "basic"
  if (is.null(prop.revjump$moody.ring)) prop.revjump$moody.ring <- c(0.5, 0.5)
  if (is.null(prop.revjump$transform.split.combine)) prop.revjump$transform.split.combine <- "richardson.green"
  if (is.null(prop.revjump$transform.split.combine.parms)) prop.revjump$transform.split.combine.parms <-  c(2, 2, 2, 2, 1, 1)
  if (is.null(prop.revjump$transform.birth.death)) prop.revjump$transform.birth.death <- "richardson.green"


  ## Algorithm
  algorithm <- pmatch(prop.revjump$algorithm, table = c("basic", "independent.av", "correlated.av"), nomatch = 0, duplicates.ok = FALSE)[1]
  if (!algorithm) stop ("Unknown algorithm to generate canonical variables for reversible jumps.")
  algorithm <- algorithm - 1    ## NOW: 0 = basic, 1 = independent.av, 2 = correlated.av


  ## Moody ring parameters
  if (algorithm == 0){
    mrparm <- c(0.5, 0.5)
  }    
  else
    if (algorithm == 1){
      if(length(prop.revjump$moody.ring) < 1) stop("Incorrect prop.revjump$moody.ring parameter.")
      mrparm <- c(prop.revjump$moody.ring[1], 0.5)
    }
    else
      if (algorithm == 2){
        if(length(prop.revjump$moody.ring) < 2) stop("Incorrect prop.revjump$moody.ring parameter.")
        mrparm <- c(prop.revjump$moody.ring[1], 0.5)      
      }      
  names(mrparm) <- c("mr.epsilon", "mr.delta")
  if (sum(is.na(mrparm))) stop("Missing moody ring parameters.")
  if (sum(mrparm < 0 | mrparm > 0.5)) stop("Moody ring parameters must lie between 0 and 0.5.")


  ## Transformation for split-combine move
  transsc <- pmatch(prop.revjump$transform.split.combine,
                    table = c("richardson.green", "brooks", "identity"), nomatch = 0, duplicates.ok = FALSE)[1]
  if (!transsc) stop ("Unknown transformation for split-combine move.")
  transsc <- transsc - 1    ## NOW: 0 = richardson.green, 1 = brooks, 2 = identity 


  ## Split-combine transformation parameters
  if (transsc == 2){
    parmssc <- rep(1.0, 6)
  }
  else
    if (transsc == 0 || transsc == 1){
      if(length(prop.revjump$transform.split.combine.parms) < 6) stop("Incorrect prop.revjump$transform.split.combine.parms.")
      parmssc <- prop.revjump$transform.split.combine.parms[1:6]
    }      
  names(parmssc) <- paste("sc", 1:6, sep = "")
  if (sum(is.na(parmssc))) stop("Missing parameters for split-combine transformation.")
  if (sum(parmssc <= 0)) stop("Parameters for split-combine transformation must be positive.")


  ## Transformation for birth-death move
  transbd <- pmatch(prop.revjump$transform.birth.death,
                    table = c("richardson.green"), nomatch = 0, duplicates.ok = FALSE)[1]
  if (!transbd) stop ("Unknown transformation for birth.death move.")
  transbd <- transbd - 1    ## NOW: 0 = richardson.green

  ## Birth-death transformation parameters
  ## in this version: only allocate space of length 6
  parmsbd <- numeric(6)
  names(parmsbd) <- paste("bd", 1:6, sep = "")  

  revjumpi <- c(algorithm, transsc, transbd)
  names(revjumpi) <- c("algorithm", "transformation.sc", "transformation.bd")
  revjumpd <- c(mrparm, parmssc, parmsbd)

  rjdi <- list(integer = revjumpi, double = revjumpd)
  attr(rjdi, "prop.revjump") <- prop.revjump

  return(rjdi)    
}  
## Subfunction for bayessurvreg1.R
##  -> just to make it more readable
##
## Write headers to files where simulated values will be stored

bayessurvreg1.writeHeaders <- function(dir, prior, store, nX, X, names.random, ncluster, nrandom,
                                       rnamesX, unique.cluster, nBetaBlocks, nbBlocks)
{   
   sink(paste(dir, "/iteration.sim", sep = ""), append = FALSE)
   cat("iteration", "\n"); sink()

   sink(paste(dir, "/loglik.sim", sep = ""), append = FALSE)
   cat("loglik", "randomloglik", "\n", sep = "  "); sink()

   sink(paste(dir, "/mweight.sim", sep = ""), append = FALSE)
   cat(paste("w", 1:prior$kmax, sep = ""), "\n", sep = "      "); sink()

   sink(paste(dir, "/mmean.sim", sep = ""), append = FALSE)
   cat(paste("mu", 1:prior$kmax, sep = ""), "\n", sep = "      "); sink()

   sink(paste(dir, "/mvariance.sim", sep = ""), append = FALSE)
   cat(paste("sigma2", 1:prior$kmax, sep = ""), "\n", sep = "      "); sink()   
   
   sink(paste(dir, "/mixmoment.sim", sep = ""), append = FALSE)
   cat("       k", "              Intercept", "           Scale", "\n", sep = "  "); sink()
   
   if (nX){ sink(paste(dir, "/beta.sim", sep = ""), append = FALSE)
            cat(colnames(X), "\n", sep = "      "); sink() }
   else
     file.remove(paste(dir, "/beta.sim", sep = ""))
   

   if (store$b){ sink(paste(dir, "/b.sim", sep = ""), append = FALSE)
                 cat(paste(rep(names.random, ncluster), ".", rep(unique.cluster, rep(nrandom, ncluster)), sep = ""),
                     "\n", sep = "    "); sink() }
   else
     file.remove(paste(dir, "/b.sim", sep = ""))

   if (nrandom){ sink(paste(dir, "/D.sim", sep = ""), append = FALSE)
                 D <- diag(nrandom)
                 rows <- row(D)[lower.tri(row(D), diag = TRUE)]
                 cols <- col(D)[lower.tri(col(D), diag = TRUE)]            
                 cat("det", paste("D.", rows, ".", cols, sep = ""), "\n", sep = "      "); sink() }
   else
     file.remove(paste(dir, "/D.sim", sep = ""))

   if (store$y){sink(paste(dir, "/Y.sim", sep = ""), append = FALSE)
                cat(paste("Y", rnamesX, sep = ""), "\n", sep = "      "); sink() }
   else
     file.remove(paste(dir, "/Y.sim", sep = ""))   

   if (store$r){sink(paste(dir, "/r.sim", sep = ""), append = FALSE)
                cat(paste("r", rnamesX, sep = ""), "\n", sep = "      "); sink() }
   else
     file.remove(paste(dir, "/r.sim", sep = ""))   

   sink(paste(dir, "/otherp.sim", sep = ""), append = FALSE)
   cat("eta", "\n", sep = "      "); sink()

   sink(paste(dir, "/MHinfo.sim", sep = ""), append = FALSE)
   cat("accept.spl.comb", "split", "accept.birth.death", "birth  ", sep = "  ")
   if (nX > 0) cat(paste("beta.block.", 1:nBetaBlocks, sep = ""), "  ", sep = "  ")
   cat("\n", sep = "  ");
   sink()

   if (store$MHb & nrandom > 0){
     sink(paste(dir, "/MHbinfo.sim", sep = ""), append = FALSE)
     cat(paste("b.block.", rep(1:nbBlocks, ncluster), ".", rep(unique.cluster, rep(nbBlocks, ncluster)), "", sep = ""), sep = "  ")
     cat("\n", sep = "  ");
     sink()
   }
   else
     file.remove(paste(dir, "/MHbinfo.sim", sep = ""))     

   
   if (store$u){ sink(paste(dir, "/u.sim", sep = ""), append = FALSE)
                 headeru <- paste(" u.", rep(1:prior$kmax, rep(3, prior$kmax)), ".", rep(1:3, prior$kmax), sep = "")
                 headeru[1] <- "    mood"; headeru[2] <- " u.0.0"; headeru[3] <- " u.0.0"
                 cat(headeru, "\n", sep = "     ");
                 sink()
               }
   else
     file.remove(paste(dir, "/u.sim", sep = ""))

   if (store$regresres){sink(paste(dir, "/regresres.sim", sep = ""), append = FALSE)
                        cat(paste("res", rnamesX, sep = ""), "\n", sep = "      "); sink() }
   else
     file.remove(paste(dir, "/regresres.sim", sep = ""))      
   
}
#########################################################
#### AUTHOR:     Arnost Komarek                      ####
####             (05/05/2004)                        ####
####                                                 ####
#### FILE:       densplot2.R                         ####
####                                                 ####
#### FUNCTIONS:  densplot2                           ####
#########################################################
##
## Slightly modified 'densplot' function of 'coda' library
##
## ========================================================
##
## argument 'plot' added to indicate whether a plot is to be created directly
## or a list with data.frames for future plotting is to be returned
## and some other plotting arguments allowed to be changed by a user
##
densplot2 <-
function (x, plot = TRUE, show.obs = FALSE, bwf, bty = "n", main = "", xlim, ylim, xlab, ylab, ...) 
{
    xx <- as.matrix(x)
    toplot <- list()
    for (i in 1:nvar(x)) {
        y <- xx[, i, drop = TRUE]
        if (missing(bwf)) 
            bwf <- function(x) {
                x <- x[!is.na(as.vector(x))]
                return(1.06 * min(sd(x), IQR(x)/1.34) * length(x)^-0.2)
            }
        bw <- bwf(y)
        width <- 4 * bw
        if (max(abs(y - floor(y))) == 0 || bw == 0) 
            hist(y, prob = TRUE, main = main, ...)
        else {
            scale <- "open"
            if (max(y) <= 1 && 1 - max(y) < 2 * bw) {
                if (min(y) >= 0 && min(y) < 2 * bw) {
                  scale <- "proportion"
                  y <- c(y, -y, 2 - y)
                }
            }
            else if (min(y) >= 0 && min(y) < 2 * bw) {
                scale <- "positive"
                y <- c(y, -y)
            }
            else scale <- "open"
            dens <- density(y, width = width)
            if (scale == "proportion") {
                dens$y <- 3 * dens$y[dens$x >= 0 & dens$x <= 
                  1]
                dens$x <- dens$x[dens$x >= 0 & dens$x <= 1]
            }
            else if (scale == "positive") {
                dens$y <- 2 * dens$y[dens$x >= 0]
                dens$x <- dens$x[dens$x >= 0]
            }
            if (missing(ylim)) ylim <- c(0, max(dens$y))
            if (missing(xlab)) xlab <- paste("N =", niter(x), "  Bandwidth =", formatC(dens$bw))
            if (missing(ylab)) ylab <- ""
            if (missing(xlim)) xlim <- NULL
            if (plot){
                plot(dens, main = main, type = "l", bty = bty, xlab = xlab, ylab = ylab, xlim = xlim, ylim = ylim, ...)
               if (show.obs) 
                   lines(y[1:niter(x)], rep(max(dens$y)/100, niter(x)), 
                     type = "h")
            }                
            else
               toplot[[i]] <- data.frame(x = dens$x, y = dens$y)
        }
        if (!is.null(varnames(x)) && is.null(list(...)$main)) 
            title(paste("Density of", varnames(x)[i]))
    }
    if (plot) return(invisible(x))
    else      return(toplot)
}
#########################################################
#### AUTHOR:     Arnost Komarek                      ####
####             (2004)                              ####
####                                                 ####
#### FILE:       files2coda.R                        ####
####                                                 ####
#### FUNCTIONS:  files2coda                          ####
#########################################################
## 11/04/2004: rewritten such that it returns directly mcmc object
##             with chains for all parameters


### ======================================
### files2coda
### ======================================
files2coda <- function(files,
                       data.frames,
                       variant = 1,
                       dir = getwd(),
                       start = 1,
                       end,
                       thin = 1,
                       header = TRUE,
                       chain)
{
  filesindir <- dir(dir)    ## character vector with available files

  if (missing(files)) misfiles <- TRUE
  else                misfiles <- FALSE

  if (missing(data.frames)) misdatfr <- TRUE
  else                      misdatfr <- FALSE
  
  ## Appropriate files, if not given
  if (misfiles && misdatfr){
    if (variant == 1)
      files <- paste(c("mixmoment", "beta", "b", "D", "Y", "r", "otherp", "u", "MHinfo", "MHinfob", "loglik"), ".sim", sep = "")
    else
      stop ("Not yet implemented variant of bayessurvreg function.")
  }
  else{
    if (!misdatfr && misfiles) files <- character(0)
  }  

  ## Indeces of iterations
  iters <- NULL
  if (sum(!is.na(match(filesindir, "iteration.sim")))){
    help <- read.table(file = paste(dir, "/", "iteration.sim", sep = ""), header = header)
    if (!dim(help)[1]) stop("Incorrect 'header' parameter.")
    iters <- as.numeric(help[,1])
    if (length(iters) == 0) iters <- NULL
  }

  ## Start and end for mcmc function of coda library
  if (!is.null(iters)){
    if (start > length(iters)) stop("start is not compatible with iteration.sim file supplied.")
    mcstart <- iters[start]
    if (missing(end))
      end <- length(iters)
    else{
      if (end > length(iters))
        stop("end is not compatible with iteration.sim file supplied.")
      else
        if (end < start) stop("start and end are not compatible.")
    }      
    mcend <- iters[end]
  }
  else{
    mcstart <- start
    if (missing(end))
      end <- numeric(0)
    else
      if (end < start) stop("start and end are not compatible.")
    mcend <- end
  }    

  ## Create a big matrix with all sampled values
  tmc <- numeric(0)

  comp <- 1
  if (length(files) > 0){
    for (i in 1:length(files)){
      if (sum(!is.na(match(filesindir, files[i])))){
        if (files[i] == "mixture.sim"){
          help <- read.table(file = paste(dir, "/", files[i], sep = ""), header = header, nrow = 1)
          if (!dim(help)[1]) stop("Incorrect 'header' parameter.")
          kncol <- dim(help)[2]
          help <- scan(file = paste(dir, "/", files[i], sep = ""), skip = 1*header)
          wanna <- seq(1, length(help), by = kncol)
          help <- data.frame(k = help[wanna])
        }
        else{        
          help <- read.table(file = paste(dir, "/", files[i], sep = ""), header = header)
          if (!dim(help)[1]) stop("Incorrect 'header' parameter.")
        }          

        if (start > nrow(help)) stop("start is not compatible with data.")
        if (!length(end))
          end <- nrow(help)
        else{
         if (end > nrow(help))
            stop("end is not compatible with data.")
          else
            if (end < start) stop("start and end are not compatible.")
        }        
          
        wanna <- seq(start, end, by = thin)
        if (ncol(help) == 1){
          cname <- colnames(help)
          help <- as.data.frame(help[wanna, ])
          colnames(help) <- cname                
        }
        else{
          help <- help[wanna, ]
        }
        if (comp == 1) tmc <- help
        else           tmc <- cbind(tmc, help)        
        comp <- comp + 1
      }  
    }
  }

  if (!misdatfr){
    for (i in 1:length(data.frames)){
        if (missing(chain)) help <- get(data.frames[i])
        else                help <- get(data.frames[i])[[chain]]
        if (!dim(help)[1]) stop("Incorrect data.frame supplied.")

        if (start > nrow(help)) stop("start is not compatible with a data.frame.")
        if (!length(end))
          end <- nrow(help)
        else{
          if (end > nrow(help))
            stop("end is not compatible with a data.frame.")
          else
            if (end < start) stop("start and end are not compatible.")
        }
        wanna <- seq(start, end, by = thin)
        if (ncol(help) == 1){
          cname <- colnames(help)
          help <- as.data.frame(help[wanna, ])
          colnames(help) <- cname                
        }
        else{
          help <- help[wanna, ]
        }
        if (comp == 1) tmc <- help
        else           tmc <- cbind(tmc, help)        
        comp <- comp + 1
    }
  }

  ## Create an mcmc object
  if (length(mcend)) mc <- mcmc(tmc, start = mcstart, end = mcend, thin = thin)
  else               mc <- mcmc(tmc, start = mcstart, thin = thin)
  
  return(mc)  
}  


predictive <- function(
     formula,
     random,
     time0 = 0,
     data = parent.frame(),
     grid,
     type,
     subset,
     na.action = na.fail,
     quantile = c(0, 0.025, 0.5, 0.975, 1),                       
     nsimul = list(niter = 10, nwrite = 10),
     predict = list(Et=TRUE, t=FALSE, Surv=TRUE, hazard=FALSE, cum.hazard=FALSE),
     store = list(Et=TRUE, t = FALSE, Surv = FALSE, hazard = FALSE, cum.hazard=FALSE),
     Eb0.depend.mix = FALSE,
     dir = getwd(),
     toler.chol = 1e-10,
     toler.qr = 1e-10,
     ...)
{
   thispackage = "bayesSurv"
  
   transform = function(t){log(t)}
   dtransform = function(t){1/t}

   typeError<- pmatch(type, table = c("mixture", "spline", "polya.tree"), nomatch = 0) - 1
   if (typeError == -1 || typeError >= 2) stop("Unknown or not yet implemented error type.")
   
   control <- predictive.control(predict, store, quantile)
   predict <- control$predict
   store <- control$store
   predictShch <- predict$Surv || predict$hazard || predict$cum.hazard

   ## Starting time for the survival model
   ## ====================================
   if (missing(time0)) time0 <- 0
   if (time0 < 0) stop("time0 must be non-negative.")

   ## Extract all the design information from the function call
   ## ==========================================================
   m <- match.call(expand.dots = FALSE)
   des <- bayessurvreg.design(m, formula, random, data, transform, dtransform)

   ## Check whether needed files are available
   ## and whether at least first row has correct number of elements
   ## ==============================================================
   filesindir <- dir(dir)    ## character vector with available files   
   if (sum(!is.na(match(filesindir, "mweight.sim")))){
     mix <- read.table(paste(dir, "/mweight.sim", sep = ""), nrows = 1)     
     kmax <- length(mix)
   }
   else
     stop("File with simulated values of mixture weights not found.")

   if (sum(!is.na(match(filesindir, "mmean.sim")))){
     mix <- read.table(paste(dir, "/mmean.sim", sep = ""), nrows = 1)
     kmax2 <- length(mix)
     if (kmax != kmax2) stop("Different kmax indicated by files mweight.sim and mmean.sim.")
   }
   else
     stop("File with simulated values of mixture means not found.")
   
   if (sum(!is.na(match(filesindir, "mvariance.sim")))){
     mix <- read.table(paste(dir, "/mvariance.sim", sep = ""), nrows = 1)     
     kmax2 <- length(mix)
     if (kmax != kmax2) stop("Different kmax indicated by files mweight.sim and mvariance.sim.")
   }
   else
     stop("File with simulated values of mixture variances not found.")

   if (sum(!is.na(match(filesindir, "mixmoment.sim"))) == 0)
     stop("File mixmoment.sim not found.")     
   
   nbeta <- des$nfixed + des$nrandom - des$randomInt
   if (nbeta){
     if (sum(!is.na(match(filesindir, "beta.sim")))){
       beta <- scan(file = paste(dir, "/beta.sim", sep = ""), nlines = 1, skip = 1)
       if (length(beta) != nbeta) stop("Incorrect 'beta.sim' file supplied.")
     }
     else
       stop("File with simulated values of regression parameters not found.")
   }

   if (des$nrandom){
     nD <- 0.5*(des$nrandom*(1 + des$nrandom))
     if (sum(!is.na(match(filesindir, "D.sim")))){
       D <- scan(file = paste(dir, "/D.sim", sep = ""), nlines = 1, skip = 1)     
       if (length(D) != 1 + nD) stop("Incorrect 'D.sim' file supplied.")
     }
     else
       stop("File with simulated values of a covariance matrix of random effects not found.")
   }

   ## nsimul
   ## =======
   if (is.null(nsimul$nwrite)) nsimul$nwrite <- nsimul$niter
   if (nsimul$nwrite > nsimul$niter) nsimul$nwrite <- nsimul$niter
   
   ## Grids
   ## ======  
   if (predictShch){
     if (is.list(grid)){
       if (length(grid) != des$n) stop("Incorrect 'grid' parameter supplied.")
       ngrid <- sapply(grid, length)
       gridall <- unlist(grid)
     }
     else{
       ngrid <- rep(length(grid), des$n)
       gridall <- rep(grid, des$n)
     }
     if (!is.numeric(ngrid)) stop("Incorrect 'grid' parameter supplied.")
     if (sum(ngrid <= 0)) stop("Incorrect 'grid' parameter supplied.")
     if (sum(gridall <= 0.0)) stop("All grid values must be positive.")
   }
   cumngrid <- c(0, cumsum(ngrid))
  
   ## Create files to store predictive values
   ## =======================================
   filesindir <- dir(dir)    ## character vector with files in dir
      
   write.headers <- function(filename, predict0, dir, filesindir, n, gridall, cumngrid, obs = TRUE, label, namesX)
   {     
     lf <- nchar(filename)
     remove <- substring(filesindir, first = 1, last = lf)
     remove <- pmatch(remove, table = filename, nomatch = 0, duplicates.ok = TRUE)
     file.remove(paste(dir, "/", filesindir[remove == 1], sep = ""))
     if (predict0){
       if (obs){
         for (i in 1:n){
           sink(paste(dir, "/", filename, i, ".sim", sep = ""), append = FALSE)
           cat(paste(gridall[(cumngrid[i] + 1):cumngrid[i+1]], sep = ""), "\n", sep = "  "); sink()
         }
       }
       else{
           sink(paste(dir, "/", filename, ".sim", sep = ""), append = FALSE)
           cat(paste(label, namesX, sep = ""), "\n", sep = "     "); sink()
       } 
     }          
   }

   write.headers("predET", store$Et, dir, filesindir, des$n, gridall, cumngrid, FALSE, "ET", des$rnamesX)
   write.headers("quantET", predict$Et, dir, filesindir, des$n, gridall, cumngrid, FALSE, "ET", des$rnamesX)
   write.headers("predT", store$t, dir, filesindir, des$n, gridall, cumngrid, FALSE, "T", des$rnamesX)
   write.headers("quantT", predict$t, dir, filesindir, des$n, gridall, cumngrid, FALSE, "T", des$rnamesX)
   write.headers("predS", store$Surv, dir, filesindir, des$n, gridall, cumngrid)
   write.headers("predhazard", store$hazard, dir, filesindir, des$n, gridall, cumngrid)
   write.headers("predcumhazard", store$cum.hazard, dir, filesindir, des$n, gridall, cumngrid)
   write.headers("quantS", predict$Surv, dir, filesindir, des$n, gridall, cumngrid)
   write.headers("quanthazard", predict$hazard, dir, filesindir, des$n, gridall, cumngrid)
   write.headers("quantcumhazard", predict$cum.hazard, dir, filesindir, des$n, gridall, cumngrid)


   ## Sample
   ## ========
   dims <- c(des$n, des$ncluster, des$nwithin, des$nY, des$nX, des$nfixed, des$nrandom, 1*des$randomInt, nsimul$nwrite)
   dims2 <- c(length(quantile) ,ngrid)
   predictV <- c(predict$Et, predict$t, predict$Surv, predict$hazard, predict$cum.hazard)
   storeV <- c(store$Et, store$t, store$Surv, store$hazard, store$cum.hazard)
   nsimul <- c(nsimul$niter, 0, nsimul$nwrite)       ## (nthin is ignored)
   tolers <- c(toler.chol, toler.qr)
   prior.pari <- c(kmax, 0, 1*Eb0.depend.mix)

   cat("Simulation started on                       ", date(), "\n", sep = "")
   fit <- .C("predictive", as.integer(typeError),
                           as.character(dir),
                           as.integer(dims),
                           as.integer(dims2),
                           X = as.double(des$X),
                           indb = as.integer(des$indb),
                           quant = as.double(quantile),
                           grid = as.double(gridall),
                           prior.pari = as.integer(prior.pari),
                           prior.pard = as.double(time0),
                           nsimul = as.integer(nsimul),
                           predict = as.integer(predictV),
                           store = as.integer(storeV),
                           tolers = as.double(tolers),
                           err = integer(1),
             PACKAGE = thispackage)

   if (fit$err != 0) warning ("Something went wrong during the simulation.")
   cat("Simulation finished on                      ", date(), "\n", sep = "")

   return(fit$err)      
}
predictive.control <- function(predict, store, quantile)
{
  if (is.null(predict$Et)) predict$Et <- FALSE 
  if (is.null(predict$t)) predict$t <- FALSE
  if (is.null(predict$Surv)) predict$Surv <- FALSE
  if (is.null(predict$hazard)) predict$hazard <- FALSE
  if (is.null(predict$cum.hazard)) predict$cum.hazard <- FALSE
  if (!(predict$Et || predict$t || predict$Surv || predict$hazard || predict$cum.hazard))
    stop("Nothing to be predicted.")


  if (is.null(store$Et)) store$Et <- FALSE
  if (is.null(store$t)) store$t <- FALSE
  if (is.null(store$Surv)) store$Surv <- FALSE
  if (is.null(store$hazard)) store$hazard <- FALSE
  if (is.null(store$cum.hazard)) store$cum.hazard <- FALSE

  if (!predict$Et) store$Et <- FALSE
  if (!predict$t) store$t <- FALSE
  if (!predict$Surv) store$Surv <- FALSE
  if (!predict$hazard) store$hazard <- FALSE
  if (!predict$cum.hazard) store$cum.hazard <- FALSE

  if (sum(quantile < 0 | quantile > 1)) stop("Quantiles must lie between 0 and 1.")

  back <- list(predict = predict, store = store)
  return(back)
}

sampleCovMat <- function(sample)
{
  dimnames <- colnames(sample)
  sample <- as.matrix(sample)

  n.sample <- nrow(sample)

  ybar <- apply(sample, 2, mean)              ## mean
  t.y.min.ybar <- t(sample) - ybar            ## (sample - mean)'
  covmat <- t.y.min.ybar %*% t(t.y.min.ybar)  ## sum (y-ybar)*(y-ybar)'

  if (n.sample > 1)  covmat <- (1/(n.sample - 1)) * covmat
  ## else return only sum of squares, which is zero

  rownames(covmat) <- dimnames
  colnames(covmat) <- dimnames  
  
  return(covmat)  
}  
#########################################################
#### AUTHOR:     Arnost Komarek                      ####
####             (2004)                              ####
####                                                 ####
#### FILE:       traceplot2.R                        ####
####                                                 ####
#### FUNCTIONS:  traceplot2                          ####
#########################################################

### ======================================
### traceplot2
### ======================================
traceplot2 <- function(x, chains, bty = "n", main, xlab, ...)
{
  if (attr(x, "class") != "mcmc") stop("This function handles only objects of class mcmc.")
  
  if (missing(chains)){
    if (is.null(attr(x, "dim"))) chains <- 1
    else                         chains <- 1:attr(x, "dim")[2]
  }
  else{
    if (is.null(attr(x, "dim")))
      chains <- 1
    else{
      ch <- chains %in% 1:attr(x, "dim")[2]
      chains <- chains[ch]
    }      
  }

  iters <- attr(x, "mcpar")
  iters <- seq(iters[1], iters[2], by = iters[3])

  if (missing(xlab)) xlab <- "Iteration"
  
  for (i in 1:length(chains)){
    name <- attr(x, "dimnames")[[2]][chains[i]]
    if (is.null(name)) name <- paste("var ", i, sep = "")
    if (missing(main)) mmain <- paste("Trace of ", name, sep = "")
    else               mmain <- main
    if (is.null(attr(x, "dim"))) sample <- x
    else                         sample <- x[,chains[i]]
    
    plot(iters, sample, type = "l", lty = 1, xlab = xlab, bty = bty)
    title(main = mmain)    
  }

  return(invisible(x))    
}  
###############################################
#### AUTHOR:    Arnost Komarek             ####
####            (2004)                     ####
####                                       ####
#### FILE:      zzz.R                      ####
####                                       ####
#### FUNCTIONS: .First.lib                 ####
###############################################

### =============================================
### .First.lib
### =============================================
.First.lib <- function(lib, pkg)
{
   require(survival)
   require(coda)
   library.dynam("bayesSurv", pkg, lib)

   invisible()
}

