.packageName <- "bayesmix"
require(coda)
jags <- ifelse(.Platform$OS.type == "windows", "jags.exe", "jags")
jags <- system.file("exec", jags, package = "bayesmix")

if (file.exists(jags)) options(jags.exe = jags)

"BMMmodel" <-
function(y, k, priors, inits = "initsFS", aprioriWeights = 1, no.empty.classes = FALSE, restrict = "none", ...) {
  if (missing(y)) {
    call <- match.call(expand.dots = TRUE)
    model <- as.list(call[-1])
    model <- lapply(model, eval)
    class(model) <- "BMMsetup"
  }
  else {
    if (!is.null(dim(y))) {
      if (dim(y)[1] == 1) y <- y[1, ]
      else if (dim(y)[2] == 1) y <- y[ ,1]
      else stop("Only univariate data allowed")
    }
    y <- as.numeric(y)
    N <- as.numeric(length(y))
    
    if (!inherits(priors, "BMMpriors") & is.list(priors)) priors <- BMMpriors(priors, y = y)
    if (!inherits(priors, "BMMpriors")) stop("Priors not specified correctly")
    if (!all(sapply(priors$var, length) %in% c(0, 1, k))) stop("Priors not specified correctly - dimension differ")
    if (is.character(inits)) {
      inits <- get(inits)(y, k, restrict, ...)
    }
    for (i in names(inits)[names(inits) %in% names(priors$var)]) {
      if (length(inits[[i]]) != length(priors$var[[i]])) stop("Priors and inits dimension differ!")
    }
    if (length(aprioriWeights) != k) e <- rep(aprioriWeights[1], k) else e <- aprioriWeights
    model <- list()
    model$inits <- priors$inits <- inits
    if (priors$name[1] == "condconjugate") priors$var$B <- rep(NA, length(inits$tau*priors$var$B0inv))
    priors$var <- priors$var[!names(priors$var) %in% names(priors$inits)]
  
    index <- sapply(priors$var, function(x) (length(x) > 0) && !is.na(x))
    const <-  priors$var[index]
    const <- c(const, k = k, N = N)
    var <- c(priors$inits, priors$var[!index], list(e = e), list(y = y), list(S = y))  
    if (no.empty.classes) {
      var <- c(list(ind = matrix(0, nrow = N, ncol = k)),
               list(tot = vector(length = k)),
               list(seg = diag(k)),
               var)
      const <- c(const, list(Itot = rep(1, k)))
      model$inits$seg = diag(k)
    }
    varlist <- varSpec(c(const, var))
    
    bugs <- paste("var \n", varlist, sep = "")
    restrict <- match.arg(restrict, c("none", "mu", "tau"))
    if (restrict == "mu") {
      bugs <- paste(bugs,
                    "model\t{\n\tfor (i in 1:N) {\n\t\ty[i] ~ dnorm(mu,tau[S[i]]);\n\t\tS[i] ~ dcat(eta[]);\n\t}\n",
                    sep = "")
    }
    else if (restrict == "tau") {
      bugs <- paste(bugs,
                    "model\t{\n\tfor (i in 1:N) {\n\t\ty[i] ~ dnorm(mu[S[i]],tau);\n\t\tS[i] ~ dcat(eta[]);\n\t}\n",
                    sep = "")
    }
    else {
      bugs <- paste(bugs,
                    "model\t{\n\tfor (i in 1:N) {\n\t\ty[i] ~ dnorm(mu[S[i]],tau[S[i]]);\n",
                    "\t\tS[i] ~ dcat(eta[]);\n\t}\n", sep = "")
    }
    bugs <- paste(bugs, modelPriors(priors, restrict), sep = "")
    if (no.empty.classes) bugs <- paste(bugs, "\tfor (i in 1:N) {\n\t\tind[i,] <- seg[S[i],];\n\t}\n",
                                        "\tfor (j in 1:k){\n\t\ttot[j] <- sum(ind[,j]);\n",
                                        "\t\tItot[j] ~ dinterval(tot[j], 0);\n\t}\n", sep = "")
    bugs <- paste(bugs,paste("\teta[] ~ ddirch(e[]);\n}\n"), sep = "")
    model$data <- c(const, list(e = e), list(y = y))
    if (length(priors$name) > 1) {
      if (priors$name[2] == "tau") {
        if(!all(c(model$data$g0G0Half,model$data$g0Half)) & !("S0" %in% names(model$inits)))
          stop("Priors not specified correctly: Need an initial value for S0 with improper hierarchical prior.")
      }
      else stop("Should not be possible to have an hierarchical prior other than tau")
    }
    
    model$bugs <- bugs
    class(model) <- c("BMMmodel", "JAGSmodel")
  }
  model
}

"modelParameters" <-
function(priors) {
  parlist <- NULL
  for (i in names(priors$var)) {
    if (!any(is.na(priors$var[[i]]))) {
      cc <- priors$var[[i]]
      if (length(cc) > 1) {
        for (j in 1:length(cc)) {
          parlist <- paste(parlist, "\t", i,"[",j,"] <- ",cc[j],";\n ", sep = "")
        }
      }
    }
  }
  parlist
}

"modelPriors" <-
function(priors, restrict) {
  variants <- c("independence", "condconjugate")
  variant <- match.arg(tolower(priors$name[1]), variants)
  var <- c(priors$var, priors$inits)
  if (restrict == "tau") var$tau <- rep(NA, 1)
  if (restrict == "mu") var$mu <- rep(NA, 1)
  for (i in 1:length(var)) {
    if (length(var[[i]]) > 1) assign(names(var[i]), paste(names(var[i]),"[j]", sep = ""))
    else assign(names(var[i]), names(var[i]))
  }
  pr <- NULL
  if (variant == "independence") {
    pr <- c(pr, paste(mu, " ~ dnorm(", b0,",", B0inv, ");\n", sep = ""))
  }
  else if (variant == "condconjugate") {
    pr <- c(pr, paste(mu, " ~ dnorm(", b0,",", B, ");\n", sep = ""))
    pr <- c(pr, paste(B, " <- ", B0inv , "*", tau, ";\n", sep = ""))
  }
  pr <- c(pr, paste(tau ," ~ dgamma(", nu0Half, ",", nu0S0Half,");\n", sep = ""))

  if(length(priors$name) > 1 && priors$name[2] == "tau") {
    pr <- c(pr, paste(S0," ~ dgamma(", g0Half, ",",g0G0Half,");\n", sep = ""))
    pr <- c(pr, paste(nu0S0Half , " <- ", nu0Half ," * ", S0, ";\n", sep = ""))
  }

  priorsSpec <- paste("\tfor (j in 1:k) {\n\t\t",
                      paste(pr[grep("j", pr)], collapse = "\t\t"), "\t}\n\t",
                      paste(pr[-grep("j", pr)], collapse = "\t"), "\n",
                      sep = "")
  priorsSpec
}

"varSpec" <-
function(var)  {
  varlist = NULL
  for (i in names(var)) {
    if (!is.null(dim(var[[i]]))) {
      cc <- dim(var[[i]])
      varlist = paste(varlist, "\t", i, "[", cc[1],",", cc[2], "],\n", sep = "")
    }
    else {
      cc <- length(var[[i]])
      if (i == "S") {
        varlist = paste(varlist, "\t", i, "[", cc,"]; \n\n", sep = "")
      }
      else if (cc <= 1) {
        varlist = paste(varlist, "\t", i,",\n ", sep = "")
      }
      else{
        varlist = paste(varlist, "\t", i, "[",cc,"], \n", sep = "")
      }
    }
  }
}

print.JAGSmodel <- function(x, ...) {
  cat("Data for nodes: ", paste(names(x$data), collapse = ", "), "\n", sep ="")
  cat("Initial values for nodes: ", paste(names(x$inits), collapse = ", "), "\n\n", sep ="")
  cat("Model specification in BUGS language:\n\n")
  cat(x$bugs)
}

priorsRaftery <- function(y) {
  para <- list()
  para$b0 <- mean(y)
  R <- diff(range(y))
  para$B0 <- 2.6/R^2
  para$nu0 <- 2.56
  para$S0 <- (length(y)-1)/length(y)*var(y)
  para
}

priorsFish <- function(y) {
  para <- list()
  para$b0 <- median(y)
  para$B0 <- 10
  para$nu0 <- 20
  para$S0 <- 0
  para
}

priorsUncertain <- function(y) {
  para <- list()
  para$b0 <- mean(y)
  para$B0 <- Inf
  para$nu0 <- 0
  para$S0 <- 0
  para
}
  
BMMpriors <- function(specification, y) {
  priors <- list()
  default <- list(kind = "independence", parameter = "priorsUncertain", hierarchical = NULL, mod = list()) 
  if (missing(specification)) specification <- default
  else {
    n <- names(specification)
    s <- names(default)
    p <- pmatch(n, s)
    if(any(is.na(p)))
      stop(paste("\nInvalid name(s) in specification :", paste(n[is.na(p)], collapse=" ")))
    names(specification) <- s[p]
    for (i in names(specification)) {
      default[[i]] <- specification[[i]]
    }
    specification <- default
  }
  priors$name <- match.arg(tolower(specification$kind), c("independence", "condconjugate"))
  y <- as.vector(y)
  parameter <- specification$parameter
  if (is.character(parameter)) {
    specification$parameter <- get(parameter)(y)
  }
  if (length(specification$mod) > 0) {
    nam <- names(specification$mod)
    for (i in 1:length(specification$mod)) {
      specification$parameter[[nam[i]]] <- specification$mod[[i]]
    }
  }
  var <- specification$parameter
  priors$var <- list(b0 = var$b0, B0inv = 1/var$B0,
                     nu0Half = var$nu0/2, nu0S0Half = var$nu0*var$S0/2)
  if (!is.null(specification$hierarchical)) {
    specification$hierarchical <- match.arg(tolower(specification$hierarchical), c("tau"))
    if (specification$hierarchical == "tau") {
      priors$name <- c(priors$name, "tau")
      names(priors$name) <- c("type", "hierarchical prior for")
      for (x in c("g0", "G0")) {
        if (!x %in% names(var)) {
          var[[x]] <- 0
        }
      }
      priors$var <- c(priors$var, list(g0Half = var$g0/2, g0G0Half = var$g0/2*var$G0))
      priors$var$S0 <- rep(NA, max(length(var$S0), length(var$g0), length(var$G0)))
      priors$var$nu0S0Half <- rep(NA, max(length(priors$var$nu0S0Half), length(priors$var$S0)))
    }
    else stop("Hierarchical method not supported")
  }
  class(priors) <- c("BMMpriors", "JAGSpriors")
  priors
}

"initsPrint" <-
function(x) {
  x$B <- NULL
  var <- initsVar(x)
  paste("list(",paste(var, collapse = ",\n "),")\n")
}

"initsVar" <-
function(x) {
  n <- length(x)
  var <- vector(length = n)
  for (i in 1:n) var[i] <- paste(names(x)[i], " = c(",paste(x[[i]],collapse = ", "),")", sep = "")
  var
}

"initsFS" <-
function(x, k, restrict, initialValues = list()) {
  if (missing(restrict)) restrict <- ""
  x <- as.matrix(x)
  if (any(names(initialValues) %in% c("mu", "eta", "tau"))) stop("initialValues are not specified correctly")
  eta <- rep(1/k,k)
  eta[k] = 1-sum(eta[-k])
  if (restrict == "mu")  mu <- mean(x)
  else mu <- quantile(x, probs = seq(1/(k+1),k/(k+1),length = k))
  names(mu) <- NULL
  R <- IQR(x)
  sigma2 <- (R/1.34)^2
  if (restrict == "tau") tau <- 1/sigma2
  else tau <- rep(1/sigma2,k)
  z <- c(list(eta = eta, mu = mu, tau = tau), initialValues)
  z
}

JAGSsetup <- function(model, y, prefix, control, ...) {
  UseMethod("JAGSsetup")
}

JAGSsetup.default <- function(model, y, prefix, control, ...) {
  if (!inherits(model, "JAGSmodel")) stop("Only for use with 'JAGSmodel' objects!")
  with(model$data, dump(names(model$data), file = paste(prefix, "-data.R", sep = "")))
  with(model$inits, dump(names(model$inits), file = paste(prefix, "-inits.R", sep = "")))     
  
  if (length(model$bugs) > 1) {
    model$bugs <- .collapse(model$bugs, prefix)
  }
  write(model$bugs, file = paste(prefix,".bug", sep = ""))
  if (!any(names(control) %in% "text")) stop("control not specified correctly!")
  if (length(control$text) > 1) {
    control$text <- .collapse(control$text, prefix)
  }
  write(control$text, file = paste(prefix, ".cmd", sep = ""))
  return(list(control = control, model = model))
}

JAGSsetup.BMMsetup <- function(model, y, prefix, control, ...) {
  dummy <- model
  model <- list(k = 2, priors = BMMpriors(y = y), inits = "initsFS",
                aprioriWeights = 1, restrict = "none", no.empty.classes = FALSE)
  n <- names(dummy)
  s <- names(model)
  p <- pmatch(n, s)
  if(any(is.na(p)))
    stop(paste("\nInvalid name(s) in model :", paste(n[is.na(p)], collapse=" ")))
  names(dummy) <- s[p]
  for (i in names(dummy)) {
    model[[i]] <- dummy[[i]]
  }
  model <- BMMmodel(y, model$k, model$priors, model$inits,
                    model$aprioriWeights, model$no.empty.classes, model$restrict, ...)
  if (!inherits(model, "BMMmodel")) stop("Model not specified correctly")
  JAGSsetup(model, y, prefix, control)
}

JAGScall <- function(prefix, jags, quiet = FALSE) {
  if (is.null(jags)) jags = "jags"
  if (.Platform$OS.type == "windows") exit <- system(paste(jags, " ", prefix,".cmd", sep = ""))
  else  exit <- system(paste(jags, "< ",prefix,".cmd > /dev/null", sep = ""), ignore.stderr = quiet)
  if (exit) stop("System call not successfull")
  if (file.info("jags.out")[1] == 0) exit <- 1
  exit
}

JAGSread <- function(exit, transform = TRUE) {
  if (!exit) {
    if(!all(paste("jags" ,c("out", "ind"), sep = ".") %in% list.files()))
      stop("Cannot read jags output: .out or .ind file is missing!")
    results <- read.jags(quiet = TRUE)
    index <- grep("tau", colnames(results))
    variables <- unique(sapply(colnames(results), function(x) strsplit(x, "\\[")[[1]][1]))
    if (transform & length(index) > 0) {
      results[,index] <- 1/results[,index]
      colnames(results) <- sub("tau","sigma2", colnames(results))
      variables <- sub("tau", "sigma2", variables)
    }
  }
  else{
    results <- ifelse(file.exists("jags.dump"), source("jags.dump")[[1]], NULL)
    warning("Jags has encountered an error. Files are not deleted! Dump of jags will be returned.")
  }
  return(list(results = results, variables = variables))
}



"JAGScontrol" <-
function(variables, draw = 1000, burnIn = 0, seed) {
  text = NULL
  if (!missing(seed)) text <- paste("seed ", seed, "\n", sep = "")
  text <- paste(text,"model in \"", sep = "")
  text[2] <- ".bug\"\ndata in \""
  text[3] <- "-data.R\"\ncompile\ninits in \""
  text[4] <- "-inits.R\"\n\initialize\n"
  if (burnIn > 0) text[4] <- paste(text[4],"update ",burnIn,"\n", sep = "")
  text[4] <- paste(text[4], paste("monitor set ",variables,"\n", sep = "", collapse = ""), sep = "")
  text[4] <- paste(text[4],"update ",draw,"\ncoda *\nexit\n", sep = "")
  z <- list()
  z$text <- text
  z$variables <- variables
  class(z) <- "JAGScontrol"
  z
}

print.JAGScontrol <- function(x, prefix = "jags", ...) {
  cat("Commands for JAGS:\n\n")
  cat(paste(x$text, collapse = prefix))
}
"JAGSrun" <-  function(y, prefix = yname,  model = BMMmodel(k = 2),
                       control = JAGScontrol(variables = c("mu", "tau", "eta")),
                       tmp = TRUE, cleanup = TRUE, jags = getOption("jags.exe"), ...) {
  yname <- deparse(substitute(y))
  if (!is.null(dim(y))) {
    if (dim(y)[1] == 1) y <- y[1,]
    else if (dim(y)[2] == 1) y <- y[,1]
    else stop("Only univariate response allowed")
  }
  y <- as.numeric(y)
  cl <- match.call()
  if (tmp) {
    dir <- getwd()
    tmpdir <- tempdir()
    if (!file.exists(tmpdir)){
      if (!dir.create(tmpdir)) stop("Error creating tmp directory")
    }
    setwd(tmpdir)
  }
  specification <- JAGSsetup(model, y, prefix, control, ...)
  exit <- JAGScall(prefix, jags)
  results <- JAGSread(exit)
  if (any(!is.finite(results$results))) {
    warning("Infinite values occured: These draws are omitted!")
    results$results <- as.mcmc(na.omit(results$results))
    if (dim(results$results)[1] == 0) results$results <- NULL
  }
  if (!exit) {
    if (cleanup) {
     unlink(c(paste(prefix, c(".cmd", ".bug","-inits.R", "-data.R", ".txt"), sep = ""), "jags.out", "jags.ind"))
    }
    if (tmp) setwd(dir)
  }
  z = list(call = cl, results = results$results, model = specification$model,
    variables = results$variables, data = y)
  class(z) <- "jags"
  z
}

"summaryShort.mcmc" <-
function (object, quantiles = c(0.025, 0.975), 
    ...) 
{
    x <- as.mcmc(object)
    statnames <- c("Mean", "SD")
    varstats <- matrix(nrow = nvar(x), ncol = length(statnames), 
        dimnames = list(varnames(x), statnames))
    if (is.matrix(x)) {
        xmean <- apply(x, 2, mean)
        xvar <- apply(x, 2, var)
        varquant <- t(apply(x, 2, quantile, quantiles))
    }
    else {
        xmean <- mean(x, na.rm = TRUE)
        xvar <- var(x, na.rm = TRUE)
        varquant <- quantile(x, quantiles)
    }
    varstats[, 1] <- xmean
    varstats[, 2] <- sqrt(xvar)
    varstats <- drop(varstats)
    varquant <- drop(varquant)
    out <- list(statistics = varstats, quantiles = varquant,
        start = start(x), end = end(x), thin = thin(x), nchain = 1)
    class(out) <- "summaryShort.mcmc"
    return(out)
}

"print.jags" <-
function(x, ...) {
  cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
  if (inherits(x$model, "BMMmodel")) {
    if (is.null(x$results)) cat("No results!\n")
    else {
      cat("Markov Chain Monte Carlo (MCMC) output:\nStart =", start(x$results), 
          "\nEnd =", end(x$results), "\nThinning interval =", thin(x$results), 
          "\n")
      for (i in x$variables) {
        y <- x$results[,grep(i, colnames(x$results)), drop = FALSE]
        if(dim(y)[2] <=  x$model$data$k) {
          yout <- summaryShort.mcmc(y)
          class(yout) <- "summaryShort.mcmc"
          cat(paste("\n Empirical mean, standard deviation and 95% CI for", i, "\n"))
          print(yout, ...)
        }
      }
    }
  }
}

"print.summaryShort.mcmc" <-
function(x, digits = max(3, .Options$digits - 3), ...) {
  if (is.matrix(x$statistics)) {
    print(cbind(x$statistics, x$quantiles), digits = digits, ...)
  }
  else print(c(x$statistics, x$quantiles), digits = digits, ...)
}

.collapse <- function(text, prefix) {
  paste(text, collapse = prefix)
}
## Function taken from e1071
"permutations" <- function (n) {
    if (n == 1) 
        return(matrix(1))
    else if (n < 2) 
        stop("n must be a positive integer")
    z <- matrix(1)
    for (i in 2:n) {
        x <- cbind(z, i)
        a <- c(1:i, 1:(i - 1))
        z <- matrix(0, ncol = ncol(x), nrow = i * nrow(x))
        z[1:nrow(x), ] <- x
        for (j in 2:i - 1) {
            z[j * nrow(x) + 1:nrow(x), ] <- x[, a[1:i + j]]
        }
    }
    dimnames(z) <- NULL
    z
}


"Sort" <- function(x, by = NULL) {
  if (!(inherits(x, "jags") && inherits(x$model, "BMMmodel"))) 
    stop("Use only with 'jags' objects with model of class 'BMMmodel'.")
  x.old <- x
  n <- dim(x$results)
  if (is.null(by)) by <- x$variables
  else by <- x$variables[pmatch(by, x$variables)]
  by <- by[1]
  if (is.na(by)) stop("by not specified correctly")
  index <- grep(by, colnames(x$results))
  nn <- length(index)
  if (nn != x$model$data$k) stop("by not specified correctly")
  dd <- order(row(x$results[,index]), x$results[,index])
  ind <- apply(x$results[, index], 1, order)
  for (name in x$variables) {
    ii <- grep(name, colnames(x$results))
    if (length(ii) == nn) {
      x$results[,ii] <- matrix(x$results[,ii][dd], nrow = n[1], byrow = TRUE)
    }
    else if (length(levels(as.factor(x$results[,ii]))) == x$model$data$k) {
      ps <- permutations(x$model$data$k)
      for (j in 1:dim(ps)[1]) {
        ps1 <- ps[j,]
        index <- apply(ind, 2, function(x) all(x == ps1))
        if (any(index)) {
          dummy <- factor(x$results[index,ii], levels = 1:x$model$data$k)
          levels(dummy) <- order(ps1)
          x$results[index,ii] <- as.numeric(levels(dummy))[as.integer(dummy)]
        }
      }
    }
    else if (length(ii) != 1) {
      warning("Sorting not successful. Original object returned!")
      return(x.old)
    }
  }
  x
}

haveJAGS <- function(jags = getOption("jags.exe")) {
  opt <- options("warn" = -1)
  on.exit(options(opt))
  dir <- getwd()
  tmpdir <- tempdir()
  if (!file.exists(tmpdir)){
    if (!dir.create(tmpdir)) stop("Error creating tmp directory")
  }
  setwd(tmpdir)
  if (file.exists("haveJAGS.cmd")) stop("Remove file 'haveJAGS.cmd' first.")
  write("exit", file = "haveJAGS.cmd")
  if (is.null(jags) || jags == "") jags <- "jags"
  if (.Platform$OS.type == "windows") exit <- system(paste(jags, "haveJAGS.cmd"))
  else  exit <- system(paste(jags, " < haveJAGS.cmd > /dev/null", sep = ""), ignore.stderr = TRUE)
  unlink("haveJAGS.cmd")
  setwd(dir)
  as.logical(!exit)
}
"BMMdiag" <-
function(object, which = 1:2, variables, ask = interactive(), fct1, fct2, 
         xlim, ylim, auto.layout = TRUE, caption = NULL, main = "", ...) {  
  if (!(inherits(object, "jags") && inherits(object$model, "BMMmodel"))) 
    stop("Use only with 'jags' objects with model of class 'BMMmodel'.")
  if (!is.numeric(which) || any(which < 1) || any(which > 2)) 
    stop("`which' must be in 1:2")
  k <- object$model$data$k
  oldpar <- NULL
  on.exit(par(oldpar))
  oldpar <- par(ask = ask)
  show <- rep(FALSE, 2)
  show[which] <- TRUE
  if (missing(variables)) variables <- object$variables
  setxlim <- ifelse(missing(xlim), TRUE, FALSE)
  setylim <- ifelse(missing(ylim), TRUE, FALSE)
  vars <- variables[sapply(variables, function(x) length(grep(x, colnames(object$results))) <= k)]
  numVars <- length(vars)
  if (is.null(caption)) {
    caption <- c(sapply(vars, function(x) paste(x, "[k] versus ", vars, "[k]", sep = ""))
                 [lower.tri(matrix(nrow = numVars, ncol = numVars))],
                 paste(vars, "[k] versus ", vars ,"[l]", sep = ""))
  }
  if (show[1]) {
    if (numVars > 1) {
      if (auto.layout) oldpar <- c(oldpar, par(mfrow = c(1, numVars*(numVars-1)/2)))
      h <- 0
      for (i in 1:(numVars-1)) {
        kvar1 <- grep(vars[i], colnames(object$results))
        var1 <- matrix(object$result[,kvar1], ncol = length(kvar1))
        vars1 <- vars[i]
        if (!missing(fct1)) {
          var1 <- get(fct1)(var1)
          vars1 <- paste(fct1, "(", vars1, ")", sep = "")
        }
        if (setxlim) xlim <- range(var1)
        for (j in (i + 1):numVars) {
          h <- h + 1
          kvar2 <- grep(vars[j], colnames(object$results))
          if (length(kvar2) <= k) {
            vars2 <- vars[j]
            var2 <- matrix(object$result[,kvar2], ncol = length(kvar2))
            if (!missing(fct2)) {
              var2 <- get(fct2)(var2)
              vars2 <- paste(fct2, "(", vars2, ")", sep = "")
            }            
            if (setylim) ylim <- range(var2)
            plot(var1[,1], var2[,1], xlim = xlim,
                 ylim = ylim, xlab = vars1, ylab = vars2, main = main, ...)
            mtext(caption[h], 3, 0.25)
            h <- max(length(kvar1), length(kvar2))
            if (h > 1) {
              for (l in 2:h) points(var1[, min(l,length(kvar1))], var2[, min(l, length(kvar2))], ...)
            }
          }
        }
      }
    }
    else warning("The first plot option requires at least two variables.")
  }
  if (show[2]) {
    if (auto.layout) oldpar <- c(oldpar, par(mfrow = c(1, numVars)))
    for (l in 1:numVars) {
      varNam <- vars[l]
      kvar <- grep(varNam, colnames(object$results))
      if (length(kvar) > 1) {
        var <- matrix(object$result[,kvar], ncol = length(kvar))
        if (setxlim) xlim <- range(var)
        plot(xlim, xlim, type = "l", xlab = paste(varNam,"[k]", sep = ""), ylab = paste(varNam, "[l]", sep = ""),
             main = main, ...)
        mtext(caption[numVars*(numVars-1)/2+l], 3, 0.25)
        for (i in 1:(length(kvar)-1)) {
          for (j in (i+1):length(kvar)) {
            points(var[,i], var[,j], ...)
            points(var[,j], var[,i], ...)
          }
        }
      }
    }
  }
}

BMMposteriori <- function(object, class, caption = NULL, plot = TRUE, auto.layout = TRUE, ...) {
  if (!(inherits(object, "jags") && inherits(object$model, "BMMmodel"))) 
    stop("Use only with 'jags' objects with model of class 'BMMmodel'.")
  k <- object$model$data$k
  if (missing(class)) class <- 1:k
  if (is.null(caption)) caption <- paste("Group", class)
  S <- object$result[,grep("S", colnames(object$results))]
  if (dim(S)[2] == 0) stop("A posteriori plot not possible. Please provide class observations!")
  uniqPoints <- unique(object$data)
  n <- dim(object$results)[1]
  tab <- sapply(uniqPoints, function(x)
                table(factor(S[,object$data == x], levels = 1:k))/(n*length(which(object$data == x))))
  x <- list()
  x$post <- tab[class,]
  x$data <- uniqPoints
  class(x) <- "BMMposteriori"
  if (plot) {
    if (auto.layout) {
      oldpar <- par(mfrow = c(length(class), 1))
      on.exit(par(oldpar))
    }
    plot(x, caption, ...)
    invisible(x)
  }
  else x
}  

plot.BMMposteriori <- function(x, caption, main = "", ...) {
  if (!is.matrix(x$post)) x$post <- matrix(x$post, nrow = 1)
  for (i in 1:dim(x$post)[1]) {
    plot(x$data, x$post[i,], type = "h", xlab = "data", ylab = "a posteriori probability",
         ylim = c(0,1), main = main, ...)
    mtext(caption[i], 3, 0.25)
    points(x$data, x$post[i,], pch = 19, ...)
  }
}

"plot.jags" <-
function (x, variables = NULL, trace = TRUE, density = TRUE, 
                       smooth = TRUE, bwf, num, xlim, auto.layout = TRUE, ask = interactive(), ...)  
{
  if (inherits(x$model, "BMMmodel")) {
    if (is.null(variables)) {
      variables <- x$variables
      variables <- variables[sapply(variables, function(y) length(grep(y, colnames(x$results))) <= x$model$data$k)]
    }
    for (name in variables) {
      u <- x$results[,grep(name,colnames(x$results)), drop = FALSE]
      if (NCOL(u) > 0) {
        if (!missing(num)) {
          if (any(num > NCOL(u))) warning("num modified.")
          num <- num[num <= NCOL(u)]
          u <- u[,num, drop = FALSE]
        }
        oldpar <- NULL
        on.exit(par(oldpar))
        if (auto.layout) {
          mfrow <- set.mfrow(Nchains = nchain(u), Nparms = nvar(u), 
                             nplots = trace + density)
          oldpar <- par(mfrow = mfrow)
        }
        oldpar <- c(oldpar, par(ask = ask))
        for (i in 1:nvar(u)) {
          y <- as.matrix(u)[, i, drop = FALSE]
          if (trace) 
            traceplot(y, smooth = smooth)
          if (density) {
            if (missing(xlim)) xl <- range(u)
            else xl <- xlim
            if (missing(bwf)) 
              densplot(y, xlim = xl, ...)
            else densplot(y, bwf = bwf, xlim = xl, ...)
          }
        }
      }
      else warning("Variable ", name, " omitted.")
    }
  }
  else plot(x$results)
}
"randomPermutation" <- function(x) {
  if (!(inherits(x, "jags") && inherits(x$model, "BMMmodel"))) 
    stop("Use only with 'jags' objects with model of class 'BMMmodel'.")
  k <- x$model$data$k
  n <- dim(x$results)
  permutedIndex <- as.vector(t(apply(matrix(1:(n[1]*k), ncol = k), 1, sample, size = k)))
  variables <- x$variables
  dropIndex <- NULL
  for (i in 1:length(variables)) {
    name <- variables[i]
    ii <- grep(name, colnames(x$results))
    if (length(ii) == k) {
      dummy <- x$results[,ii]
      dummy <- dummy[permutedIndex]
      x$results[,ii] <- dummy
    }
    else if (length(ii) > 1) {
      x$results <- x$results[,-ii]
      dropIndex <- c(dropIndex, i)
    }
  }
  if (length(dropIndex)) {
    x$variables <- x$variables[-dropIndex]
    warning("Variables have been dropped.")
  }
  x
}
