.packageName <- "geoRglm"

"prior.glm.control" <- 
  function(beta.prior = c("flat", "normal", "fixed"), beta = NULL, beta.var.std = NULL, 
           sigmasq.prior = c("uniform", "sc.inv.chisq", "reciprocal", "fixed"), sigmasq = NULL, df.sigmasq = NULL,
           phi.prior = c("uniform","exponential", "fixed", "squared.reciprocal", "reciprocal"), phi = NULL, 
           phi.discrete = NULL, tausq.rel = 0)
{
  beta.prior <- match.arg(beta.prior)
  if(beta.prior == "fixed" & is.null(beta))
    stop("argument \"beta\" must be provided with fixed value for this parameter")
  if(beta.prior == "normal"){
    if(is.null(beta) | is.null(beta.var.std))
      stop("arguments \"beta\" and \"beta.var.std\" must be provided when using normal prior for the parameter beta")
    if((length(beta))^2 != length(beta.var.std))
      stop(" beta and beta.var.std have incompatible dimensions")
    if(any(beta.var.std != t(beta.var.std)))
      stop(" non symmetric matrix in beta.var.std")
    if(inherits(try(chol(beta.var.std)), "try-error"))
      stop(" matrix in beta.var.std is not positive definit")
  }
  ##
  sigmasq.prior <- match.arg(sigmasq.prior)
  if(sigmasq.prior == "fixed" & is.null(sigmasq))
    stop("argument \"sigmasq\" must be provided when the parameter sigmaq is fixed")
  if(sigmasq.prior == "sc.inv.chisq")
    if(is.null(sigmasq) | is.null(df.sigmasq))
      stop("arguments \"sigmasq\" and \"df.sigmasq\" must be provided for inverse chisq prior")
  if(!is.null(sigmasq))
    if(sigmasq < 0) stop("negative values not allowed for \"sigmasq\"")
  if(sigmasq.prior == "reciprocal"){
    warning("This choice of sigmasq.prior gives an improper posterior !!!!!!! \n")
    sigmasq <- 0
    df.sigmasq <- 0
  }
  if(sigmasq.prior == "uniform"){
    sigmasq <- 0
    df.sigmasq <- -2
  }
  ##
  if(!is.null(phi) && length(phi) > 1)
    stop("prior.glm.control: length of phi must be one. ")
  if(is.numeric(phi.prior)){
    phi.prior.probs <- phi.prior
    phi.prior <- "user"
    if(is.null(phi.discrete))
      stop("prior.glm.control: argument phi.discrete with support points for phi must be provided\n")
    if(length(phi.prior.probs) != length(phi.discrete))
      stop("prior.glm.control: user provided phi.prior and phi.discrete have incompatible dimensions\n")
    if(round(sum(phi.prior.probs), dig=6) != 1)
      stop("prior.glm.control: prior probabilities provided for phi do not sum up to 1")
  }
  else phi.prior <- match.arg(phi.prior)
  if(phi.prior == "fixed"){
    if(is.null(phi)){
      stop("argument \"phi\" must be provided with fixed prior for this parameter")
    }
    phi.discrete <- phi
  }
  else{
    if(phi.prior == "exponential" & (is.null(phi) | (length(phi) > 1)))
      stop("argument \"phi\" must be provided when using the exponential prior for the parameter phi")
    if(any(phi.prior == c("reciprocal", "squared.reciprocal")) & any(phi.discrete == 0)){
      warning("degenerated prior at phi = 0. Excluding value phi.discrete[1] = 0")
      phi.discrete <- phi.discrete[phi.discrete > 1e-12]
    }
    if(!is.null(phi.discrete)){
      discrete.diff <- diff(phi.discrete)
      if(round(max(1e08 * discrete.diff)) != round(min(1e08 * discrete.diff)))
        stop("The current implementation requires equally spaced values in the argument \"phi.discrete\"\n")
    }
    if(phi.prior != "exponential") phi <- NULL
    if(sigmasq.prior == "fixed") stop("option for fixed sigmasq and random phi not implemented")
  }
  if(any(phi.discrete < 0))
    stop("negative values not allowed for parameter phi")
  ##
  if(is.null(tausq.rel)) stop("argument \"tausq.rel\" must be provided")
  ##
  ip <- list(beta=list(), sigmasq=list(), phi=list())
  ##
  if(beta.prior == "fixed"){
    ip$beta$status <- "fixed"
    ip$beta$fixed.value <- beta 
  }
  else{
    ip$beta <- list(dist = beta.prior)
    if(beta.prior == "flat")
      ip$beta$pars <- c(0, +Inf)
    if(beta.prior == "normal"){
      if(length(beta) == 1)
        ip$beta$pars <- c(mean=beta, var.std=beta.var.std)
      else
        ip$beta$pars <- list(mean=beta, var.std=beta.var.std)
    }
  }
  ##
  if(sigmasq.prior == "fixed"){
    ip$sigmasq$status <- "fixed"
    ip$sigmasq$fixed.value <- sigmasq 
  }
  else{
    ip$sigmasq <- list(dist = sigmasq.prior)
    if(sigmasq.prior == "reciprocal")
      ip$sigmasq$pars <- c(df=0, var=+Inf)
    if(sigmasq.prior == "uniform")
      ip$sigmasq$pars <- c(df=-2, var=+Inf)
    if(sigmasq.prior == "sc.inv.chisq")
      ip$sigmasq$pars <- c(df=df.sigmasq, var=sigmasq)
  }
  ##
  if(phi.prior == "fixed"){
    ip$phi$status <- "fixed"
    ip$phi$fixed.value <- phi
  }
  else{
    ip$phi$dist <- phi.prior
    if(is.null(phi.discrete))
      stop("phi.discrete must be given when parameter phi is random")
    else{
      ip$phi$probs <- switch(phi.prior,
                             uniform = rep(1/length(phi.discrete), length(phi.discrete)),
                             exponential = (1/phi) * exp(- phi.discrete/phi),
                             squared.reciprocal = (1/(phi.discrete^2))/sum(1/(phi.discrete^2)),
                             reciprocal = (1/phi.discrete)/sum(1/phi.discrete),
                             user = phi.prior.probs)
      names(ip$phi$probs) <- phi.discrete
    }
    if(phi.prior == "exponential") ip$phi$pars <- c(ip$phi$pars, exp.par=phi)
  }
  ##
  ip$tausq.rel <- list(status = "fixed", fixed.value = tausq.rel)
  ##
  res <- list(beta.prior = beta.prior, beta = beta, beta.var.std = beta.var.std, sigmasq.prior = sigmasq.prior, sigmasq = sigmasq, 
              df.sigmasq = df.sigmasq, phi.prior = phi.prior, phi = phi, 
              phi.discrete = phi.discrete, tausq.rel = tausq.rel, priors.info = ip)
  class(res) <- "prior.geoRglm"
  return(res)
}


"prior.glm.check.aux" <-
  function(prior, fct)
{
  if(class(prior) != "prior.geoRglm"){
    if(!is.list(prior))
      stop(paste(fct,": argument prior only takes a list or an output of the function prior.glm.control"))
    else{
      prior.names <- c("beta.prior", "beta", "beta.var.std", "sigmasq.prior",
                       "sigmasq", "df.sigmasq", "phi.prior", "phi", "phi.discrete", "tausq.rel") 
      prior <- object.match.names(prior,prior.names)
      if(is.null(prior$beta.prior)) prior$beta.prior <- "flat"
      if(is.null(prior$sigmasq.prior)) prior$sigmasq.prior <- "uniform"
      if(is.null(prior$phi.prior)) prior$phi.prior <- "uniform"
      if(is.null(prior$tausq.rel)) prior$tausq.rel <- 0
      prior <- prior.glm.control(beta.prior = prior$beta.prior,
                                 beta = prior$beta, beta.var.std = prior$beta.var.std,
                                 sigmasq.prior = prior$sigmasq.prior,
                                 sigmasq = prior$sigmasq,  df.sigmasq = prior$df.sigmasq,
                                 phi.prior = prior$phi.prior,
                                 phi = prior$phi, phi.discrete = prior$phi.discrete, 
                                 tausq.rel = prior$tausq.rel)
    }
  }
  return(prior)
}


"image.glm.krige.bayes" <-
  function (x, locations, borders, 
            values.to.plot = c("median", "uncertainty",
              "quantiles", "probabilities", "simulation"),
            number.col, coords.data, xlim, ylim,
            x.leg, y.leg, ...) 
{
  ldots <- match.call(expand.dots = FALSE)$...
  ldots[!is.na(match(names(ldots), "offset.leg"))] <- NULL
  if(is.null(ldots[!is.na(match(names(ldots), "xlab"))])) ldots$xlab <- "X Coord"
  if(is.null(ldots[!is.na(match(names(ldots), "ylab"))])) ldots$ylab <- "Y Coord"
  if(missing(x)) x <- NULL
  attach(x)
  on.exit(detach(x))
  if(missing(locations))
    locations <-  eval(attr(x, "prediction.locations"))
  if(is.null(locations)) stop("prediction locations must be provided")
  if(ncol(locations) != 2)
    stop("locations must be a matrix or data-frame with two columns")
  if(!is.numeric(values.to.plot))
    values.to.plot <-
      match.arg(values.to.plot,
                choices = c("median", "uncertainty",
                  "quantiles", "probabilities", "simulation"))


  if(missing(borders)){
    if(!is.null(attr(x, "borders"))) borders.arg <- borders <- eval(attr(x, "borders"))
    else borders.arg <- borders <- NULL
  }
  else{
    borders.arg <- borders
    if(is.null(borders)) borders <- eval(attr(x, "borders"))
  }
  
  if(missing(borders)){
    if(!is.null(attr(x, "borders"))) borders <- eval(attr(x, "borders"))
    else borders <- NULL
  }
  if(missing(number.col)) number.col <- NULL
  if(missing(coords.data)) coords.data <- NULL
  if(missing(xlim)) xlim <- NULL
  if(missing(ylim)) ylim <- NULL
  if(missing(x.leg)) x.leg <- NULL
  if(missing(y.leg)) y.leg <- NULL
  if(!is.null(attr(x, 'sp.dim')) && attr(x, 'sp.dim') == '1D')
    plot.1d(values, xlim=xlim, ylim = ylim,
            x1vals = unique(round(locations[,1], dig=12)), ...)
  else{
    locations <- prepare.graph.krige.bayes(obj=x, locations=locations,
                                           borders=borders,
                                           borders.obj = eval(attr(x,"borders")),
                                           values.to.plot=values.to.plot,
                                           number.col = number.col,
                                           xlim = xlim, ylim = ylim)
    pty.prev <- par()$pty
    par(pty = "s")
    do.call("image", c(list(x=locations$x, y=locations$y, z=locations$values,
                          xlim = locations$coords.lims[,1], ylim = locations$coords.lims[,2]), ldots))
    if(!is.null(coords.data)) points(coords.data)
    if(!is.null(borders.arg)) polygon(borders, lwd=2)
    dots.l <- list(...)
    if(is.null(dots.l$col)) dots.l$col <- heat.colors(12)
    if(!is.null(x.leg) & !is.null(y.leg)){
      legend.krige(x.leg=x.leg, y.leg=y.leg,
                   values=locations$values,
                   vertical = vertical, cex=cex.leg,
                   col=dots.l$col, ...)
    }
  }
  par(pty=pty.prev)
  return(invisible())
}

"persp.glm.krige.bayes" <-
  function (x, locations, borders, 
            values.to.plot = c("median", "uncertainty",
              "quantiles", "probabilities", "simulation"), number.col, ...) 
{
  if(missing(x)) x <- NULL
  attach(x)
  on.exit(detach(x))
  if(missing(locations)) locations <-  eval(attr(x, "prediction.locations"))
  if(is.null(locations)) stop("prediction locations must be provided")
  if(ncol(locations) != 2) stop("locations must be a matrix or data-frame with two columns")
  if(!is.numeric(values.to.plot)){
    values.to.plot <- match.arg(values.to.plot,
                                choices = c("median", "uncertainty",
                                  "quantiles", "probabilities", "simulation"))
  }
  if(missing(borders)) borders <- NULL
  if(missing(number.col)) number.col <- NULL
  if(!is.null(attr(x, 'sp.dim')) && attr(x, 'sp.dim') == '1D')
    plot.1d(values, xlim=xlim, ylim = ylim,
            x1vals = unique(round(locations[,1], dig=12)), ...)
  else{
    locations <- prepare.graph.krige.bayes(obj=x, locations=locations,
                                         borders=borders,
                                         borders.obj = eval(attr(x,"borders")),
                                         values.to.plot=values.to.plot,
                                         number.col = number.col)
    persp(locations$x, locations$y, locations$values, ...)
  }
  return(invisible())
}


### The cases where some of the parameters are fixed are not implemented yet.

"hist.glm.krige.bayes" <-
  function(x, pars, density.est = TRUE,
           histogram = TRUE, ...)
{
  ## we need to construct the object such that it fits into geoR function hist.krige.bayes
  kb <- list(posterior=list(),prior=list())
  kb$prior$tausq.rel <- list(status="fixed") 
  kb$prior$phi <- list(status=ifelse(!is.null(x$prior$phi$status), x$prior$phi$status, "random"))
  kb$prior$sigmasq <- list(status=ifelse(!is.null(x$prior$sigmasq$status),x$prior$sigmasq$status, "random"))
  if(is.vector(x$posterior$beta$sample)){
    n.sim <- length(x$posterior$beta$sample)
    if(kb$prior$phi$status =="random"){
      kb$posterior$sample <- as.data.frame(cbind(x$posterior$beta$sample,
                                                 x$posterior$sigmasq$sample, x$posterior$phi$sample)) 
    }
    else{
      if(kb$prior$sigmasq$status =="random"){
        kb$posterior$sample <- as.data.frame(cbind(x$posterior$beta$sample,
                                                   x$posterior$sigmasq$sample, rep(9999,n.sim))) 
      }
      else kb$posterior$sample <- as.data.frame(cbind(x$posterior$beta$sample,
                                                      rep(-99,n.sim), rep(9999,n.sim))) 
    }
    names(kb$posterior$sample) <- c("beta", "sigmasq", "phi")
  }
  else{
    n.sim <- nrow(x$posterior$beta$sample)
    if(kb$prior$phi$status =="random"){
      kb$posterior$sample <- as.data.frame(cbind(t(x$posterior$beta$sample),
                                                 x$posterior$sigmasq$sample, x$posterior$phi$sample))
    }
    else{
      if(kb$prior$sigmasq$status =="random"){
        kb$posterior$sample <- as.data.frame(cbind(t(x$posterior$beta$sample),
                                                   x$posterior$sigmasq$sample, rep(9999,n.sim)))
      }
      else kb$posterior$sample <- as.data.frame(cbind(t(x$posterior$beta$sample),
                                                      rep(-99,n.sim), rep(9999,n.sim)))
    }
    names(kb$posterior$sample) <- c(names(x$posterior$beta$mean), "sigmasq", "phi")
  }
  kb$posterior$sample$tausq.rel <- rep(-99,n.sim)
  hist.krige.bayes(x=kb, pars=pars, density.est = density.est, histogram = histogram)
  return(invisible())
}


"mcmc.bayes.binom.logit" <- 
  function(data, units.m, trend, mcmc.input, messages.screen, cov.model, kappa, tausq.rel, coords, ss.sigma, df, phi.prior, phi.discrete)
{
  ##
  ## This is the MCMC engine for the Bayesian analysis of a spatial binomial logit Normal model
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  if(any(mcmc.input$S.start=="default")) {
    S <- as.vector(ifelse(data > 0, log(data), -1) - ifelse(units.m-data > 0, log(units.m-data), -1) )
  }
  else{
    if(!any(mcmc.input$S.start=="random")){
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        S <- as.vector(mcmc.input$S.start)
      }
      else  stop(" S.start must be a vector of same dimension as data ")
    }
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
  if(any(mcmc.input$phi.start=="default")) phi <- median(phi.discrete)
  else  phi <- mcmc.input$phi.start
  nmphi <-  length(phi.discrete)
  if(is.null(mcmc.input$phi.scale)) {
    if(nmphi > 1) stop("mcmc.input$phi.scale not given ")
    else phi.scale <- 0
  }
  else {
    phi.scale <- mcmc.input$phi.scale
    if(nmphi > 1 && pnorm((phi.discrete[nmphi] - phi.discrete[1])/(nmphi - 1), sd = sqrt(phi.scale)) > 0.975)
      warning("Consider making the grid in phi.discrete more dense. The algorithm may have problems moving.")
  }
  messages.C <- ifelse(messages.screen,1,0) ## for C-code
  ##
  ##                                                                      
  ## ---------- sampling ----------- ###### 
  cov.model.number <- cor.number(cov.model)
  if(is.vector(trend)) beta.size <- 1
  else beta.size <- ncol(trend)
  n.sim <- floor(n.iter/thin)
  ## remeber this rather odd coding for telling that S.start is from the prior !!!
  if(any(mcmc.input$S.start=="random")) Sdata <- as.double(as.vector(c(rep(0, n.sim*n - 1),1)))
  else Sdata <- as.double(as.vector(c(S, rep(0, (n.sim - 1) * n))))
  result <-  .C("mcmcrun4binom",
                as.integer(n),
                as.double(data),
                as.double(units.m),
                as.double(as.vector(t(trend))),
                as.integer(beta.size),
                as.integer(cov.model.number),
                as.double(kappa),
                as.double(tausq.rel),
                as.double(coords[,1]),
		as.double(coords[,2]),
                as.double(S.scale),
                as.double(phi.scale),
                as.integer(n.iter),
                as.integer(thin),
                as.integer(burn.in),
                as.integer(messages.C),
                as.double(ss.sigma),
                as.integer(df),
                as.double(phi.prior),
                as.double(phi.discrete),
                as.integer(nmphi),
                Sdata = Sdata,
                phi.sample = as.double(rep(phi, n.sim)),
		acc.rate = rep(0,floor(n.iter/1000)+1), 
		acc.rate.phi = rep(0,floor(n.iter/1000)+1), DUP=FALSE, PACKAGE = "geoRglm")[c("Sdata", "phi.sample","acc.rate","acc.rate.phi" )]
  attr(result$Sdata, "dim") <- c(n, n.sim)
  if(nmphi>1) result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate,result$acc.rate.phi))
  else result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate))
  result$acc.rate.phi <- NULL
  if(burn.in==0) result$acc.rate <- result$acc.rate[-1,]
  if(nmphi>1) names(result$acc.rate) <- c("iter.numb", "Acc.rate", "Acc.rate.phi")
  else names(result$acc.rate) <- c("iter.numb", "Acc.rate")
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n")) 
  return(result)
}

"mcmc.bayes.conj.binom.logit" <- 
  function(data, units.m, meanS, ttvbetatt, mcmc.input, messages.screen, cov.model, kappa, tausq.rel, coords, ss.sigma, df, phi.prior, phi.discrete)
{
  ##
  ## This is the MCMC engine for the Bayesian analysis (with normal prior for beta) of a spatial binomial logit Normal model
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  if(any(mcmc.input$S.start=="default")) {
    S <- as.vector(ifelse(data > 0, log(data), -1) - ifelse(units.m-data > 0, log(units.m-data), -1) ) - meanS
  }
  else{
    if(!any(mcmc.input$S.start=="random")){
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        S <- as.vector(mcmc.input$S.start)
      }
      else  stop(" S.start must be a vector of same dimension as data ")
    }
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
  if(any(mcmc.input$phi.start=="default")) phi <- median(phi.discrete)
  else  phi <- mcmc.input$phi.start
  nmphi <-  length(phi.discrete)
  if(is.null(mcmc.input$phi.scale)) {
    if(nmphi > 1) stop("mcmc.input$phi.scale not given ")
    else phi.scale <- 0
  }
  else {
    phi.scale <- mcmc.input$phi.scale
    if(nmphi > 1 && pnorm((phi.discrete[nmphi] - phi.discrete[1])/(nmphi - 1), sd = sqrt(phi.scale)) > 0.975)
      warning("Consider making the grid in phi.discrete more dense. The algorithm may have problems moving.")
  }
  messages.C <- ifelse(messages.screen,1,0) ## for C-code
  ##
  ## ---------- sampling ----------- ###### 
  cov.model.number <- cor.number(cov.model)
  if(is.null(ttvbetatt)) ttvbetatt <- matrix(0,beta.size,beta.size)
  n.sim <- floor(n.iter/thin)
  ## remeber this rather odd coding for telling that S.start is from the prior !!!
  if(any(mcmc.input$S.start=="random")) Sdata <- as.double(as.vector(c(rep(0, n.sim*n - 1),1)))
  else Sdata <- as.double(as.vector(c(S, rep(0, (n.sim - 1) * n))))
  result <-  .C("mcmcrun5binom",
                as.integer(n),
                as.double(data),
                as.double(units.m),
                as.double(as.vector(meanS)),
                as.double(as.vector(ttvbetatt)),
                as.integer(cov.model.number),
                as.double(kappa),
                as.double(tausq.rel),
                as.double(coords[,1]),
		as.double(coords[,2]),
                as.double(S.scale),
                as.double(phi.scale),
                as.integer(n.iter),
                as.integer(thin),
                as.integer(burn.in),
                as.integer(messages.C),
                as.double(ss.sigma),
                as.integer(df),
                as.double(phi.prior),
                as.double(phi.discrete),
                as.integer(nmphi),
                Sdata = Sdata,
                phi.sample= as.double(rep(phi, n.sim)),
		acc.rate = rep(0,floor(n.iter/1000)+1), 
		acc.rate.phi = rep(0,floor(n.iter/1000)+1), DUP=FALSE, PACKAGE = "geoRglm")[c("Sdata", "phi.sample","acc.rate","acc.rate.phi" )]
  attr(result$Sdata, "dim") <- c(n, n.sim)
  if(nmphi>1) result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate,result$acc.rate.phi))
  else result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate))
  result$acc.rate.phi <- NULL
  if(burn.in==0) result$acc.rate <- result$acc.rate[-1,]
  if(nmphi>1) names(result$acc.rate) <- c("iter.numb", "Acc.rate", "Acc.rate.phi")
  else names(result$acc.rate) <- c("iter.numb", "Acc.rate")
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n"))
  return(result)
}


"binom.krige.bayes" <- 
  function(geodata, coords = geodata$coords, data = geodata$data, units.m = "default", locations = "no", model, prior, mcmc.input, output){
###########
  if(missing(geodata))
    geodata <- list(coords=coords, data=data, units.m=units.m)
  call.fc <- match.call()
  seed <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
  do.prediction <- ifelse(all(locations == "no"), FALSE, TRUE)
  ##
  ## Checking data configuration
  ##
  if(is.vector(coords)) {
    coords <- cbind(coords, 0)
    warning("vector of coordinates: one spatial dimension assumed")
  }
  coords <- as.matrix(coords)
  n <- length(data)
  if(nrow(coords) != n) stop("number of data is different from number of data locations (coordinates)")
  if(any(units.m == "default")){
    if(!is.null(geodata$units.m)) units.m <- geodata$units.m
    else units.m <- rep(1, n)
  }
  ##
  ## reading input
  ##
  if(missing(model)) model <- model.glm.control()
  else model <- model.glm.check.aux(model, fct = "binom.krige.bayes")
  cov.model <- model$cov.model
  kappa <- model$kappa
  tausq.rel <- prior$tausq.rel
  ## reading prior input
  ##
  if(missing(prior))  stop("binom.krige.bayes: argument prior must be given ")
  else prior <- prior.glm.check.aux(prior, fct = "binom.krige.bayes")
  beta.prior <- prior$beta.prior
  beta <- prior$beta
  beta.var <- prior$beta.var.std
  sigmasq.prior <- prior$sigmasq.prior
  if(sigmasq.prior == "fixed") sigmasq <- prior$sigmasq
  else{
    df.sigmasq <- prior$df.sigmasq
    S2.prior <- prior$sigmasq
  }
  phi.prior <- prior$phi.prior 
  phi <- prior$phi
  if(phi.prior != "fixed") phi.discrete <- prior$phi.discrete
  else phi.discrete <- phi
  ##
  ## reading output options
  ##
  if(missing(output)) output <- output.glm.control()
  else output <- output.glm.check.aux(output, fct = "binom.krige.bayes")
  quantile.estimator <- output$quantile.estimator
  probability.estimator <- output$probability.estimator
  inference <- output$inference
  messages.screen <- output$messages.screen
  ## check == here
  data.dist <- as.vector(dist(coords))
  if(1000000000000. * min(data.dist) == 0) stop("Two coords are identical; not allowed.")
  ##
  trend.d <- model$trend.d
  if(messages.screen) {
    cat(switch(as.character(trend.d)[1],
                 "cte" = "binom.krige.bayes: model with mean being constant",
                 "1st" = "binom.krige.bayes: model with mean given by a 1st order polynomial on the coordinates",
                 "2nd" = "binom.krige.bayes: model with mean given by a 2nd order polynomial on the coordinates",
                 "binom.krige.bayes: model with mean defined by covariates provided by the user"))
    cat("\n")
  }
  trend.data <- unclass(trend.spatial(trend=trend.d, geodata = geodata))
  dimnames(coords) <- list(NULL, NULL)
  dimnames(trend.data) <- list(NULL, NULL)
  beta.size <- ncol(trend.data)
  if(nrow(trend.data) != n) stop("length of trend is different from the length of the data")
  if(beta.size > 1)
    beta.names <- paste("beta", (0:(beta.size-1)), sep="")
  else beta.names <- "beta"
  ##
  if(beta.prior == "normal" |  beta.prior == "fixed"){
    if(beta.size != length(beta))
      stop("binom.krige.bayes: size of beta incompatible with the trend model (covariates)")
  }
  if(sigmasq.prior == "uniform"){
    if(max(units.m) == 1) warning("This choice of sigmasq.prior gives an improper posterior !!!!!!! \n")
    if(sum(ifelse(units.m>1,1,0)) < beta.size + 3) warning("This choice of sigmasq.prior may give an improper posterior !!!!!!! \n")
  }
  aniso.pars <- model$aniso.par
  if(!is.null(aniso.pars)) coords.transf <- coords.aniso(coords = coords, aniso.pars = aniso.pars)
  else coords.transf <- coords
  ##
  # checking prediction locations
  ##
  if((inference) & (do.prediction)){
    ## Checking the consistency between coords, locations, and trends
    trend.l <- model$trend.l
    if(is.vector(locations)){
      if(length(locations) == 2) {
        locations <- t(as.matrix(locations))
        warning("only one location to be predicted (in two-dimensional space) \n")
      }
      else locations <- as.matrix(cbind(locations, 0))
    }
    else locations <- as.matrix(locations)
    ni <- nrow(locations)
    ## Checking for 1D prediction 
    if(length(unique(locations[,1])) == 1 | length(unique(locations[,2])) == 1)
      krige1d <- TRUE
    else krige1d <- FALSE
    ##
    if(is.null(trend.l)) stop("trend.l needed for prediction")
    if(inherits(trend.d, "formula") | inherits(trend.l, "formula")){
      if((!inherits(trend.d, "formula")) | (!inherits(trend.l, "formula")))
        stop("trend.d and trend.l must have similar specification\n")
    }
    else{
      if((class(trend.d)=="trend.spatial") & (class(trend.l)=="trend.spatial")){
        if(ncol(trend.d) != ncol(trend.l))
          stop("trend.d and trend.l do not have the same number of columns")
      }
      else if(trend.d != trend.l) stop("trend.l is different from trend.d")
    }
    if(nrow(unclass(trend.spatial(trend=model$trend.l, geodata = list(coords = locations)))) != ni)
      stop("binom.krige.bayes: number of points to be estimated is different of the number of trend locations")
    kb.results <- list(posterior = list(), predictive = list())
  }
  else {
    if(do.prediction & messages.screen) cat(paste("need to specify inference=TRUE to make predictions \n"))
    kb.results <- list(posterior = list(), predictive = paste("prediction not performed"))
    do.prediction <- FALSE
  }
  ##
  ## ##### preparing for MCMC -------------------------------------------------------
  ##
  if(missing(mcmc.input)) stop("binom.krige.bayes: argument mcmc.input must be given")
  mcmc.input <- mcmc.check.aux(mcmc.input, fct="binom.krige.bayes")
  ##
  if(beta.prior == "fixed" | beta.prior == "normal") mean.d <- as.vector(trend.data %*% beta)
  else mean.d <- rep(0, n)
  if(sigmasq.prior != "fixed"){
    if(beta.prior == "flat") df.model <- n - beta.size + df.sigmasq
    else df.model <- n + df.sigmasq
  }
  else df.model <- Inf
  if(beta.prior == "normal"){
    if(beta.size > 1) ttvbetatt <- trend.data%*%beta.var%*%t(trend.data)
    else ttvbetatt <- crossprod(t(trend.data))*beta.var
  }  
  else ttvbetatt <- NULL
  if(sigmasq.prior == "fixed"){     ### implies that phi is fixed !
    invcov <- varcov.spatial(coords = coords, cov.model = cov.model, kappa = kappa, nugget = tausq.rel*sigmasq,
                               cov.pars = c(sigmasq,phi), inv = TRUE, func.inv = "cholesky",
                               try.another.decomposition = FALSE)$inverse
    if(beta.prior != "fixed"){
      ivtt <- invcov%*%trend.data
      if(beta.prior == "normal") invcov <- invcov-ivtt%*%solve.geoR(crossprod(trend.data, ivtt) + solve(beta.var), t(ivtt))
      else invcov <- invcov-ivtt%*%solve.geoR(crossprod(trend.data, ivtt), t(ivtt))
    }
  }
  if((phi.prior == "fixed") && (sigmasq.prior != "fixed")){
    phi.prior.prob <- 1
    phi.discrete <- phi
  }
  else phi.prior.prob <-  prior$priors.info$phi$probs
  ##
############----------PART 2 ------------##############################
############-----------MCMC -------------##############################
  ##
  if(sigmasq.prior == "fixed"){ 
    log.odds <- mcmc.binom.logit(data = data, units.m = units.m, meanS = mean.d, invcov=invcov, mcmc.input = mcmc.input, messages.screen=messages.screen)
  }
  else {
    kb.results$posterior$phi <- list()
    ## take care re-using log.odds !
    if(beta.prior == "flat"){
      log.odds <- mcmc.bayes.binom.logit(data=data, units.m=units.m, trend=trend.data, mcmc.input=mcmc.input, messages.screen=messages.screen, cov.model=cov.model, 
                                         kappa=kappa, tausq.rel = tausq.rel, coords=coords.transf, 
                                         ss.sigma = df.sigmasq*S2.prior, df = df.model, phi.prior = phi.prior.prob,
                                         phi.discrete = phi.discrete)
    }
    else{     
      log.odds <- mcmc.bayes.conj.binom.logit(data=data, units.m=units.m, meanS = mean.d, ttvbetatt = ttvbetatt, mcmc.input=mcmc.input, messages.screen=messages.screen,
                                              cov.model=cov.model, kappa=kappa, tausq.rel = tausq.rel,
                                              coords=coords.transf, ss.sigma = df.sigmasq*S2.prior, df = df.model,
                                              phi.prior = phi.prior.prob, phi.discrete = phi.discrete)
    }
    kb.results$posterior$phi$sample <- log.odds$phi.sample
  }
  kb.results$posterior$acc.rate  <- log.odds$acc.rate
  log.odds <- log.odds$Sdata
  ##
##############-------------PART 3----------######################
##############------------prediction-------######################
  ##
  n.sim <- ncol(log.odds)
  if(inference) {
    if(phi.prior=="fixed") phi.posterior <- list(phi.prior=phi.prior, phi=phi)
    else  phi.posterior <- list(phi.prior=phi.prior, phi.discrete=phi.discrete, sample=kb.results$posterior$phi$sample)
    predict.temp <- pred.aux(S=log.odds, coords=coords, locations=locations, model=model, prior=prior, output=output, phi.posterior=phi.posterior, link="logit")
    temp.post <- predict.temp$temp.post
    if(do.prediction) {
      temp.pred <- predict.temp$temp.pred
      kb.results$predictive$simulations <- predict.temp$pred.simulations
    }
    ##
    if(do.prediction) {
      ##
      d0mat <- loccoords(coords, locations)
      loc.coincide <- (colSums(d0mat < 1e-10) == 1)

      if((is.logical(quantile.estimator) && (quantile.estimator)) || (is.numeric(quantile.estimator))){
        predi.q <- pred.quan.aux(temp.pred, loc.coincide, df.model, ni, quantile.estimator)
        kb.results$predictive$median <- plogis(predi.q$median)
        kb.results$predictive$uncertainty <- (plogis(predi.q$upper) - plogis(predi.q$lower))/4      
        if(is.data.frame(predi.q$quantiles)){
          names.q <- names(predi.q$quantiles)
          kb.results$predictive$quantiles <- as.data.frame(plogis(as.matrix(predi.q$quantiles)))
          names(kb.results$predictive$quantiles) <- names.q
        }
        else kb.results$predictive$quantiles <- plogis(predi.q$quantiles)
      }
      ##
      ## ------ probability estimators
      ##
      if(!is.null(probability.estimator)) {
        logit.probab <- ifelse(probability.estimator < 1, log(probability.estimator) - log(1-probability.estimator), 1e+17)
        logit.probab <- ifelse(probability.estimator > 0, logit.probab, 1e-17)
        len.p <- length(probability.estimator)
        if(len.p== 1){
          kb.results$predictive$probability <- round(pmixed(logit.probab, temp.pred, df.model), digits = 3)
        }
        else{
          kb.results$predictive$probability <- matrix(NA,ni,len.p)
          for(ii in seq(length=len.p)){
            kb.results$predictive$probability[,ii] <- round(pmixed(logit.probab[ii], temp.pred, df.model), digits = 3)
          }
        }
      }
      remove("temp.pred")
      ## 
      if(messages.screen) cat("binom.krige.bayes: Prediction performed \n")
    }
    else {
      kb.results$predictive <- "no locations to perform prediction were provided"
      if(messages.screen) cat(paste("Only Bayesian estimation of model parameters "))
    }
    ##
    ##----- calculating posterior summaries ----------------##
    ##
    if(beta.prior == "fixed") kb.results$posterior$beta <- paste("provided by user: ", beta)
    else {
      kb.results$posterior$beta <- list()
      kb.results$posterior$beta$mean <- rowMeans(temp.post$beta.mean)
      names(kb.results$posterior$beta$mean) <- beta.names
      kb.results$posterior$beta$var <- rowMeans(temp.post$beta.var, dims = 2) + var(t(temp.post$beta.mean))
      dimnames(kb.results$posterior$beta$var) <- list(beta.names,beta.names)
    }
    if(sigmasq.prior == "fixed") kb.results$posterior$sigmasq <- paste("provided by user: ", sigmasq) 
    else{
      kb.results$posterior$sigmasq <- list()
      kb.results$posterior$sigmasq$mean <- mean(temp.post$S2)*df.model/(df.model-2)
      kb.results$posterior$sigmasq$var <- (mean(temp.post$S2)*2/(df.model-4) + var(temp.post$S2))*df.model^2/(df.model-2)^2
    }
    if(phi.prior == "fixed") kb.results$posterior$phi <- paste("provided by user: ", phi) 
    else{
      kb.results$posterior$phi$mean <- mean(kb.results$posterior$phi$sample)
      kb.results$posterior$phi$var <- var(kb.results$posterior$phi$sample)
    }
    ##
    ## Simulations from the posterior of parameters.
    ##
    if(output$sim.posterior){
      if(beta.size == 1) {
        if(sigmasq.prior == "fixed") {
          if(beta.prior != "fixed")
            kb.results$posterior$beta$sample <- rnorm(n.sim) * as.vector(sqrt(temp.post$beta.var)) + as.vector(temp.post$beta.mean)
        }
        else{
          kb.results$posterior$sigmasq$sample <- rinvchisq(n.sim, df.model, temp.post$S2)
          if(beta.prior != "fixed"){
            cond.beta.sd <- sqrt((as.vector(temp.post$beta.var) * kb.results$posterior$sigmasq$sample)/temp.post$S2)
            kb.results$posterior$beta$sample <- rnorm(n.sim) * cond.beta.sd + as.vector(temp.post$beta.mean)
          }
        }
      }
      else {
        if(sigmasq.prior == "fixed"){
          if(beta.prior != "fixed")
            kb.results$posterior$beta$sample <- array(apply(temp.post$beta.var,3,multgauss),dim=c(beta.size, n.sim))+temp.post$beta.mean
        }
        else {
          kb.results$posterior$sigmasq$sample <- rinvchisq(n.sim, df.model, temp.post$S2)
          if(beta.prior != "fixed"){
            if(is.R()) cond.beta.var <- temp.post$beta.var *rep(kb.results$posterior$sigmasq$sample/temp.post$S2,rep(beta.size^2,n.sim))
            else cond.beta.var <- temp.post$beta.var *rep(kb.results$posterior$sigmasq$sample/temp.post$S2,each = beta.size^2)
            kb.results$posterior$beta$sample <- array(apply(cond.beta.var,3,multgauss),dim=c(beta.size, n.sim)) + temp.post$beta.mean
          }
        }
      }
    }
    remove("temp.post")
  }
  if(output$keep.mcmc.sim) kb.results$posterior$simulations <- plogis(log.odds)  
  model$lambda <- NULL
  kb.results$model <- model
  kb.results$prior <- prior$priors.info
  kb.results$mcmc.input <- mcmc.input
  kb.results$.Random.seed <- seed
  kb.results$call <- call.fc
  attr(kb.results, "prediction.locations") <- call.fc$locations
  if(do.prediction) attr(kb.results, 'sp.dim') <- ifelse(krige1d, "1d", "2d")
  if(!is.null(call.fc$borders)) attr(kb.results, "borders") <- call.fc$borders
  class(kb.results) <- "glm.krige.bayes"
  return(kb.results)
}



"mcmc.binom.aux" <- function(z, data, units.m, meanS, QQ, S.scale, nsim, thin, QtivQ)
{
  ##
###### ------------------------ doing the mcmc-steps -----------
###### 
  ##
  n <- length(data)
  result <- .C("mcmc1binom",
               as.integer(n),
               z = as.double(z),
               S = as.double(rep(0, nsim * n)),
               as.double(data),
	       as.double(units.m),
               as.double(meanS),
               as.double(as.vector(t(QQ))),
               as.double(as.vector(QtivQ)),
               as.double(rnorm(n * nsim * thin) * sqrt(S.scale)),
               as.double(runif(nsim * thin)),
               as.double(S.scale),
               as.integer(nsim),
               as.integer(thin),
               acc.rate = as.double(1), DUP=FALSE, PACKAGE = "geoRglm")[c("z", "S", "acc.rate")]
  attr(result$S, "dim") <- c(n, nsim) 
  return(result)
}

"mcmc.binom.logit" <- function(data, units.m, meanS, invcov, mcmc.input, messages.screen)
{
####
  ## This is the MCMC engine for the spatial Binomial logit Normal model ----
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  fisher.l <- data*(1-data/units.m)
  ## gives a singular matrix when : beta.prior="flat", data are binary and all units.m=1. Therefore a fix.
  if(all(round(fisher.l)==0)) fisher.l <- data*0.025
  QQ <- t(chol(solve(invcov + diag(fisher.l))))
  sqrtfiQ <- sqrt(fisher.l)*QQ 
  QtivQ <- diag(n)-crossprod(sqrtfiQ)  
  if(any(mcmc.input$S.start=="default")){
    signum <- round(2*data/units.m) -1
    S <- as.vector(ifelse(data > 0 & data < units.m, qlogis(data/units.m), signum*1.96) - meanS)
    z <- as.vector(solve(QQ,S))
  }
  else{
    if(any(mcmc.input$S.start=="random")) z <- rnorm(n)
    else{
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        else z <- as.vector(solve(QQ,mcmc.input$S.start))
      }
      else stop(" S.start must be a vector of same dimension as data ")
    }
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
## ---------------- burn-in ----------------- ######### 
  if(burn.in > 0) {
    mcmc.output <- mcmc.binom.aux(z, data, units.m, meanS, QQ, S.scale, 1, burn.in, QtivQ)
    if(messages.screen) cat(paste("burn-in = ", burn.in, " is finished. Acc.-rate = ", round(mcmc.output$acc.rate, digits=3), "\n"))
    acc.rate.burn.in <- c(burn.in, mcmc.output$acc.rate)
  }
  else mcmc.output <- list(z = z)
##### ---------- sampling periode ----------- ###### 
  if(n.iter <= 1000) {
    n.temp <- round(n.iter/thin)
    n.turn <- 1
  }
  else {
    n.temp <- round(1000/thin)
    n.turn <- round(n.iter/1000)
  }
  n.sim <- n.turn * n.temp
  Sdata <- matrix(NA, n, n.sim)
  acc.rate <- matrix(NA, n.turn, 2)
  for(i in seq(length=n.turn)) {
    mcmc.output <- mcmc.binom.aux(mcmc.output$z, data, units.m, meanS, QQ, S.scale, n.temp, thin, QtivQ)
    Sdata[, seq((n.temp * (i - 1) + 1),(n.temp * i))] <- mcmc.output$S+meanS    
    if(messages.screen) cat(paste("iter. numb.", i * n.temp * thin+burn.in, " : Acc.-rate = ", round(mcmc.output$acc.rate, digits=3), "\n"))
    acc.rate[i,1] <-  i * n.temp * thin
    acc.rate[i,2] <- mcmc.output$acc.rate
  }
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n"))
  if(burn.in > 0) acc.rate <- as.data.frame(rbind(acc.rate.burn.in,acc.rate))
  else acc.rate <- as.data.frame(acc.rate)
  names(acc.rate) <- c("iter.numb", " Acc.rate")
#########
  return(list(Sdata=Sdata, acc.rate=acc.rate))
}


"binom.krige" <- function(geodata, coords = geodata$coords, data = geodata$data, units.m = "default", locations = NULL, borders = NULL, mcmc.input, krige, output)
{
  if(missing(geodata))
    geodata <- list(coords=coords, data=data, units.m=units.m)
  call.fc <- match.call()
  n <- length(data)
  if(any(units.m == "default")){
    if(!is.null(geodata$units.m)) units.m <- geodata$units.m
    else units.m <- rep(1, n)
  }
  if(missing(krige)) stop("must provide object krige")
  krige <- krige.glm.check.aux(krige,fct="binom.krige")
  cov.model <- krige$cov.model
  kappa <- krige$kappa
  beta <- krige$beta
  cov.pars <- krige$cov.pars
  nugget <- krige$nugget
  micro.scale <- krige$micro.scale
  aniso.pars <- krige$aniso.pars
  trend.d <- krige$trend.d
  trend.l <- krige$trend.l
  dist.epsilon <- krige$dist.epsilon
  if(krige$type.krige == "ok") beta.prior <- "flat"
  if(krige$type.krige == "sk") beta.prior <- "deg"
  if(missing(output)) output <- output.glm.control()
  else output <- output.glm.check.aux(output, fct = "binom.krige")
  sim.predict <- output$sim.predict
  messages.screen <- output$messages.screen
  ##
  if(is.vector(coords)){
    coords <- cbind(coords, 0)
    warning("vector of coordinates: one spatial dimension assumed")
  }
  coords <- as.matrix(coords)
  dimnames(coords) <- list(NULL, NULL)
  ## Checking for 1D prediction 
  if(length(unique(locations[,1])) == 1 | length(unique(locations[,2])) == 1)
    krige1d <- TRUE
  else krige1d <- FALSE
  ##
  if(is.null(locations)) {
    if(messages.screen) cat(paste("locations need to be specified for prediction; prediction not performed \n"))
  }
  else {
    if(is.null(trend.l))
      stop("trend.l needed for prediction")
  }
  trend.data <- unclass(trend.spatial(trend=trend.d, geodata = geodata))
  beta.size <- ncol(trend.data)
  if(nrow(trend.data) != n) stop("length of trend is different from the length of the data")
  if(beta.prior == "deg")
    if(beta.size != length(beta))
      stop("size of mean vector is incompatible with trend specified") 
  if(beta.size > 1)
    beta.names <- paste("beta", (0:(beta.size-1)), sep="")
  else beta.names <- "beta"
  ##
  ## preparing for MCMC 
  ##
  if(missing(mcmc.input)) stop("binom.krige: argument mcmc.input must be given")
  mcmc.input <- mcmc.check.aux(mcmc.input, fct="binom.krige")
  ##
  if(beta.prior == "deg") mean.d <-  as.vector(trend.data %*% beta)
  else mean.d <- rep(0,n)
  if(!is.null(aniso.pars)) {
    invcov <- varcov.spatial(coords = coords.aniso(coords = coords, aniso.pars = aniso.pars), cov.model = cov.model, kappa = kappa, 
                             nugget = nugget, cov.pars = cov.pars, inv = TRUE, func.inv = "cholesky",
                             try.another.decomposition = FALSE)$inverse
  }
  else {
    invcov <- varcov.spatial(coords = coords, cov.model = cov.model, kappa = kappa, nugget = nugget, cov.pars = cov.pars,
                             inv = TRUE, func.inv = "cholesky", try.another.decomposition = FALSE)$inverse
  }
  ##
########################----- MCMC ------#####################
  ##
  if(beta.prior == "flat") {
    ivtt <- invcov%*%trend.data
    invcov <- invcov-ivtt%*%solve.geoR(crossprod(trend.data, ivtt),t(ivtt))
  }
  res.mcmc <- mcmc.binom.logit(data = data, units.m = units.m, meanS= mean.d, invcov=invcov, mcmc.input = mcmc.input, messages.screen=messages.screen)
  acc.rate <- res.mcmc$acc.rate
  ##
  ##------------------------------------------------------------
######################## ---- prediction ----- #####################
  if(!is.null(locations)) {
    if(!is.null(borders)){
      locations <- locations.inside(locations, borders)
      if(nrow(locations) == 0) stop(" binom.krige : there are no prediction locations inside the borders")
      if(messages.screen) cat(" binom.krige: results will be returned only for prediction locations inside the borders\n")
    }
    krige <- list(type.krige = krige$type.krige, beta = beta, trend.d = trend.d, trend.l = trend.l, cov.model = cov.model, 
                  cov.pars = cov.pars, kappa = kappa, nugget = nugget, micro.scale = micro.scale, dist.epsilon = dist.epsilon, 
                  aniso.pars = aniso.pars, link = "logit")
    kpl.result <- glm.krige.aux(data = res.mcmc$Sdata, coords = coords, locations = locations, krige = krige,
					output = list(n.predictive = ifelse(sim.predict,1,0),
					      signal = TRUE, messages=FALSE))			   
    remove(list = c("res.mcmc"))
    kpl.result$krige.var <- rowMeans(kpl.result$krige.var) + apply(kpl.result$predict, 1, var)
    if(nrow(locations) > 1) kpl.result$mcmc.error <- sqrt(asympvar(kpl.result$predict)/ncol(kpl.result$predict))
    else kpl.result$mcmc.error <- sqrt(asympvar(as.vector(kpl.result$predict), messages = FALSE)/length(as.vector(kpl.result$predict)))
    kpl.result$predict <- rowMeans(kpl.result$predict)
    if(beta.prior == "flat") {
      kpl.result$beta.est <- rowMeans(kpl.result$beta)
      names(kpl.result$beta.est) <- beta.names
    }
    kpl.result$beta <- NULL
  }
  else{
    if(beta.prior == "flat") {
      ## GLS
      beta.est <- solve.geoR(crossprod(trend.data, ivtt),t(ivtt))%*%rowMeans(res.mcmc$Sdata)
      kpl.result <- list(prevalence=plogis(res.mcmc$Sdata), beta.est = beta.est, acc.rate=acc.rate)
    }
    else kpl.result <- list(prevalence=plogis(res.mcmc$Sdata), acc.rate=acc.rate)
  }
  kpl.result$call <- call.fc
#######################################
  attr(kpl.result, "prediction.locations") <- call.fc$locations
  if(!is.null(locations)) attr(kpl.result, 'sp.dim') <- ifelse(krige1d, "1d", "2d")
  if(!is.null(call.fc$borders)) attr(kpl.result, "borders") <- call.fc$borders
  class(kpl.result) <- "kriging"
  return(kpl.result)
}

"glm.krige.aux" <- 
function(data, coords, locations, krige, output)
{
  krige$lambda <- 1
  krige$link <- NULL 
  kc.result <- krige.conv.extnd(data = data, coords = coords, locations = locations, krige = krige, output = output)	
  ##
  ##################### Back-transforming predictions
  ## using second order taylor-expansion + facts for N(0,1) [third moment = 0 ; fourth moment = 12].
  ivlogit2 <- ifelse(kc.result$predict<700, exp(kc.result$predict)*(-expm1(kc.result$predict))/(1+exp(kc.result$predict))^3, 0)
  kc.result$predict <- plogis(kc.result$predict) + 0.5*ivlogit2*kc.result$krige.var
  kc.result$krige.var <- ifelse(kc.result$predict<700, exp(kc.result$predict)/(1+exp(kc.result$predict))^2, 0)^2*kc.result$krige.var+(11/4)*ivlogit2^2*kc.result$krige.var^2
  remove(list = c("ivlogit2"))
  if(output$n.predictive > 0) {
    kc.result$simulations <- plogis(kc.result$simulations)
  }
  return(kc.result)
}
"krige.bayes.extnd" <- 
  function(geodata, coords=geodata$coords, data=geodata$data,
           locations = "no", model, prior, output)
{
  ##
  ## ======================= PART 1 ==============================
  ##                Reading and Checking Input
  ## =============================================================
  ##
  ## setting output object and environments
  ##
  if(missing(geodata))
    geodata <- list(coords=coords, data=data)
  call.fc <- match.call()
  seed <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
  do.prediction <- ifelse(all(locations == "no"), FALSE, TRUE)
  base.env <- sys.frame(sys.nframe())
  message.prediction <- character()
  ##
  ## reading model input
  ##
  if(missing(model))
    model <- model.control()
  else{
    if(class(model) != "model.geoR"){
      if(!is.list(model))
        stop("krige.bayes.extnd: the argument model only takes a list or an output of the function model.control")
      else{
        model.names <- c("trend.d", "trend.l", "cov.model", "kappa", "aniso.pars", "lambda") 
        model <- object.match.names(model,model.names)
        if(is.null(model$trend.d)) model$trend.d <- "cte"  
        if(is.null(model$trend.l)) model$trend.l <- "cte"  
        if(is.null(model$cov.model)) model$cov.model <- "matern"  
        if(is.null(model$kappa)) model$kappa <- 0.5
        if(is.null(model$aniso.pars)) model$aniso.pars <- NULL 
        if(is.null(model$lambda)) model$lambda <- 1
        model <- model.control(trend.d = model$trend.d,
                               trend.l = model$trend.l,
                               cov.model = model$cov.model,
                               kappa = model$kappa,
                               aniso.pars = model$aniso.pars,
                               lambda = model$lambda)
      }
    }
  }
  cov.model <- model$cov.model
  cov.model.number <- cor.number(cov.model)
  kappa <- model$kappa
  ##
  ## reading prior input
  ##
  if(missing(prior))
    prior <- prior.control()
  else{
    if(class(prior) != "prior.geoR"){
      if(!is.list(prior))
        stop("krige.bayes.extnd: the argument prior only takes a list or an output of the function prior.control")
      else{
        prior.names <- c("beta.prior", "beta", "beta.var.std", "sigmasq.prior",
                         "sigmasq", "df.sigmasq", "phi.prior", "phi", "phi.discrete",
                         "tausq.rel.prior", "tausq.rel", "tausq.rel.discrete") 
        prior <- object.match.names(prior,prior.names)
        ## DO NOT CHANGE ORDER OF THE NEXT 3 LINES
        if(is.null(prior$beta)) prior$beta <-  NULL
        if(is.null(prior$beta.prior)) prior$beta.prior <-  c("flat", "normal", "fixed")
        if(is.null(prior$beta.var.std)) prior$beta.var.std <-  NULL
        ## DO NOT CHANGE ORDER OF THE NEXT 3 LINES
        if(is.null(prior$sigmasq)) prior$sigmasq <- NULL
        if(is.null(prior$sigmasq.prior))
          prior$sigmasq.prior <- c("reciprocal",  "uniform", "sc.inv.chisq",  "fixed") 
        if(is.null(prior$df.sigmasq)) prior$df.sigmasq <- NULL
        ## DO NOT CHANGE ORDER OF THE NEXT 3 LINES
        if(is.null(prior$phi)) prior$phi <- NULL
        if(is.null(prior$phi.prior))
          prior$phi.prior <- c("uniform", "exponential", "fixed", "squared.reciprocal","reciprocal")
        if(is.null(prior$phi.discrete)) prior$phi.discrete <- NULL
        ## DO NOT CHANGE ORDER OF THE NEXT 3 LINES
        if(is.null(prior$tausq.rel)) prior$tausq.rel <- 0
        if(is.null(prior$tausq.rel.prior))
          prior$tausq.rel.prior <- c("fixed", "uniform")
        if(is.null(prior$tausq.rel.discrete)) prior$tausq.rel.discrete <- NULL 
        prior <- prior.control(beta.prior = prior$beta.prior,
                               beta = prior$beta, beta.var.std = prior$beta.var.std,
                               sigmasq.prior = prior$sigmasq.prior,
                               sigmasq = prior$sigmasq,  df.sigmasq = prior$df.sigmasq,
                               phi.prior = prior$phi.prior,
                               phi = prior$phi, phi.discrete = prior$phi.discrete, 
                               tausq.rel.prior = prior$tausq.rel.prior,
                               tausq.rel = prior$tausq.rel,
                               tausq.rel.discrete = prior$tausq.rel.discrete)
      } 
    }
  }
  ##
  kb <- list(posterior = list(beta=list(), sigmasq=list(), phi=list(), tausq.rel=list()),
             predictive=list(mean = NULL, variance = NULL, distribution = NULL),
             prior = prior$priors.info, model = model) 
  ##class(kb$posterior) <- "krige.bayes.posterior"
  ##class(kb$predictive) <- "krige.bayes.predictive"
  ##class(kb$prior) <- "krige.bayes.prior"
  pred.env <- new.env()
  ##
  beta <- prior$beta
  if(prior$beta.prior == "fixed") beta.fixed <- beta
  if(prior$beta.prior == "normal"){
    beta.var.std <- prior$beta.var.std
    inv.beta.var.std <- solve.geoR(beta.var.std)
    betares <- list(iv = inv.beta.var.std, ivm = drop(solve.geoR(beta.var.std, beta)),
                    mivm = drop(crossprod(beta, solve.geoR(beta.var.std, beta))))
  }
  if(prior$sigmasq.prior == "fixed") sigmasq.fixed <- prior$sigmasq
  else S2.prior <- prior$sigmasq
  df.sigmasq <- prior$df.sigmasq
  ##
  if(prior$phi.prior != "fixed") stop("krige.bayes.extnd: only phi fixed is allowed.")
  ##
  tausq.rel.fixed <- tausq.rel <- prior$tausq.rel
  if(prior$tausq.rel.prior != "fixed") stop("krige.bayes.extnd: only tausq fixed is allowed.")
  ##
  ## checking data configuration
  ##
  data <- as.matrix(data)
  n.datasets <- ncol(data)
  n <- nrow(data)
  if(n.datasets == 1) stop("krige.bayes.extnd: this function is for multiple datasets. Use krige.bayes instead.")
  ##
  if(is.vector(coords)){
    coords <- cbind(coords, 0)
    warning("krige.bayes.extnd: vector of coordinates: assuming one spatial dimension (transect)")
  }
  coords <- as.matrix(coords)
  dists.env <- new.env()
  assign("data.dist", as.vector(dist(coords)), envir=dists.env)
  data.dist.range <- range(get("data.dist", envir=dists.env))
  data.dist.min <- data.dist.range[1]
  data.dist.max <- data.dist.range[2]
  if(1e12*data.dist.min < 0.5) stop("krige.bayes.extnd: this function does not allow two data at same location")
  ##
  ## reading output options
  ##
  if(missing(output))
    output <- output.control()
  else{
    if(class(output) != "output.geoR"){
      if(!is.list(output))
        stop("krige.bayes.extnd: the argument output only takes a list or an output of the function output.control")
      else{
        output.names <- c("n.posterior","n.predictive","moments","n.back.moments","simulations.predictive",
                          "mean.var","quantile","threshold","signal","messages.screen")
        output <- object.match.names(output,output.names)
        if(is.null(output$n.posterior)) output$n.posterior <- 1000 
        if(is.null(output$n.predictive)) output$n.predictive <- NULL
        if(is.null(output$moments)) output$moments <- TRUE
        if(is.null(output$n.back.moments)) output$n.back.moments <- 1000 
        if(is.null(output$simulations.predictive)){
          if(is.null(output$n.predictive)) output$simulations.predictive <- NULL
          else
            output$simulations.predictive <- ifelse(output$n.predictive > 0, TRUE, FALSE)
        }
        if(is.null(output$mean.var)) output$mean.var <- NULL
        if(is.null(output$quantile)) output$quantile <- NULL
        if(is.null(output$threshold)) output$threshold <- NULL
        if(is.null(output$signal)) output$signal <- NULL
        if(is.null(output$messages.screen)) output$messages.screen <- TRUE
        output <- output.control(n.posterior = output$n.posterior,
                                 n.predictive = output$n.predictive,
                                 moments = output$moments,
                                 n.back.moments = output$n.back.moments, 
                                 simulations.predictive = output$simulations.predictive,
                                 mean.var = output$mean.var, quantile = output$quantile,
                                 threshold = output$threshold, signal = output$signal,
                                 messages = output$messages.screen)
      }
    }
  }
  n.posterior <- output$n.posterior
  messages.screen <- output$messages.screen
  if(do.prediction){
    n.predictive <- as.integer(output$n.predictive)
    if(is.null(n.predictive)) n.predictive <- as.integer(0)
    simulations.predictive <- output$simulations.predictive
    if(is.null(simulations.predictive)) simulations.predictive <- FALSE
    if(!is.null(output$signal) && output$signal) stop("krige.bayes.extnd: prediction of the signal is not implemented")
    if(!is.null(output$probability.estimator)) stop("krige.bayes.extnd: probability.estimator not implemented\n")
    if(!is.null(output$quantile.estimator)) stop("krige.bayes.extnd: quantile.estimator not implemented\n")
    if(simulations.predictive && n.predictive == 0) n.predictive <- as.integer(1)
  }
  ##
  ## Box-Cox transformation
  ##
  if(abs(model$lambda-1)>0.0001) stop("krige.bayes.extnd: Box-Cox transformation not allowed \n")    
  ##
  ## Building trend (covariates/design) matrices:   
  ##
  dimnames(coords) <- list(NULL, NULL)
  if(nrow(coords) != n) stop("krige.bayes.extnd: number of data is different of number of data locations (coordinates)")
  trend.data <- unclass(trend.spatial(trend=model$trend.d, geodata = geodata))
  beta.size <- ncol(trend.data)
  if(nrow(trend.data) != n) stop("length of trend is different from the length of the data")
  if((prior$beta.prior == "normal") && (beta.size != length(beta)) )
    stop("krige.bayes.extnd: size of beta incompatible with the trend model (covariates)")
  ##
  if(do.prediction) {
    ##
    ## Checking the spatial dimension for prediction
    ##  1 (data/prediction on a transect) or 2 (data/prediction on an area)
    ##
    if(is.vector(locations)) {
      if(length(locations) == 2) {
        locations <- t(as.matrix(locations))
        warning("krige.bayes.extnd: FUNCTION IS CONSIDERING YOU HAVE ENTERED WITH 1 LOCATION TO BE PREDICTED IN A 2-DIM. REGION\n")
      }
      else locations <- as.matrix(cbind(locations, 0))
    }
    else locations <- as.matrix(locations)
    ##
    ## Checking trend specification
    ##
    if(inherits(model$trend.d, "formula") | inherits(model$trend.l, "formula")){
      if((inherits(model$trend.d, "formula") == FALSE) | (inherits(model$trend.l, "formula") == FALSE))
        stop("krige.bayes.extnd: model$trend.d and model$trend.l must have similar specification\n")
    }
    else{
      if((class(model$trend.d) == "trend.spatial") & (class(model$trend.l) == "trend.spatial")){
        if(ncol(model$trend.d) != ncol(model$trend.l))
          stop("krige.bayes.extnd: trend.d and trend.l do not have the same number of columns")
      }
      else
        if(model$trend.d != model$trend.l)
          stop("krige.bayes.extnd: especification of model$trend.l and model$trend.d must be similar")
    }
    ##
    if(messages.screen){
      cat(switch(model$trend.d,
                 "cte" = "krige.bayes.extnd: model with mean being constant",
                 "1st" = "krige.bayes.extnd: model with mean given by a 1st order polynomial on the coordinates",
                 "2nd" = "krige.bayes.extnd: model with mean given by a 2nd order polynomial on the coordinates",
                 "krige.bayes.extnd: model with mean defined by covariates provided by the user"))
      cat("\n")
    }
    ##
    dimnames(locations) <- list(NULL, NULL)
    assign("trend.loc", unclass(trend.spatial(trend=model$trend.l, geodata = list(coords = locations))), envir=pred.env)
    ni <- nrow(get("trend.loc", envir=pred.env))
    if(nrow(locations) != ni)
      stop("krige.bayes.extnd: number of points to be estimated is different of the number of trend locations")
  }
  ##
  ## Anisotropy correction
  ##   (warning: this must be placed here, AFTER trend matrices be defined)
  ##
  if(!is.null(model$aniso.pars)) {
    if((abs(model$aniso.pars[1]) > 0.001) & (abs(model$aniso.pars[2] - 1) > 0.001)){
      if(messages.screen) cat("krige.bayes.extnd: anisotropy parameters provided and assumed to be constants\n")
      coords <- coords.aniso(coords = coords, aniso.pars = model$aniso.pars)
      if(do.prediction) locations <- coords.aniso(coords = locations, aniso.pars = model$aniso.pars)
      remove("dists.env")
      dists.env <- new.env()
      assign("data.dist", as.vector(dist(coords)), envir=dists.env)
    }
  }
  ##
  ## Distances between data and prediction locations
  ## Must be here AFTER anisotropy be checked
  ##
  if(do.prediction){
    assign("d0", loccoords(coords = coords, locations = locations), envir=pred.env)
    ##
    ## checking coincident data and prediction locations
    ##
    loc.coincide <- apply(get("d0", envir=pred.env), 2, function(x){any(x < 1e-10)})
    if(any(loc.coincide))
      loc.coincide <- which(loc.coincide)
    else
      loc.coincide <- NULL
    if(!is.null(loc.coincide)){
      temp.f <- function(x, data){return(data[x < 1e-10,])}
      data.coincide <- t(apply(get("d0", envir=pred.env)[,loc.coincide, drop=FALSE],
                               2,temp.f, data=data))
    }
    else data.coincide <- NULL
    n.loc.coincide <- length(loc.coincide)    
  }
  ##
  ## Preparing prior information on beta and sigmasq
  ##
  beta.info <-
    switch(prior$beta.prior,
           fixed = list(mivm = 0, ivm = 0, iv = Inf, beta.fixed = beta.fixed, p = 0),
           flat = list(mivm = 0, ivm = 0, iv = 0, p = beta.size),
           normal = list(mivm = betares$mivm, ivm = betares$ivm, iv = betares$iv, p= 0))
  sigmasq.info <-
    switch(prior$sigmasq.prior,
           fixed = list(df.sigmasq = Inf, n0S0 = 0, sigmasq.fixed = sigmasq.fixed),
           reciprocal = list(df.sigmasq = 0, n0S0 = 0),
           uniform = list(df.sigmasq = -2, n0S0 = 0),
           sc.inv.chisq = list(df.sigmasq = df.sigmasq, n0S0 = df.sigmasq*S2.prior))
  ##
  ## ====================== PART 2 =============================
  ##                 FIXED PHI AND TAUSQ.REL
  ## ===========================================================
  ##
  phi.fixed <- prior$phi
  ##
  ## Computing parameters of the posterior for $\(\beta, \sigma^2)$ 
  ## and variables to be used for prediction (if applies)
  ##
  iR <- varcov.spatial(dists.lowertri = get("data.dist", envir=dists.env),
                       cov.model = model$cov.model,
                       kappa = model$kappa, nugget = tausq.rel.fixed,
                       cov.pars = c(1, phi.fixed), inv = TRUE,
                       only.inv.lower.diag = TRUE)
  yiRy <- diagquadraticformXAX(data, iR$lower.inverse, iR$diag.inverse)
  xiRy.x <- bilinearformXAY(X = trend.data, lowerA = iR$lower.inverse,
                            diagA = iR$diag.inverse, Y = cbind(data, trend.data))
  if(!is.matrix(xiRy.x)) xiRy.x <- is.matrix(xiRy.x, 1, n.datasets+beta.size)
  ind.datasets <- seq(length=n.datasets)
  xiRx <- xiRy.x[,-ind.datasets, drop = FALSE]
  ## 1. Computing parameters of posterior for beta
  ##
  if(any(beta.info$iv == Inf)){ 
    beta.post <- beta.info$beta.fixed
    beta.var.std.post <- matrix(0, ncol = beta.size, nrow = beta.size)
    inv.beta.var.std.post <- Inf
  }
  else{
    inv.beta.var.std.post <- as.matrix(beta.info$iv + xiRx)    
    beta.var.std.post <- solve.geoR(inv.beta.var.std.post)
    beta.post <- beta.var.std.post %*% (beta.info$ivm + xiRy.x[,ind.datasets])
  }
  ##
  ## 2. Computing parameters of posterior for sigmasq
  ##
  if(sigmasq.info$df.sigmasq == Inf){
    S2.post <- sigmasq.info$sigmasq.fixed
    df.post <- Inf
  }
  else{
    df.post <- n + sigmasq.info$df.sigmasq - beta.info$p
    ##
    if(any(beta.info$iv == Inf)){
      S2.post <- sigmasq.info$n0S0 + yiRy - 2*crossprod(beta.post,xiRy.x[,ind.datasets,drop = FALSE]) + as.vector(t(beta.post)%*%xiRx%*%beta.post)
    }
    else{
      S2.post <- sigmasq.info$n0S0 + beta.info$mivm + yiRy - diagquadraticformXAX(beta.post, inv.beta.var.std.post[lower.tri(inv.beta.var.std.post)], diag(inv.beta.var.std.post))     
    }
    S2.post <- drop(S2.post/df.post)
  } 
  ##
  ## Preparing output of the posterior distribution
  ##
  if(prior$beta.prior == "fixed") kb$posterior$beta <- list(status = "fixed", fixed.value = beta.fixed)
  else {
    if(prior$sigmasq.prior == "fixed") kb$posterior$beta <- list(distribution = "normal")
    else kb$posterior$beta <- list(distribution = "t", conditional = "normal")
    kb$posterior$beta$pars <- list(mean = beta.post, var = beta.var.std.post%o%S2.post)
  }
  if(prior$sigmasq.prior == "fixed") kb$posterior$sigmasq <- list(status="fixed", fixed.value=sigmasq.fixed)
  else kb$posterior$sigmasq <- list(distribution = "sc.inv.chisq", pars = list(df = df.post, S2 = S2.post))
  kb$posterior$phi<- list(status= "fixed", fixed.value = phi.fixed)
  kb$posterior$tausq.rel <- list(status= "fixed", fixed.value = tausq.rel.fixed)
  ##
  ## Preparing output of the predictive distribution
  ##
  if(do.prediction){
    v0 <- cov.spatial(obj = get("d0", envir=pred.env),
                      cov.model = cov.model, kappa = kappa,
                      cov.pars = c(1, phi.fixed))
    ## care here, reusing object b
    b <- bilinearformXAY(X = cbind(data, trend.data),
                         lowerA = as.vector(iR$lower.inverse),
                         diagA = as.vector(iR$diag.inverse), 
                         Y = as.vector(v0))    
    tv0ivdata <- t(b[ind.datasets, , drop=FALSE])
    b <- t(get("trend.loc", envir=pred.env)) - b[-ind.datasets, , drop=FALSE]
    ##
    if(any(beta.info$iv == Inf)) kb$predictive$mean <- tv0ivdata + as.vector(crossprod(b, beta.post))
    else kb$predictive$mean <- tv0ivdata + crossprod(b, beta.post)
    if((tausq.rel.fixed < 1e-12) & (!is.null(loc.coincide))){
      kb$predictive$mean[loc.coincide,] <- data.coincide
    }
    ##
    R.riRr.bVb <- 1 - diagquadraticformXAX(X = v0, lowerA = iR$lower.inverse, diagA = iR$diag.inverse)
    if(all(beta.info$iv != Inf))
      R.riRr.bVb <- R.riRr.bVb + diagquadraticformXAX(X = t(b), lowerA=beta.var.std.post[lower.tri(beta.var.std.post)],
                                                      diagA = diag(beta.var.std.post))
    ##
    kb$predictive$variance <- as.vector(tausq.rel.fixed + R.riRr.bVb)%*%t(S2.post)
    if(((tausq.rel.fixed < 1e-12) ) & !is.null(loc.coincide)) kb$predictive$variance[loc.coincide] <- 0
    kb$predictive$variance[kb$predictive$variance < 1e-16] <- 0
    if(sigmasq.info$df.sigmasq != Inf)
      kb$predictive$variance <- (df.post/(df.post-2))*kb$predictive$variance
    kb$predictive$distribution <- ifelse(prior$sigmasq.prior == "fixed", "normal", "t")
    remove("R.riRr.bVb","yiRy","xiRy.x","xiRx")  
  }
  ##
  ## ======================= PART 4 ==============================
  ##                Sampling from the predictive
  ## =============================================================
  ##
  if(do.prediction && simulations.predictive){
    if(is.R()){
      if(cov.model.number > 12)
        stop("simulation in krige.bayes.extnd not implemented for the choice of correlation function")
    }
    else
      if(cov.model.number > 10)
        stop("simulation in krige.bayes.extnd not implemented for the chosen correlation function")
    if(messages.screen){
      cat("krige.bayes.extnd: sampling from the predictive\n")
    }
    tmean <- kb$predictive$mean
    tv0ivdata <- NULL        ### se efter om den kan fjernes foer
    Dval <-  1.0 + tausq.rel
    coincide.cond <- any(loc.coincide)
    nloc <- ni - n.loc.coincide
    if(coincide.cond){
      ind.not.coincide <- (-loc.coincide)
      v0 <- v0[, ind.not.coincide, drop=FALSE]
      tmean <- tmean[ind.not.coincide, , drop=FALSE]
      b <- b[,ind.not.coincide, drop=FALSE]
    }
    else ind.not.coincide <- TRUE
    if(n.predictive > 1){
      warning("n.predictive > 1 is not implemented, n.predictive = 1")
      n.predictive <- 1
    }
    kb$predictive$simulations <- matrix(NA, nrow=ni, ncol=n.datasets)
    if(nloc>0){
      kb$predictive$simulations[ind.not.coincide,] <- cond.sim(env.loc = base.env, env.iter = base.env,
                                                               loc.coincide = loc.coincide,
                                                               coincide.cond = coincide.cond, tmean = tmean,
                                                               Rinv = list(lower=iR$lower.inverse, diag=iR$diag.inverse),
                                                               mod = list(beta.size = beta.size, nloc = nloc, Nsims = n.datasets, n = n,
                                                                 Dval = Dval, df.model = df.post, s2 = S2.post,
                                                                 cov.model.number = cov.model.number, phi = phi.fixed, kappa = kappa),
                                                               vbetai = beta.var.std.post,
                                                               fixed.sigmasq = (sigmasq.info$df.sigmasq == Inf))
    }
    if(coincide.cond) kb$predictive$simulations[loc.coincide,] <- data.coincide
  }
  if(!do.prediction) kb$predictive <- "no prediction locations provided"
  kb$.Random.seed <- seed
  kb$max.dist <- data.dist.max
  kb$call <- call.fc
  attr(kb, "prediction.locations") <- call.fc$locations
  return(kb)
}


"krige.conv.extnd" <- 
function(geodata, coords = geodata$coords, data = geodata$data, locations, krige, output)
{
##############################################################################
  ##     An extended version of krige.conv (geoR) allowing for multivariate data 
  ##     (useful for performing kriging on simulations: output from MCMC)
##############################################################################
  ##     data   : n*m-matrix of values of observations (m datasets)
  ##     n.predictive  : only 0 and 1 is allowed
##############################################################################
  base.env <- sys.frame(sys.nframe())
  if(missing(geodata))
    geodata <- list(coords=coords, data=data)
  cl <- match.call()
  data <- as.matrix(data)
  n.datasets <- ncol(data)
  n <- nrow(data)
  ##
  ## reading input
  ##
  if(missing(krige))
    krige <- krige.control()
  else{
    if(class(krige) != "krige.geoR"){
      if(!is.list(krige))
        stop("krige.conv.extnd: the argument krige only takes a list or an output of the function krige.control")
      else{
        krige.names <-c("type.krige","trend.d","trend.l","obj.model","beta","cov.model",
"cov.pars","kappa","nugget","micro.scale","dist.epsilon","lambda","aniso.pars")
        krige <- object.match.names(krige,krige.names)
        if(is.null(krige$type.krige)) krige$type.krige <- "ok"  
        if(is.null(krige$trend.d)) krige$trend.d <-  "cte"
        if(is.null(krige$trend.l)) krige$trend.l <-  "cte"
        if(is.null(krige$obj.model)) krige$obj.model <-  NULL
        if(is.null(krige$beta)) krige$beta <- NULL 
        if(is.null(krige$cov.model)) krige$cov.model <- "matern"  
        if(is.null(krige$cov.pars))
          stop("covariance parameters (sigmasq and phi) should be provided in cov.pars")
        if(is.null(krige$kappa)) krige$kappa <-  0.5
        if(is.null(krige$nugget)) krige$nugget <-  0
        if(is.null(krige$micro.scale)) krige$micro.scale <- 0  
        if(is.null(krige$dist.epsilon)) krige$dist.epsilon <-  1e-10
        if(is.null(krige$aniso.pars)) krige$aniso.pars <- NULL  
        if(is.null(krige$lambda)) krige$lambda <- 1 
        krige <- krige.control(type.krige = krige$type.krige,
                               trend.d = krige$trend.d, trend.l = krige$trend.l,
                               obj.model = krige$obj.model,
                               beta = krige$beta, cov.model = krige$cov.model,
                               cov.pars = krige$cov.pars, kappa = krige$kappa,
                               nugget = krige$nugget, micro.scale = krige$micro.scale,
                               dist.epsilon = krige$dist.epsilon, 
                               aniso.pars = krige$aniso.pars,
                               lambda = krige$lambda)
      }
    }
  }
  cov.model <- krige$cov.model
  kappa <- krige$kappa
  lambda <- krige$lambda
  beta <- krige$beta
  cov.pars <- krige$cov.pars
  nugget <- krige$nugget
  ##
  ## reading output options
  ##
  if(missing(output))
    output <- output.control()
  else{
    if(class(output) != "output.geoR"){
      if(!is.list(output))
        stop("krige.conv.extnd: the argument output only takes a list or an output of the function output.control")
      else{
        output.names <- c("n.posterior","n.predictive","moments","n.back.moments","simulations.predictive",
                          "mean.var","quantile","threshold","signal","messages.screen")
        output <- object.match.names(output,output.names)
        if(is.null(output$n.posterior)) output$n.posterior <- 1000 
        if(is.null(output$n.predictive)) output$n.predictive <- NULL
        if(is.null(output$moments)) output$moments <- TRUE
        if(is.null(output$n.back.moments)) output$n.back.moments <- 1000 
        if(is.null(output$simulations.predictive)){
          if(is.null(output$n.predictive)) output$simulations.predictive <- NULL
          else
            output$simulations.predictive <- ifelse(output$n.predictive > 0, TRUE, FALSE)
        }
        if(is.null(output$mean.var)) output$mean.var <- NULL
        if(is.null(output$quantile)) output$quantile <- NULL
        if(is.null(output$threshold)) output$threshold <- NULL
        if(is.null(output$signal)) output$signal <- NULL
        if(is.null(output$messages.screen)) output$messages.screen <- TRUE
        output <- output.control(n.posterior = output$n.posterior,
                                 n.predictive = output$n.predictive,
                                 moments = output$moments,
                                 n.back.moments = output$n.back.moments, 
                                 simulations.predictive = output$simulations.predictive,
                                 mean.var = output$mean.var, quantile = output$quantile,
                                 threshold = output$threshold, signal = output$signal,
                                 messages = output$messages.screen)
      }
    }
  }
  ##
  signal <- ifelse(is.null(output$signal), FALSE, output$signal)
  messages.screen <- output$messages.screen
  n.predictive <- output$n.predictive
  n.back.moments <- output$n.back.moments
  ##
  ## checking input
  ##
  if(krige$type.krige == "ok") beta.prior <- "flat"
  else beta.prior <- "deg"
  if(n.predictive > 1) {
    warning("n.predictive redefined: n.predictive=1")
    n.predictive <- 1
  }
  if(n.datasets == 1)
    stop("Only one dataset; please use krige.conv instead")
  ##
  if(is.vector(coords)) {
    coords <- cbind(coords, 0)
    warning("vector of coordinates: one spatial dimension assumed")
  }
  coords <- as.matrix(coords)
  if(is.vector(locations)) {
    if(length(locations) == 2) {
      locations <- t(as.matrix(locations))
      warning("assuming that there is 1 prediction point")
    }
    else {
      warning("vector of locations: one spatial dimension assumed")
      locations <- as.matrix(cbind(locations, 0))
    }
  }
  else locations <- as.matrix(locations)
  dimnames(coords) <- list(NULL, NULL)
  dimnames(locations) <- list(NULL, NULL)
  ##
  if(messages.screen){
    if(is.numeric(krige$trend.d))
      cat("krige.conv.extnd: model with covariates matrix provided by the user")
    else
      cat(switch(as.character(krige$trend.d)[1],
                 "cte" = "krige.conv.extnd: model with mean being constant",
                 "1st" = "krige.conv.extnd: model with mean given by a 1st order polynomial on the coordinates",
                 "2nd" = "krige.conv.extnd: model with mean given by a 2nd order polynomial on the coordinates",
                 "krige.conv.extnd: model with mean defined by covariates provided by the user"))
    cat("\n")
  }
  trend.data <- unclass(trend.spatial(trend=krige$trend.d, geodata = geodata))
  beta.size <- ncol(trend.data)
  if(nrow(trend.data) != n) stop("length of trend is different from the length of the data")
  trend.l <- unclass(trend.spatial(trend=krige$trend.l, geodata = list(coords = locations)))
  ni <- nrow(trend.l)
  if(nrow(locations) != ni) stop("length of trend is different from the number of locations for prediction")
  ##
  ## Anisotropy correction (should be placed AFTER trend.d/trend.l
  ##
  if(!is.null(krige$aniso.pars)){
    if(messages.screen) cat("krige.conv.extnd: anisotropy correction performed\n")
    coords <- coords.aniso(coords = coords, aniso.pars = krige$aniso.pars)
    locations <- coords.aniso(coords = locations, aniso.pars = krige$aniso.pars)
  }
  ##
  ## Box-Cox transformation
  ##
  if(lambda != 1) {
    if(messages.screen) cat("krige.conv.extnd: Data transformation (Box-Cox) performed.\n")
    if(lambda == 0)
      data <- log(data)
    else data <- ((data^lambda) - 1)/lambda
  }
  ##
  ## setting covariance parameters
  ##
  tausq <- nugget
  if(is.vector(cov.pars)) {
    sigmasq <- cov.pars[1]
    phi <- cov.pars[2]
  }
  else stop("covariance parametershould be given as a vector")
  sill.partial <- krige$micro.scale + sigmasq
  sill.tot <- tausq + sigmasq
  ##
  ## starting kriging calculations
  ##
  kc.result <- list()
  invcov <- varcov.spatial(coords = coords, cov.model = cov.model, kappa = kappa, nugget = nugget, cov.pars = cov.pars, 
                           inv = TRUE, only.inv.lower.diag = TRUE)
  ## From a numerical perspective there quite a few ugly bits here. See R-help e-mails from Bates and Ripley about LS/GLS. However, changing is not simple.
  ittivtt <- solve.geoR(bilinearformXAY(trend.data, invcov$lower.inverse, invcov$diag.inverse, trend.data))
  if(beta.prior == "flat") {
    beta.flat <- ittivtt %*% bilinearformXAY(trend.data, invcov$lower.inverse, invcov$diag.inverse, data)
  }
  temp <- NULL
  d0mat <- loccoords(coords, locations)
  if(krige$micro.scale != 0) {
    v0 <- ifelse(d0mat < krige$dist.epsilon, sill.partial, 
                 cov.spatial(obj = d0mat, cov.model = cov.model, kappa = kappa, cov.pars = cov.pars))
  }
  else {
    v0 <- cov.spatial(obj = d0mat, cov.model = cov.model, kappa = kappa, cov.pars = cov.pars)
  }
  tv0ivv0 <- diagquadraticformXAX(v0, invcov$lower.inverse, invcov$diag.inverse)
  b <- trend.l - bilinearformXAY(v0, invcov$lower.inverse, invcov$diag.inverse, trend.data)
  tv0ivdata <- bilinearformXAY(v0, invcov$lower.inverse, invcov$diag.inverse, data)
  if(n.predictive == 0) {
    remove(list = c("v0", "invcov"))
  }
  if(beta.prior == "deg") {
    kc.result$predict <- array(tv0ivdata + rep(as.vector(b %*% beta), n.datasets), dim = c(ni, n.datasets))
    if(n.predictive == 0) remove(list = c("b", "tv0ivdata"))
    else remove("tv0ivdata")
    if(output$signal) kc.result$krige.var <- as.vector(sill.partial - tv0ivv0)
    else kc.result$krige.var <- as.vector(sill.tot - tv0ivv0)
    beta.est <- paste("krige.conv.extnd: Simple kriging (beta provided by user)\n")
  }
  if(beta.prior == "flat") {
    kc.result$predict <- array(tv0ivdata, dim = c(ni, n.datasets)) + b %*% beta.flat
    remove(list = c("tv0ivdata"))
    if(beta.size == 1) bitb <- as.vector(b^2) * as.vector(ittivtt)
    else bitb <- diagquadraticformXAX(t(b), ittivtt[lower.tri(ittivtt)], diag(ittivtt))
    if(output$signal) kc.result$krige.var <- as.vector(sill.partial - tv0ivv0 + bitb)
    else kc.result$krige.var <- as.vector(sill.tot - tv0ivv0 + bitb)
    kc.result$beta.est <- beta.flat
    remove("beta.flat")
  }
  if(any(round(kc.result$krige.var, dig=12) < 0))
    warning("krige.conv.extnd: negative kriging variance found! Investigate why this is happening.\n")
  out.message <- "krige.conv.extnd: Kriging performed using global neighbourhood"
  if(messages.screen) cat(paste(out.message, "\n"))
############## Sampling from the resulting distribution #####################
  if(n.predictive > 0) {
    ## checking coincident data points and prediction locations
    loc.coincide <- (colSums(d0mat < krige$dist.epsilon) == 1)
    if(any(loc.coincide)) {
      if(all(loc.coincide)) stop("locations is a subset of coords; prediction not performed")
      loc.coincide <- which(loc.coincide)
    }
    else loc.coincide <- NULL
    d0mat <- NULL
    if(messages.screen) cat("krige.conv.extnd: sampling from the predictive distribution (conditional simulations)\n")    
    if(signal) Dval <- 1. + (krige$micro.scale/sigmasq)
    else Dval <- 1. + (nugget/sigmasq)
    if(beta.prior == "deg") vbetai <- matrix(0, ncol = beta.size, nrow = beta.size)
    else vbetai <- matrix(ittivtt, ncol = beta.size, nrow = beta.size)
    coincide.cond <- (((round(1e12 * nugget) == 0) | !signal) & (!is.null(loc.coincide)))
    nloc <- ni - length(loc.coincide)
    if(coincide.cond){
      ind.not.coincide <- -(loc.coincide) 
      v0 <- v0[,ind.not.coincide, drop=FALSE]
      b <- b[,ind.not.coincide, drop=FALSE]
    }
    else ind.not.coincide <- TRUE
    kc.result$simulations <- matrix(0, nrow = ni, ncol = n.datasets)
    if(nloc>0){
      kc.result$simulations[ind.not.coincide,  ] <- cond.sim(env.loc = base.env, env.iter = base.env,  loc.coincide = loc.coincide,
                                                             coincide.cond = coincide.cond,
                                                             tmean = kc.result$predict[ind.not.coincide, , drop = FALSE],
                                                             Rinv = invcov,
                                                             mod = list(beta.size = beta.size, nloc = nloc,
                                                               Nsims = n.datasets, n = n, Dval = Dval,
                                                               df.model = NULL, s2 = sigmasq,
                                                               cov.model.number = cor.number(cov.model),
                                                               phi = phi, kappa = kappa),
                                                             vbetai = vbetai, fixed.sigmasq = TRUE)
    }
    if(coincide.cond)
      kc.result$simulations[loc.coincide,  ] <- kc.result$predict[loc.coincide,, drop = FALSE]         
    remove(list = c("v0", "invcov", "b"))
  }
######################	Back-transforming predictions ############
  if(lambda != 1) {
    if(lambda == 0){
      predict.transf <- kc.result$predict
      if(messages.screen) cat("krige.conv.extnd: back-transforming the predictions using formula for EXP() \n")
      kc.result$predict <- exp(predict.transf + 0.5 * kc.result$krige.var)
      kc.result$krige.var <- (exp(2 * predict.transf - kc.result$krige.var))*expm1(kc.result$krige.var)
      remove("predict.transf")
    }
    if(lambda > 0){
      ## using second order taylor-expansion + facts for N(0,1) [third moment = 0 ; fourth moment = 12].
      if(messages.screen) cat("krige.conv.extnd: back-transforming predictions using 2. order Taylor expansion for g^{-1}() \n")
      ivBC <- BC.inv(kc.result$predict,lambda)
      kc.result$predict <- ivBC + 0.5 * ((1-lambda)*ivBC^(1-2*lambda))*kc.result$krige.var
      kc.result$krige.var <- (ivBC^(1-lambda))^2*kc.result$krige.var + (11/4)*((1-lambda)*ivBC^(1-2*lambda))^2*kc.result$krige.var^2
      remove("ivBC")
    }
    if(lambda < 0){
      if(messages.screen) cat("krige.conv.extnd: resulting distribution has no mean for lambda < 0 - back transformation not performed\n")
    }
    if(n.predictive > 0) {
      kc.result$simulations <- BC.inv(kc.result$simulations,lambda)
    }
  }
  else{
    temp <- kc.result$krige.var
    kc.result$krige.var <- matrix(0, nrow = ni, ncol = n.datasets)
    kc.result$krige.var[] <- temp
  }
  kc.result <- c(kc.result, list(message = message, call = cl))
#######################################
  return(kc.result)
}

"pmixed" <- 
  function(value, parms, df)
{
  if(ncol(parms$var)==1) parms$var <- as.vector(parms$var)
  if(df == Inf){
    sd.error <- sqrt(parms$var)
    if(any(sd.error < 1e-12)) sd.error[sd.error < 1e-12] <- 1e-12
    temp <- array(pnorm((value - parms$mean)/sd.error), dim = c(nrow(parms$mean), ncol(parms$mean)))
  }
  else{
    sd.error <- sqrt(parms$var*(df - 2)/df)
    if(any(sd.error < 1e-12)) sd.error[sd.error < 1e-12] <- 1e-12
    temp <- array(pt((value - parms$mean)/sd.error, df = df), dim = c(nrow(parms$mean), ncol(parms$mean)))
  }
  temp2 <- rowMeans(temp)
  return(temp2)
}

"model.glm.control" <- 
  function(trend.d = "cte", trend.l = "cte", cov.model = "matern", kappa = 0.5, aniso.pars = NULL, lambda = 0)
{
  cov.model <- match.arg(cov.model,
              choices = c("matern", "exponential", "gaussian",
                "spherical", "circular", "cubic", "wave", "power",
                "powered.exponential", "cauchy", "gneiting",
                "gneiting.matern", "pure.nugget"))
  if(cov.model == "powered.exponential" & (kappa <= 0 | kappa > 2))
    stop("model.glm.control: for power exponential correlation model the parameter kappa must be in the interval \(0,2\]")
  if(cov.model == "power") stop("model.glm.control: correlation function does not exist for the power variogram")
  if(!is.null(aniso.pars)){ 
    if(length(aniso.pars) != 2 | !is.numeric(aniso.pars))
      stop("anisotropy parameters must be a vector with two elements: rotation angle (in radians) and anisotropy ratio (a number > 1)")
  }
  res <- list(trend.d = trend.d, trend.l = trend.l, cov.model = cov.model, kappa = kappa, aniso.pars = aniso.pars, lambda = lambda)
  class(res) <- "model.geoRglm"
  return(res)
}


"model.glm.check.aux" <-
  function(model, fct)
{
  if(class(model) != "model.geoRglm"){
    if(!is.list(model))
      stop(paste(fct,": the argument model only takes a list or an output of the function model.glm.control"))
    else{
      model.names <- c("trend.d", "trend.l", "cov.model", "kappa", "aniso.pars", "lambda")
      model <- object.match.names(model,model.names)
      if(is.null(model$trend.d)) model$trend.d <- "cte"  
      if(is.null(model$trend.l)) model$trend.l <- "cte"  
      if(is.null(model$cov.model)) model$cov.model <- "matern"  
      if(is.null(model$kappa)) model$kappa <- 0.5
      if(is.null(model$lambda)){
        if(fct=="pois.krige.bayes") model$lambda <- 0
        if(fct=="binom.krige.bayes") model$lambda <- NULL
      }
      model <- model.glm.control(trend.d = model$trend.d,
                                 trend.l = model$trend.l,
                                 cov.model = model$cov.model,
                                 kappa = model$kappa,
                                 aniso.pars = model$aniso.pars,
                                 lambda = model$lambda)
    }
  }
  return(model)
}


"multgauss" <- 
  function(cov)
{
  if(is.R())
    return(crossprod(chol(cov), rnorm(n=ncol(cov))))
  else
    return(rmvnorm(ncol(cov), cov = cov))
}

"mcmc.control" <- 
  function(S.scale, Htrunc="default", S.start, burn.in=0, thin=10, n.iter=1000*thin, phi.start="default",  phi.scale=NULL)
{
  if(missing(S.scale)) stop("S.scale parameter must to be provided for MCMC-proposal")
  if(missing(S.start) || is.null(S.start)) S.start<- "default"
  if(is.null(Htrunc)) Htrunc <- "default"
  if(is.null(burn.in)) burn.in <- 0
  if(is.null(thin)) thin <- 10
  if(is.null(n.iter)) n.iter <- 1000*thin
  res <- list(S.scale = S.scale, Htrunc = Htrunc, S.start = S.start, burn.in = burn.in, thin = thin, n.iter = n.iter,
              phi.start = phi.start, phi.scale=phi.scale)
  class(res) <- "mcmc.geoRglm"
  return(res)
}

"mcmc.check.aux" <-
  function(mcmc.input, fct)
{
  if(class(mcmc.input) != "mcmc.geoRglm"){
    if(!is.list(mcmc.input))
      stop(paste(fct,": the argument mcmc.input only takes a list or an output of the function mcmc.control"))
    else{
      mcmc.input.names <- c("S.scale", "Htrunc", "S.start", "burn.in", "thin", "n.iter", "phi.start", "phi.scale")  
      mcmc.input <- object.match.names(mcmc.input,mcmc.input.names)
      if(fct=="pois.krige.bayes" | fct=="binom.krige.bayes"){
        if(is.null(mcmc.input$phi.start)) mcmc.input$phi.start <- "default"
      }
      mcmc.input <- mcmc.control(S.scale = mcmc.input$S.scale,Htrunc=mcmc.input$Htrunc,S.start=mcmc.input$S.start,
                                 burn.in=mcmc.input$burn.in,thin=mcmc.input$thin,n.iter=mcmc.input$n.iter,
                                 phi.start=mcmc.input$phi.start,phi.scale=mcmc.input$phi.scale)
    }
  }
  return(mcmc.input)
}


"BC.inv" <- 
  function(z,lambda)
{
  return(BCtransform(z, lambda = lambda, inverse=TRUE)$data)
}


"output.glm.control" <-
  function(sim.posterior, sim.predict, keep.mcmc.sim, quantile, threshold, inference, messages)
{
  ##
  ## Assigning default values
  ##
  if(missing(sim.posterior) || is.null(sim.posterior)) sim.posterior <- TRUE
  if(missing(sim.predict) || is.null(sim.predict)) sim.predict <- FALSE
  if(missing(keep.mcmc.sim) || is.null(keep.mcmc.sim)) keep.mcmc.sim <- TRUE
  if(missing(quantile) || is.null(quantile)) quantile.estimator <- TRUE
  else quantile.estimator <- quantile
  if(missing(threshold)) probability.estimator <- NULL 
  else probability.estimator <- threshold
  if(missing(inference) || is.null(inference)) inference <- TRUE
  if(missing(messages))
    messages.screen <- ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages"))
  else messages.screen <- messages
  ##
  if(is.null(quantile.estimator)) quantile.estimator <- TRUE
  else{
    if(is.numeric(quantile.estimator))
      if(any(quantile.estimator) < 0 | any(quantile.estimator) > 1)
        stop("quantiles indicators must be numbers in the interval [0,1]\n")
    if(!inference) {
      warning("prediction not performed; quantile.estimator is set to NULL \n")
      quantile.estimator <- NULL
    }
  }
  if(!is.null(probability.estimator)){
    if(!is.numeric(probability.estimator))
      stop("threshold must be a numeric value (or vector) of cut-off value(s)\n")
    if(!inference) {
      warning("prediction not performed; probabilitites above a threshold cannot be computed\n")
      probability.estimator <- NULL
    }
  }
  res <- list(sim.posterior = sim.posterior, sim.predict = sim.predict, keep.mcmc.sim = keep.mcmc.sim,
              quantile.estimator = quantile.estimator, probability.estimator = probability.estimator,
              inference = inference, messages.screen = messages.screen)
  class(res) <- "output.geoRglm"
  return(res)
}


"output.glm.check.aux" <-
  function(output, fct)
{
  if(class(output) != "output.geoRglm"){
    if(!is.list(output))
      stop(paste(fct,": the argument output only takes a list or an output of the function output.glm.control"))
    else{
      output.names <- c("sim.posterior","sim.predict", "keep.mcmc.sim","quantile","threshold","inference","messages.screen")
      output <- object.match.names(output,output.names)
      if(is.null(output$sim.predict)) output$sim.predict <- FALSE
      if(is.null(output$messages.screen)) output$messages.screen <- TRUE        
      output <- output.glm.control(sim.posterior = output$sim.posterior,
                                   sim.predict = output$sim.predict,
                                   keep.mcmc.sim = output$keep.mcmc.sim, quantile = output$quantile,
                                   threshold = output$threshold, inference = output$inference,
                                   messages = output$messages.screen)
    }
  }
  return(output)
}


"object.match.names" <- function(obj,names)
{
  obj.names <- names(obj)
  new.obj <- list()
  for(i in seq(along=obj)){
    n.match <- match.arg(obj.names[i], names)
    new.obj[[n.match]] <- obj[[i]]
  }
  return(new.obj)
}
"geoRglmdefunct" <- function()
{
  cat("\n")
  cat("The following functions are no longer used in geoRglm:")
  cat("---------------------------------------------------")
  cat(" pois.log.krige: use pois.krige instead")
  cat(" y50: use p50 instead")
  cat("\n")
}

cite.geoRglm <- function()
{
  cat("\n")
  cat("To cite geoR in publications, use\n\n")
  msg <- "CHRISTENSEN, O.F. & RIBEIRO Jr., P.J. (2002) geoRglm: A package for generalised linear spatial models. R-NEWS, Vol 2, No 2, 26-28. ISSN 1609-3631."
  writeLines(strwrap(msg, prefix = "  "))
  cat("\n")
  msg <- paste("Please cite geoRglm when using it for data analysis!")
  writeLines(strwrap(msg))
  cat("\nA BibTeX entry for LaTeX users is\n\n")
  cat("  @Article{,\n")
  cat("     title	   = {Christensen, O.F. and Ribeiro Jr., P.J.},\n")
  cat("     author        = {{geoRglm}: A package for generalised linear spatial models},\n")
  cat("     journal       = {R-NEWS},\n")
  cat("     year	   = {2002},\n")
  cat("     volume	   = {2},\n")
  cat("     number	   = {2},\n")
  cat("     pages	   = {26--28},\n")
  cat("     issn          = {1609-3631},\n")
  cat("     url           = {http://cran.R-project.org/doc/Rnews}\n")
  cat("   }\n\n")
}

"glsm.krige" <- function(mcmc.output, locations, borders = NULL, trend.l="cte", micro.scale=NULL, dist.epsilon= 1e-10,  output)
{
  call.fc <- match.call()
  coords <- mcmc.output$geodata$coords
  cov.model <- mcmc.output$model$cov.model
  kappa <- mcmc.output$model$kappa
  beta <- mcmc.output$model$beta
  cov.pars <- mcmc.output$model$cov.pars
  nugget <- mcmc.output$model$nugget
  aniso.pars <- mcmc.output$model$aniso.pars
  trend.d <- mcmc.output$model$trend
  lambda <- mcmc.output$model$lambda
  if(is.null(micro.scale)) micro.scale <- nugget
  else if(!is.numeric(micro.scale) || length(micro.scale)>1) stop("micro.scale must be a numeric number ")
  ##
  if(missing(output)) output <- output.glm.control()
  else output <- output.glm.check.aux(output, fct = "glsm.krige")
  sim.predict <- output$sim.predict
  messages.screen <- output$messages.screen
  if(messages.screen) cat(" glsm.krige: Prediction for a generalised linear spatial model \n")
  ##
  ## Checking for 1D prediction
  if(missing(locations))
    stop("locations need to be specified for prediction; prediction not performed")
  else {
    if(is.null(trend.l))
      stop("trend.l needed for prediction")
  } 
  if(length(unique(locations[,1])) == 1 | length(unique(locations[,2])) == 1)
    krige1d <- TRUE
  else krige1d <- FALSE
  ##
  beta.size <- length(beta)
  if(beta.size > 1) beta.names <- paste("beta", (0:(beta.size-1)), sep="")
  else beta.names <- "beta"
  ##
  ##------------------------------------------------------------
######################## ---- prediction ----- #####################
  if(!is.null(borders)){
    locations <- locations.inside(locations, borders)
    if(nrow(locations) == 0)
      stop(" glsm.krige : there are no prediction locations inside the borders")
    if(messages.screen)
      cat(" glsm.krige: results will be returned only for prediction locations inside the borders\n")
  }
  if(mcmc.output$model$family=="binomial"){
    krige <- list(type.krige = "sk", beta = beta, trend.d = trend.d, trend.l = trend.l, cov.model = cov.model, 
                  cov.pars = cov.pars, kappa = kappa, nugget = nugget, micro.scale = micro.scale, dist.epsilon = dist.epsilon, 
                  aniso.pars = aniso.pars, link = mcmc.output$model$link)
    kpl.result <- glm.krige.aux(data = mcmc.output$simulations, coords = coords, locations = locations, krige = krige,
                                output = list(n.predictive = ifelse(sim.predict,1,0),
                                  signal = TRUE, messages=FALSE))
  }
  else{
    if(mcmc.output$model$family!="poisson") stop("only poisson and binomial are allowed as error distribution")
    krige <- list(type.krige = "sk", beta = beta, trend.d = trend.d, trend.l = trend.l, cov.model = cov.model, 
                  cov.pars = cov.pars, kappa = kappa, nugget = nugget, micro.scale = micro.scale, dist.epsilon = dist.epsilon, 
                  aniso.pars = aniso.pars, lambda = lambda)
    kpl.result <- krige.conv.extnd(data = BC.inv(mcmc.output$simulations, lambda), coords = coords, locations = locations, krige = krige,
                                   output = list(n.predictive = ifelse(sim.predict,1,0), signal = TRUE, messages = FALSE))
    
  }
  ##
  kpl.result$krige.var <- rowMeans(kpl.result$krige.var) + apply(kpl.result$predict, 1, var)
  kpl.result$mcmc.error <- sqrt(asympvar(kpl.result$predict,messages=FALSE)/ncol(kpl.result$predict))
  kpl.result$predict <- rowMeans(kpl.result$predict)
  kpl.result$beta <- NULL
  kpl.result$call <- call.fc
#######################################
  attr(kpl.result, "prediction.locations") <- call.fc$locations
  if(!is.null(locations)) attr(kpl.result, 'sp.dim') <- ifelse(krige1d, "1d", "2d")
  if(!is.null(call.fc$borders)) attr(kpl.result, "borders") <- call.fc$borders
  class(kpl.result) <- "kriging"
  return(kpl.result)
}

"prepare.likfit.glsm" <-
  function(mcmc.output, use.intensity = FALSE)
{
### this is for the glsm.mcmc() function
  ##
  if(class(mcmc.output) != "glsm.mcmc") stop("mcmc.output must be an object of class ``glsm.mcmc'' ")
  n.dat <- nrow(mcmc.output$simulations)
  n.sim <- ncol(mcmc.output$simulations)
  if(use.intensity){
    if(mcmc.output$model$family != "poisson") stop("use.intensity = TRUE is only allowed for the Poisson error distribution")
    if(any(mcmc.output$geodata$data == 0)) stop("use.intensity = TRUE is only allowed when all data are positive ")
  }
  lambda <- mcmc.output$model$lambda
  S <- mcmc.output$simulations
  trend.data <- unclass(trend.spatial(trend = mcmc.output$model$trend, geodata=mcmc.output$geodata))
  beta.size <- ifelse(is.matrix(trend.data), ncol(trend.data), 1)
  cov.model <- mcmc.output$model$cov.model
  kappa <- mcmc.output$model$kappa
  beta <- mcmc.output$model$beta
  sigmasq <- mcmc.output$model$cov.pars[1]
  phi <- mcmc.output$model$cov.pars[2]
  nugget.rel <- mcmc.output$model$nugget/sigmasq
  if(is.null(mcmc.output$model$aniso.pars)) coords <- mcmc.output$geodata$coords
  else coords <- coords.aniso(coords = mcmc.output$geodata$coords, aniso.pars = mcmc.output$model$aniso.pars)
  invcov <- varcov.spatial(coords = coords, cov.model = cov.model, kappa = kappa, nugget = nugget.rel, 
                           cov.pars = c(1,phi), det = TRUE, only.inv.lower.diag = TRUE)
  SivS <- diagquadraticformXAX(S, invcov$lower.inverse, invcov$diag.inverse) 
  if(beta.size == 1){
    DivD <- as.vector(bilinearformXAY(trend.data,invcov$lower.inverse,invcov$diag.inverse,trend.data)) 
    SivD <- as.vector(bilinearformXAY(S, invcov$lower.inverse, invcov$diag.inverse, trend.data))
    log.f.sim <- -invcov$log.det.to.half - 0.5*n.dat*log(sigmasq) - 0.5*(SivS-2*SivD*beta+DivD*beta^2)/sigmasq
  }
  else{
    DivD <- bilinearformXAY(trend.data,invcov$lower.inverse,invcov$diag.inverse,trend.data)
    SivD <- bilinearformXAY(S, invcov$lower.inverse, invcov$diag.inverse, trend.data)
    log.f.sim <-  -invcov$log.det.to.half - 0.5*n.dat*log(sigmasq) - 0.5*as.vector(SivS-2*as.vector(SivD%*%beta)+t(beta)%*%DivD%*%beta)/sigmasq
  }
  if(use.intensity){
    if(lambda==0){
      logJ <- -apply(S,2,sum)
      mu <- exp(S)
    }
    else {
      logJ <- apply(log(S*lambda+1),2,sum)*(lambda-1)/lambda
      mu <- (S*lambda+1)^(1/lambda)
    }
    return(list(family=mcmc.output$model$family, link=mcmc.output$model$link, mu = mu, trend = mcmc.output$model$trend, coords = mcmc.output$geodata$coords,
                aniso.pars = mcmc.output$model$aniso.pars, lambda = lambda, log.f.sim = log.f.sim + logJ))
  }
  else{
    return(list(family=mcmc.output$model$family, link=mcmc.output$model$link, S = S, coords = mcmc.output$geodata$coords, trend = mcmc.output$model$trend,
                aniso.pars = mcmc.output$model$aniso.pars, lambda = lambda, log.f.sim = log.f.sim))
  }
}

"func.val" <-
  function(SivS, SivD, DivD, beta, sigmasq, log.f)
{
  beta.size <- length(beta)
  if(beta.size == 1) ff <- exp(-0.5*(SivS-2*SivD*beta+DivD*beta^2)/sigmasq-log.f)
  else ff <- exp(-0.5*as.vector(SivS-2*as.vector(SivD%*%beta)+t(beta)%*%DivD%*%beta)/sigmasq-log.f)
  return(ff)
}

"NewtonRhapson.step" <-
  function(SivS, SivD, DivD, SivDi2, ff, n.dat, beta, sigmasq, steplen)
{
  beta.size <- length(beta)
  n.sim <- length(SivS)
  n.dat <- length(SivS)
  meanff <-  mean(ff)
  ## removed 1/(sigmasq^(n.dat/2) everywhere
  if(beta.size == 1){
    F1 <- mean((SivD-DivD*beta)*ff)/sigmasq
    SivDbeta <- SivD*beta
    betaDivDbeta <- DivD*beta^2
    F11 <- (mean(SivDi2*ff)-(beta*DivD)^2*meanff)/(sigmasq^2)-2*beta*DivD*F1/sigmasq-meanff*DivD/sigmasq
    F12 <- -(n.dat/2+1)*F1/sigmasq+0.5*mean((SivD-DivD*beta)*(SivS-2*SivDbeta+betaDivDbeta)*ff)/(sigmasq^3)
  }
  else{
    F1 <- (colMeans(SivD*ff)-(DivD%*%beta)*meanff)/sigmasq
    SivDbeta <- as.vector(SivD%*%beta)
    betaDivDbeta <- as.vector(t(beta)%*%DivD%*%beta)
    F11 <- colMeans(SivDi2*ff)/(sigmasq^2)-(DivD%*%beta)%*%t(DivD%*%beta)*meanff/(sigmasq^2)-((DivD%*%beta)%*%t(F1)+F1%*%t(DivD%*%beta))/sigmasq-meanff*DivD/sigmasq
    F12 <- -(n.dat/2+1)*F1/sigmasq+0.5*(colMeans(SivD*(SivS-2*SivDbeta+betaDivDbeta)*ff)-(DivD%*%beta)*mean((SivS-2*SivDbeta+betaDivDbeta)*ff))/(sigmasq^3)
  }
  F2 <- -(n.dat/2)*meanff/sigmasq + 0.5*mean((SivS-2*SivDbeta+betaDivDbeta)*ff)/(sigmasq^2)
  F22 <- (n.dat/2)*(n.dat/2+1)*meanff/(sigmasq^2)-0.5*(n.dat+2)*mean((SivS-2*SivDbeta+betaDivDbeta)*ff)/(sigmasq^3)+0.25*mean((SivS-2*SivDbeta+betaDivDbeta)^2*ff)/(sigmasq^4)
  Delta2 <- rbind(cbind(F11,F12),c(t(F12),F22))
  if(det(Delta2) != 0){
    stepvec <- solve(Delta2,c(F1,F2))
    betanew <- beta - steplen*stepvec[seq(length=beta.size)]
    sigmasqnew <- sigmasq - steplen*stepvec[beta.size+1]
    flat.message <- FALSE
  }
  else{
    flat.message <- TRUE
    betanew <- beta 
    sigmasqnew <- sigmasq 
  } 
  return(list(betanew=betanew, sigmasqnew=sigmasqnew, flat = flat.message))  
}

"maxim.aux1" <-
  function(S, invcov, trend, log.f.sim, messages.screen=FALSE)
{
  n.sim <- ncol(S)
  n.dat <- nrow(S)  
  if(is.matrix(trend)) beta.size <- ncol(trend)
  else beta.size <- 1
  SivS <- diagquadraticformXAX(S,invcov$lower.inverse,invcov$diag.inverse)
  if(beta.size == 1){
    SivD <- as.vector(bilinearformXAY(S,invcov$lower.inverse,invcov$diag.inverse,trend))
    DivD <- as.vector(bilinearformXAY(trend,invcov$lower.inverse,invcov$diag.inverse,trend))
    beta.hat <- SivD/DivD
    sigmasq.hat <- diagquadraticformXAX(S-trend%*%t(beta.hat),invcov$lower.inverse,invcov$diag.inverse)/n.dat
    SivDi2 <- SivD^2
    corr1 <- mean(-0.5*(SivS-2*SivD*mean(beta.hat)+DivD*mean(beta.hat)^2)/mean(sigmasq.hat)-log.f.sim)
    beta <- beta.hat[1]
    sigmasq <- sigmasq.hat[1]
  }
  else{
    SivD <- bilinearformXAY(S,invcov$lower.inverse,invcov$diag.inverse,trend)
    DivD <- bilinearformXAY(trend,invcov$lower.inverse,invcov$diag.inverse,trend)
    beta.hat <- t(solve(DivD,t(SivD))) ####### might be improved (GLS)
    sigmasq.hat <- diagquadraticformXAX(S-trend%*%t(beta.hat),invcov$lower.inverse,invcov$diag.inverse)/n.dat
    "cp" <- function(x){return(x%*%t(x))}
    SivDi2 <- array(t(apply(SivD,1,cp)),dim=c(n.sim,beta.size,beta.size))
    corr1 <- mean(-0.5*(SivS-2*as.vector(SivD%*%colMeans(beta.hat))+t(colMeans(beta.hat))%*%DivD%*%colMeans(beta.hat))/mean(sigmasq.hat)-log.f.sim)
    beta <- beta.hat[1,]
    sigmasq <- sigmasq.hat[1]
  }
  log.f.c <- log.f.sim + corr1
  log.hh <- (-Inf)
  for(ll in seq(length=ncol(S))){
    if(beta.size == 1) ff <- func.val(SivS, SivD, DivD, beta.hat[ll], sigmasq.hat[ll], log.f.c)
    else ff <- func.val(SivS, SivD, DivD, beta.hat[ll,], sigmasq.hat[ll], log.f.c)
    log.hhnew <- log(mean(ff))-(n.dat/2)*log(sigmasq.hat[ll])
    if(log.hhnew>log.hh){
      if(beta.size == 1) beta <-beta.hat[ll]
      else beta <- beta.hat[ll,]
      sigmasq <- sigmasq.hat[ll]
      log.hh <- log.hhnew 
    }
  }
  ## new correction (again for numerical purposes).
  if(beta.size == 1) corr2 <- mean(-0.5*(SivS-2*SivD*beta+DivD*beta^2)/sigmasq-log.f.sim)
  else corr2 <- mean(-0.5*(SivS-2*as.vector(SivD%*%beta)+t(beta)%*%DivD%*%beta)/sigmasq-log.f.sim)
  log.f.c <- log.f.sim + corr2
  ff <- func.val(SivS, SivD, DivD, beta, sigmasq, log.f.c)
  log.hh <- log(mean(ff))-(n.dat/2)*log(sigmasq)
  ##
  test <- 1
  test2 <- 1
  steplen <- 1
  while(test>0.0000000000001 | test2 > 0 ){
    New <- NewtonRhapson.step(SivS, SivD, DivD, SivDi2, ff, n.dat, beta, sigmasq, steplen)
    ffnew <- func.val(SivS, SivD, DivD, New$beta, New$sigmasq, log.f.c)
    log.hhnew <- log(mean(ffnew))-(n.dat/2)*log(New$sigmasq)
    if(New$flat & messages.screen){
      cat(paste("Problems when optimising w.r.t. beta and sigmasq: likelihood is very flat \n"))
      cat(paste("likelihood value at this stage is = ",log.hhnew+corr2,"\n"))
    }
    test <- sum((beta-New$beta)^2)+(sigmasq-New$sigmasq)^2
    test2 <- log.hhnew-log.hh
    if(test2>0){
      log.hh <- log.hhnew   
      beta <- New$beta 
      sigmasq <- New$sigmasq
      ff <- ffnew  
    }
    else{
      steplen <- steplen/2 
    }
  }
  return(list(beta = beta, sigmasq = sigmasq, logh = log.hh+corr2)) 
}


"lik.sim" <- function(pars, fp, ip, temp.list)
{ 
  ## Obligatory parameter:
  phi <- pars[1]
  if(ip$f.tausq.rel) tausq.rel <- fp$tausq.rel
  else tausq.rel <- pars[2]
  messages.screen <- ifelse(is.null(temp.list$messages.screen), TRUE,temp.list$messages.screen)
  if(messages.screen) cat(paste("phi = ",phi, "tausq.rel = ",tausq.rel,"\n"))
  ##
  ## Computing likelihood
  ##
  iv <- varcov.spatial(dists.lowertri = as.vector(dist(temp.list$coords)), cov.model = temp.list$cov.model, kappa = temp.list$kappa,
                       nugget = tausq.rel, cov.pars = c(1, phi), only.inv.lower.diag = TRUE, det = TRUE)
  negloglik <- (iv$log.det.to.half - maxim.aux1(S=temp.list$z, invcov=iv, trend = temp.list$xmat, log.f.sim = temp.list$log.f.sim, messages.screen=messages.screen)$logh)
  if(messages.screen) cat(paste("log-likelihood = ",-negloglik,"\n"))
  return(negloglik)
}

"lik.sim.boxcox" <-
  function(pars, fp, ip, temp.list)
{ 
### Function for finding m.l.e. for a given phi based on samples from mu=g^{-1}(S) ###############
### This function is only valid when all observations are positive.
  ##
  ## Obligatory parameters:
  phi <- pars[1]
  if(ip$f.tausq.rel) tausq.rel <- fp$tausq.rel
  else tausq.rel <- pars[2]
  if(ip$f.lambda) lambda <- fp$lambda
  else lambda <- pars[length(pars)]
  ##
  mu <- temp.list$mu
  ##
  messages.screen <- ifelse(is.null(temp.list$messages.screen), TRUE,temp.list$messages.screen)
  if(messages.screen) cat(paste("phi = ",phi, "tausq.rel = ",tausq.rel, "lambda= ", lambda,"\n"))
  ##
  ## computing the determinant of the transformation
  ##
  log.J.lambda <- apply(log(mu),2,sum)*(lambda-1)
  if(lambda ==0) mu <- log(mu)
  else mu <- (mu^lambda-1)/lambda
  ##
  ## Computing likelihood
  ##
  iv <- varcov.spatial(dists.lowertri = as.vector(dist(temp.list$coords)), cov.model = temp.list$cov.model, kappa = temp.list$kappa,
                       nugget = tausq.rel, cov.pars = c(1, phi), only.inv.lower.diag = TRUE, det = TRUE)
  negloglik <- (iv$log.det.to.half - maxim.aux1(S=mu, invcov=iv, trend = temp.list$xmat, log.f.sim = temp.list$log.f.sim-log.J.lambda, messages.screen=messages.screen)$logh)
  if(messages.screen) cat(paste("log-likelihood = ",-negloglik,"\n"))
  return(negloglik)
}


"likfit.glsm" <-
  function (mcmc.obj, trend = mcmc.obj$trend,
            cov.model = "matern", 
            kappa = 0.5, ini.phi, fix.nugget.rel = FALSE, nugget.rel = 0, aniso.pars = NULL, 
            fix.lambda = TRUE, lambda = NULL, limits = pars.limits(), messages, ...)
{
  ##
  ## Checking input
  ##
  geodata <- list(coords=mcmc.obj$coords)
  call.fc <- match.call()
  temp.list <- list()
  if(missing(messages))
    messages.screen <- ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages"))
  else messages.screen <- messages
  ##
  cov.model <- match.arg(cov.model,
                         choices = c("matern", "exponential","gaussian",
                           "spherical", "circular", "cubic",
                           "wave", "power",
                           "powered.exponential", "cauchy", "gneiting",
                           "gneiting.matern", "pure.nugget"))
  if(!is.null(kappa)){
    if(cov.model == "matern" & kappa == 0.5) cov.model <- "exponential"
  }
  ##
  if(is.null(mcmc.obj$S)){
    if(is.null(mcmc.obj$mu)) stop("mcmc.obj should include either an object mu or an object S.")
    if(is.null(lambda)) stop("must specify lambda")
    n <- temp.list$n <- nrow(mcmc.obj$mu)
    temp.list$mu <- mcmc.obj$mu
    if(mcmc.obj$family == "poisson") another.boxcox <- TRUE
    else stop("estimation of lambda is only possible for Poisson error distribution and boxcox link function")
  }
  else{
    if(!is.null(lambda)) warning("cannot use argument lambda with the given objects in mcmc.obj")
    if(!fix.lambda){
      warning("cannot estimate lambda from the given objects in mcmc.obj")
      fix.lambda <- TRUE
    }
    n <- temp.list$n <- nrow(mcmc.obj$S)
    temp.list$z <- mcmc.obj$S 
    another.boxcox <- FALSE
  }
  if (n != nrow(mcmc.obj$coords)) stop("Number of locations does not match with number of data")
  if(!is.null(aniso.pars)){
    if(length(aniso.pars) != 2 | !is.numeric(aniso.pars))
      stop("anisotropy parameters must be a vector with two elements: rotation angle (in radians) and anisotropy ratio (a number > 1)")
    coords <- coords.aniso(coords = as.matrix(mcmc.obj$coords), aniso.pars = aniso.pars)
  }
  else{
    if(!is.null(mcmc.obj$aniso.pars)){
      coords <- coords.aniso(coords = as.matrix(mcmc.obj$coords), aniso.pars = mcmc.obj$aniso.pars)
      aniso.pars <- mcmc.obj$aniso.pars
    }
    else coords <- as.matrix(mcmc.obj$coords)
  }
  temp.list$xmat <- unclass(trend.spatial(trend = trend, geodata=geodata))
  beta.size <- temp.list$beta.size <- dim(temp.list$xmat)[2]
  ##
  temp.list$coords <- coords
  temp.list$cov.model <- cov.model
  temp.list$kappa <- kappa
  if(cov.model=="pure.nugget"){
    ini.phi <- 1
    if(!fix.nugget.rel) nugget.rel <- 0
    fix.nugget.rel <- TRUE
  }
  else if(missing("ini.phi")) stop("likfit.glsm : must specify ini.phi ")
  ini <- ini.phi
  lower.optim <- c(limits$phi["lower"])
  upper.optim <- c(limits$phi["upper"])
  fixed.values <- list()
  ##
  if(fix.nugget.rel) {
    fixed.values$tausq.rel <- nugget.rel
  }
  else {
    ini <- c(ini, nugget.rel)
    lower.optim <- c(lower.optim, limits$tausq.rel["lower"])
    upper.optim <- c(upper.optim, limits$tausq.rel["upper"])
  }
  if(another.boxcox){
    if(fix.lambda) {
      fixed.values$lambda <- lambda
    }
    else {
      ini <- c(ini, lambda)
      lower.optim <- c(lower.optim, limits$lambda["lower"])
      upper.optim <- c(upper.optim, limits$lambda["upper"])
    }
    ip <- list(f.tausq.rel = fix.nugget.rel, f.lambda = fix.lambda)
  }
  else ip <- list(f.tausq.rel = fix.nugget.rel)
  names(ini) <- NULL
  temp.list$log.f.sim <- mcmc.obj$log.f.sim
  temp.list$messages.screen <- messages.screen
  ##
  npars <- beta.size + 1 + ifelse(cov.model == "pure.nugget",0,1) + sum(!unlist(ip))
  if(cov.model != "pure.nugget" | !fix.lambda){
    if(messages.screen){
      cat("--------------------------------------------------------------------\n")
      cat("likfit.glsm: likelihood maximisation using the function optim.\n") 
    }
    if(another.boxcox){
      lik.optim <- optim(par = ini, fn = lik.sim.boxcox, method = "L-BFGS-B",lower = lower.optim, upper = upper.optim,
                         fp = fixed.values, ip = ip, temp.list = temp.list)
    }
    else{
      lik.optim <- optim(par = ini, fn = lik.sim, method = "L-BFGS-B",lower = lower.optim, upper = upper.optim,
                         fp = fixed.values, ip = ip, temp.list = temp.list)
    }
    ##
    if(messages.screen) 
      cat("likfit.glsm: end of numerical maximisation.\n")
    par.est <- lik.optim$par
    phi <- par.est[1]
    ##
    ## Values of the maximised likelihood
    ##
    loglik.max <-  - lik.optim$value
    ##
    ## Assigning values for estimated parameters
    ##
    if(!fix.nugget.rel){
      nugget.rel <- par.est[2]
    }
    if(!fix.lambda){
      lambda <- par.est[length(par.est)]
    }
  }
  if(cov.model == "pure.nugget") phi <- NA
  ##
  gc(verbose = FALSE)  
  ##
  ## Computing estimated beta and sigmasq
  if((is.na(phi) | phi < 1e-12))
    siv <- list(diag.inverse = rep(1/(1+nugget.rel), n), lower.inverse = rep(0,n*(n-1)/2))
  else{
    siv <- varcov.spatial(coords = coords, cov.model = cov.model, kappa = kappa, nugget = nugget.rel, cov.pars = c(1, phi),
                          only.inv.lower.diag = TRUE)
  }  
  if(another.boxcox){
    if(lambda == 0){
      result <- maxim.aux1(S = log(mcmc.obj$mu), invcov = siv, trend = temp.list$xmat, log.f.sim = temp.list$log.f.sim - apply(log(mcmc.obj$mu),2,sum)*(lambda-1), messages.screen=messages.screen)
    }
    else{
      result <- maxim.aux1(S = (mcmc.obj$mu^lambda-1)/lambda, invcov = siv, trend = as.vector(temp.list$xmat),
                           log.f.sim = temp.list$log.f.sim - apply(log(mcmc.obj$mu),2,sum)*(lambda-1), messages.screen=messages.screen)
    }
    if(cov.model == "pure.nugget") loglik.max <- result$logh
    results <- list(family=mcmc.obj$family, link="boxcox", cov.model = cov.model, beta = result$beta, cov.pars=c(result$sigmasq, phi), nugget.rel = nugget.rel, 
                    kappa = kappa, aniso.pars = aniso.pars, lambda = lambda, trend = trend, npars=npars, loglik = loglik.max, call = call.fc)
  }
  else{
    result <- maxim.aux1(S = mcmc.obj$S, invcov = siv, trend = temp.list$xmat,log.f.sim = temp.list$log.f.sim, messages.screen=messages.screen)
    if(cov.model == "pure.nugget") loglik.max <- result$logh
    results <- list(family=mcmc.obj$family, link=mcmc.obj$link, cov.model = cov.model, beta = result$beta, cov.pars = c(result$sigmasq, phi), nugget.rel = nugget.rel,
                    kappa = kappa, aniso.pars = aniso.pars, lambda = mcmc.obj$lambda, trend = trend, npars=npars, loglik = loglik.max, call = call.fc)
  }
  par.su <- data.frame(status=rep(-9, beta.size + 4))
  par.su$status <- c(rep("estimated", beta.size+2), ifelse(c(fix.nugget.rel,fix.lambda),"fixed", "estimated"))
  if(cov.model == "pure.nugget"){
    par.su$status[beta.size+2] <- ""
  }
  par.su$values <- round(c(results$beta, results$cov.pars, results$nugget.rel, results$lambda), dig=4)
  if(beta.size == 1) beta.name <- "beta"
  else beta.name <- paste("beta", 0:(beta.size-1), sep="")
  row.names(par.su) <- c(beta.name, "sigmasq", "phi", "tausq.rel", "lambda")
  results$parameters.summary <- par.su
  class(results) <- "likGLSM"
  return(results)
}

"print.likGLSM" <-
  function(x, digits = max(3, getOption("digits") - 3), ...)
{
  est.pars <- as.vector(x$parameters.summary[x$parameters.summary[,1] == "estimated",2])
  names.est.pars <- dimnames(x$parameters.summary[x$parameters.summary[,1] == "estimated",])[[1]]
  names(est.pars) <- names.est.pars
  cat("likfit.glsm: estimated model parameters:\n")
  print.default(format(est.pars, digits=digits), ...)
  cat("\n likfit.glsm : maximised log-likelihood = ")
  cat(format(x$loglik, digits=digits))
  cat("\n")
  return(invisible())
}

"summary.likGLSM" <-
  function(object, ...)
{
  names.pars <- dimnames(object$parameters.summary)[[1]]
  summ.lik <- list()
  summ.lik$method.lik <- "maximum likelihood"
  summ.lik$family <- object$family
  summ.lik$link <- object$link
  summ.lik$mean.component <- object$beta
  names(summ.lik$mean.component) <- names.pars[seq(along=object$beta)]
  summ.lik$cov.model <- object$cov.model
  summ.lik$kappa <- object$kappa
  summ.lik$aniso.pars <- object$aniso.pars
  summ.lik$spatial.component <- object$parameters.summary[c("sigmasq", "phi"),]
  summ.lik$nugget.component <- object$parameters.summary[c("tausq.rel"),, drop=FALSE]
  summ.lik$transformation  <- object$parameters.summary[c("lambda"),, drop=FALSE]
  summ.lik$likelihood <- list(log.L = object$loglik, n.params = as.integer(object$npars))
  summ.lik$estimated.pars <- dimnames(object$parameters.summary[object$parameters.summary[,1] == "estimated",])[[1]]
  summ.lik$call <- object$call
  class(summ.lik) <- "summary.likGLSM"
  return(summ.lik)
}

"print.summary.likGLSM" <-
  function(x, digits = max(3, getOption("digits") - 3), ...)
{
  if(length(class(x)) == 0 || all(class(x) != "summary.likGLSM"))
    stop("object is not of the class \"summary.likGLSM\"")
  cat("Summary of the maximum likelihood parameter estimation\n")
  cat("-----------------------------------\n")
  cat(paste("Family = ", x$family, ", Link = ", x$link, "\n" ))
  cat("\n")
  cat("Parameters of the mean component (trend):")
  cat("\n")
  print.default(format(x$mean.component, digits=digits), ...)
  cat("\n")
  ##
  cat("Parameters of the spatial component:")
  cat("\n")
  cat(paste("   correlation function:", x$cov.model))
  if(x$cov.model == "matern" | x$cov.model == "powered.exponential" |
     x$cov.model == "cauchy" | x$cov.model == "gneiting.matern"){
    cat(paste("\n          kappa = ", x$kappa))
    if(x$cov.model == "matern" & (round(x$kappa, digits=digits)  == 0.5)) cat(" (exponential)")
  }
  cat("\n")
  if(!is.null(x$aniso.pars)){
    cat(paste("\n (fixed) anisotropy parameters (angle, ratio) = (", x$aniso.pars, "( \n"))
  }
  cat(paste("\n      (estimated) variance parameter sigmasq (partial sill) = ", format(x$spatial.component[1,2], dig=digits)))
  if(x$cov.model != "pure.nugget") cat(paste("\n      (estimated) cor. fct. parameter phi (range parameter)  = ", format(x$spatial.component[2,2], dig=digits)))
  cat("\n")
  if(x$nugget.component[,1] == "estimated")
    cat(paste("\n (estimated) relative nugget = ", format(x$nugget.component[,2], dig=digits)))
  else
    cat(paste("\n (fixed) relative nugget =", x$nugget.component[,2]))
  cat("\n")
  cat("\n")
  if(x$family == "poisson" && x$link == "boxcox"){
    cat("\n")
    cat("Transformation parameter:")
    cat("\n")
    lambda <- x$transformation[,2]
    if(x$transformation[,1] == "estimated")
      cat(paste("      (estimated) Box-Cox parameter =", format(lambda, dig=digits)))
    else{
      cat(paste("      (fixed) Box-Cox parameter =", lambda))
      if(abs(lambda - 1) <  0.0001) cat(" (no transformation)")
      if(abs(lambda) < 0.0001) cat(" (log-transformation)")
    }
  }
  cat("\n")
  cat("\n")
  cat("Maximised Likelihood:")
  cat("\n")
  print(format(x$likelihood, digits=digits))
  cat("\n")
  cat("Call:")
  cat("\n")
  print(x$call)
  cat("\n")
  invisible(x)
}


"glsm.mcmc" <- function(geodata, coords = geodata$coords, data = geodata$data, units.m = "default", model, mcmc.input, messages)
{
  if(missing(geodata))
    geodata <- list(coords=coords, data=data, units.m=units.m)
  call.fc <- match.call()
  if(missing(messages))
    messages.screen <- ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages"))
  else messages.screen <- messages
  n <- length(data)
  if(any(units.m == "default")){
    if(!is.null(geodata$units.m)) units.m <- geodata$units.m
    else units.m <- rep(1, n)
  }
  ##
  if(is.vector(coords)){
    coords <- cbind(coords, 0)
    warning("vector of coordinates: one spatial dimension assumed")
  }
  coords <- as.matrix(coords)
  dimnames(coords) <- list(NULL, NULL)
  ##
  model <- model.glsm.mcmc.check.aux(model)
  cov.model <- model$cov.model
  kappa <- model$kappa
  beta <- model$beta
  cov.pars <- model$cov.pars
  nugget <- model$nugget
  aniso.pars <- model$aniso.pars
  trend <- model$trend
  family <- model$family
  link <- model$link
  lambda <- model$lambda
  ##
  trend.data <- unclass(trend.spatial(trend=trend, geodata = geodata))
  if(nrow(trend.data) != n) stop("length of trend is different from the length of the data")
  if(ncol(trend.data) != length(beta)) stop("size of beta is incompatible with trend specified") 
  ##
  ## preparing for MCMC 
  ##
  if(missing(mcmc.input)) stop("glsm.mcmc: argument mcmc.input must be given")
  mcmc.input <- mcmc.check.aux(mcmc.input, fct="glsm.mcmc")
  ##
  mean.d <- as.vector(trend.data %*% beta)
  if(!is.null(aniso.pars)) {
    invcov <- varcov.spatial(coords = coords.aniso(coords = coords, aniso.pars = aniso.pars), cov.model = cov.model, kappa = kappa, 
                             nugget = nugget, cov.pars = cov.pars, inv = TRUE, func.inv = "cholesky",
                             try.another.decomposition = FALSE)$inverse
  }
  else {
    invcov <- varcov.spatial(coords = coords, cov.model = cov.model, kappa = kappa, nugget = nugget, cov.pars = cov.pars,
                             inv = TRUE, func.inv = "cholesky", try.another.decomposition = FALSE)$inverse
  }
  ##
########################----- MCMC ------#####################
  ##
  if(model$family == "binomial"){
    simulations <- mcmc.binom.logit(data = data, units.m = units.m, meanS = mean.d, invcov=invcov, mcmc.input = mcmc.input, messages.screen=messages.screen)
  }
  else{
    if(lambda == 0){
      simulations <- mcmc.pois.log(data = data, units.m = units.m, meanS = mean.d, invcov=invcov, mcmc.input = mcmc.input, messages.screen=messages.screen)
    }
    else{
      simulations <- mcmc.pois.boxcox(data = data, units.m = units.m, meanS = mean.d, invcov=invcov, mcmc.input = mcmc.input, messages.screen=messages.screen, lambda=lambda)
    }
  }
  kpl.result <- list(simulations=simulations$Sdata, acc.rate=simulations$acc.rate, model=model, geodata=geodata)
  kpl.result$call <- call.fc
#######################################
  class(kpl.result) <- "glsm.mcmc"
  return(kpl.result)
}


"model.glsm.mcmc.check.aux" <-
  function(model)
{
  if(class(model)=="likGLSM"){
    model$nugget <- model$nugget.rel*model$cov.pars[1]
    return(model)
  }
  else{
    if(!is.list(model)) stop("glsm.mcmc : the argument model only takes a list ")
    model.names <- c("beta", "cov.pars", "trend", "cov.model", "kappa", "aniso.pars", "nugget", "nugget.rel", "family", "link", "lambda")
    model <- object.match.names(model,model.names)
    if(is.null(model$beta) | !is.numeric(model$beta)) stop("glsm.mcmc : need to provide beta in the model argument")
    if(is.null(model$cov.pars)) stop("glsm.mcmc : need to provide cov.pars in the model argument")
    if(is.null(model$trend)) model$trend <- "cte"   
    if(is.null(model$cov.model)) model$cov.model <- "matern"
    model$cov.model <- match.arg(model$cov.model,
                                 choices = c("matern", "exponential","gaussian",
                                   "spherical", "circular", "cubic",
                                   "wave", "power",
                                   "powered.exponential", "cauchy", "gneiting",
                                   "gneiting.matern", "pure.nugget"))
    if(model$cov.model == "power") stop("krige.glm.control: correlation function does not exist for the power variogram")
    if(is.null(model$kappa)) model$kappa <- 0.5
    if(!is.null(model$aniso.pars))
      if(length(model$aniso.pars) != 2 | !is.numeric(model$aniso.pars))
        stop("glsm.mcmc : anisotropy parameters must be provided as a numeric vector with two elements: the rotation angle (in radians) and the anisotropy ratio (a number greater than 1)")
    if(is.null(model$nugget)){
      if(!is.null(model$nugget.rel)) model$nugget <- model$nugget.rel**model$cov.pars[1]
      else model$nugget <- 0
    }
    if(is.null(model$family)) stop("glsm.mcmc : need to provide family in the model argument")
    family <- match.arg(model$family, choices = c("poisson", "binomial"))
    if(family=="poisson"){
      if(is.null(model$lambda)){
        if(is.null(model$link)){
          model$lambda <- 0
          model$link <- "log" 
        }
        if(model$link == "canonical") model$link <- "log"
        if(model$link == "boxcox") stop("glsm.mcmc : need to provide lambda in the model argument")
        if(model$link == "log") model$lambda <- 0
        if(model$link == "id") model$lambda <- 1
      }
      if(!is.null(model$lambda)){
        if(is.null(model$link)){
          if(model$lambda > 0) model$link <- "boxcox"
          if(model$lambda == 0) model$link <- "log"
        }
        if(!is.null(model$link)){
          if(model$link == "canonical"){
            if(model$lambda > 0) model$link <- "boxcox"
            if(model$lambda == 0) model$link <- "log"
          }
          if(model$link == "boxcox" & model$lambda == 0) model$link <- "log"
          if(model$link == "log" & model$lambda > 0) warning("glsm.mcmc : value of argument lambda will be ignored since it is inconsistent with log-link ")
          if(model$link == "id" & model$lambda < 1) warning("glsm.mcmc : value of argument lambda will be ignored since it is inconsistent with identity-link ")
        }
      }
    }
    if(family=="binomial"){
      if(!is.null(model$link)){
        if(model$link != "logit" & model$link != "canonical") stop("glsm.mcmc : only the canonical logit link function is implemented ")
      }
      model$link <- "logit"
      model$lambda <- NULL
    }
    return(model)
  }
}


"create.mcmc.coda" <- function(x, mcmc.input)
{
  require(coda)
  if(exists("mcmc.input")){
    if(!is.list(mcmc.input)) stop(" mcmc.input must be given as a list ")
    if(is.null(mcmc.input$thin)) thin <- 10
    else thin <- mcmc.input$thin
    if(is.null(mcmc.input$burn.in)) st.val <- 1
    else st.val <- mcmc.input$burn.in + 1
  }
  else {
    thin <- 10
    st.val <- 1
  }
  if(class(x)=="glsm.mcmc"){
    n <- nrow(x$simulations)
    S.names <- rep(NA,n)
    for(i in 1:n) S.names[i] <- paste("S[",i,"]",sep="")
    temp <- t(x$simulations)
    colnames(temp) <- S.names
    return(mcmc(data=temp, start = st.val, thin=thin))
  }
  else{
    if(class(x)=="glm.krige.bayes"){
      n <- nrow(x$posterior$simulations)
      S.names <- rep(NA,n)
      for(i in 1:n) S.names[i] <- paste("S[",i,"]",sep="")
      if(!is.null(x$prior$phi$status) && x$prior$phi$status =="fixed"){
        temp <- t(x$posterior$simulations)
        colnames(temp) <- S.names
      }
      else{
        temp <- cbind(x$posterior$phi$sample,t(x$posterior$simulations))
        colnames(temp) <- c("phi", S.names)
      }
      return(mcmc(data=temp, start = st.val, thin=thin))
    }
    else{
      if(is.matrix(x)) return(mcmc(data=t(x), start = st.val, thin=thin))
      if(is.vector(x)) return(mcmc(data=x, start = st.val,thin=thin))
    }
  }
}

 

"asympvar" <- 
  function(timeseries, type = "mon", lag.max = 100, messages)
{
  if(missing(messages))
    messages.screen <- ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages"))
  else messages.screen <- messages
  ##
  if(is.vector(timeseries)) n.series <- 1
  else n.series <- nrow(timeseries) 
  if(type == "mon" | type == "all" | type == "pos") {
    if(messages.screen & type == "mon")
      cat(paste("calculating the initial monotone sequence estimate \n"))
    if(messages.screen & type == "pos") 
      cat(paste("calculating the initial positive sequence estimate \n"))
    if(messages.screen & type == "all") 
      cat(paste("calculating the initial positive sequence estimate, and the initial monotone sequence estimate \n"))      
  }
  else stop("Must specify type as either: mon, pos or all")
  len.Gamma <- floor(lag.max/2)-1
  if(n.series == 1){
     asy.gamma <- acf(timeseries, type = "covariance", plot = FALSE, lag.max = lag.max)$acf
     asy.gamma1 <- c(asy.gamma[(1 + 2 * c(0:len.Gamma))])
     asy.gamma2 <- c(asy.gamma[(2 + 2 * c(0:len.Gamma))])
     asy.Gamma <- asy.gamma1 + asy.gamma2
     ##--------- initial monotone sequence estimate -----------#
     kmaxpos <- min(c(which(asy.Gamma<0)-1, len.Gamma))
     if(type == "all" | type =="mon"){
       kmax <- min(c(which(diff(asy.Gamma)>0),kmaxpos))
       monvarest <- 2*sum(asy.Gamma[seq(length=kmax)])-asy.gamma[1]
       if(kmax == len.Gamma) warning("value of argument lag.max is not suffiently long")
     }
     ##--------- initial positive sequence estimate -----------#
     if(type == "pos" | type == "all"){
       posvarest <- 2*sum(asy.Gamma[seq(length=kmaxpos)])-asy.gamma[1]
       if (kmaxpos == len.Gamma) warning("value of argument lag.max is not suffiently long")
     }   
  }
  else{
     if(type == "all" | type == "pos") posvarest <- rep(1,n.series)
     if(type == "all" | type == "mon") monvarest <- rep(1,n.series)
     for(i in seq(length=n.series)){     
        asy.gamma <- acf(timeseries[i,], type = "covariance", plot = FALSE, lag.max = lag.max)$acf
        asy.gamma1 <- c(asy.gamma[(1 + 2 * c(0:len.Gamma))])
        asy.gamma2 <- c(asy.gamma[(2 + 2 * c(0:len.Gamma))])
        asy.Gamma <- asy.gamma1 + asy.gamma2
        ##--------- initial monotone sequence estimate -----------#
        kmaxpos <- min(c(which(asy.Gamma<0)-1, len.Gamma))
        if(type == "all" | type =="mon"){
          kmax <- min(c(which(diff(asy.Gamma)>0),kmaxpos))
          monvarest[i] <- 2*sum(asy.Gamma[seq(length=kmax)])-asy.gamma[1]
          if(kmax == len.Gamma) warning("value of argument lag.max is not suffiently long")
        }
        ##--------- initial positive sequence estimate -----------#
        if(type == "pos" | type == "all"){
          posvarest[i] <- 2*sum(asy.Gamma[seq(length=kmaxpos)])-asy.gamma[1]
          if (kmaxpos == len.Gamma) warning("value of argument lag.max is not suffiently long")
        }
     }
  }
  if(type == "pos") return(posvarest)
  if(type == "all") return(list(posvarest = posvarest, monvarest = monvarest))
  if(type == "mon") return(monvarest)
}

#### consider vectorising the while loops.



"covariog" <-  function(geodata, coords = geodata$coords, data = geodata$data, units.m = "default", uvec = "default", bins.lim = "default",
                        estimator.type = c("poisson", "not-poisson"), max.dist = NULL, pairs.min = 2)
{
  call.fc <- match.call()
  estimator.type <- match.arg(estimator.type)
  coords <- as.matrix(coords)
  data <- as.matrix(data)
  n.data <- nrow(coords)
  if(n.data != nrow(data)) stop("Dimension of data and coords do not match ")
  n.datasets <- ncol(data)
  if(any(units.m == "default")){
    if(!is.null(geodata$units.m)) units.m <- geodata$units.m
    else units.m <- rep(1, n.data)
  }
  data1 <- as.matrix(data/units.m)
  if(n.datasets == 1)
    data1 <- as.vector(data1)
  u <- as.vector(dist(as.matrix(coords)))
  if(!is.null(max.dist))
    umax <- max(u[u < max.dist])
  else umax <- max(u)
  if(all(bins.lim == "default")) {
    if(!is.null(max.dist))
      umax <- max(u[u < max.dist])
    else umax <- max(u)
    if(all(uvec == "default"))
      uvec <- seq(0, umax, l = 15)
    nvec <- length(uvec)
    d <- 0.5 * diff(uvec[seq(length=nvec)])
    bins.lim <- c(0, (uvec[seq(length=(nvec - 1))] + d), (d[nvec - 1] + uvec[nvec]))
    if(uvec[1] == 0)
      uvec[1] <- (bins.lim[1] + bins.lim[2])/2
  }
  else {
    if(!all(uvec == "default")) warning(" Both uvec and bins.lim are specified; using values in bins.lim ")
    if(!is.null(max.dist))
      if(max.dist != max(bins.lim)) stop("conflict between max.dist and max(bins.lim)")
  }
  nbins <- length(bins.lim) - 1
  if(is.null(max.dist))
    max.dist <- max(bins.lim)
  temp.list <- list(n.data = n.data, coords = coords, data = data1, nbins = nbins, bins.lim = bins.lim, max.dist = max.dist)
  bin.f <- function(data, temp.list)
    {
      result <- .C("binitprod",
                   as.integer(temp.list$n.data),
                   as.double(as.vector(temp.list$coords[, 1])),
                   as.double(as.vector(temp.list$coords[, 2])),
                   as.double(as.vector(data)),
                   as.integer(temp.list$nbins),
                   as.double(as.vector(temp.list$bins.lim)),
                   as.double(temp.list$max.dist),
                   cbin = as.integer(rep(0, temp.list$nbins)),
                   vbin = as.double(rep(0, temp.list$nbins)), DUP=FALSE, PACKAGE = "geoRglm")[c("vbin", "cbin")]
      return(result)
    }
  result <- array(unlist(lapply(as.data.frame(data1), bin.f, temp.list = temp.list)), dim = c(nbins, 2, n.datasets))
  mm2 <- colMeans(as.matrix(data1))^2
  sigma2 <- log(t(t((apply(as.matrix(data1), 2, var) * (n.data - 1))/n.data - colMeans(as.matrix(data1/units.m)))/mm2) + 1)
  if(estimator.type == "poisson")
    calcresult <- array(rbind(as.vector(sigma2), log(t(t(as.matrix(result[, 1,  ]))/mm2))), dim = c(nbins + 1, n.datasets))
  else calcresult <- array(log(t(t(as.matrix(result[, 1,  ]))/mm2)), dim = c(nbins, n.datasets))
  indp <- (result[, 2, 1] >= pairs.min)
  if(estimator.type == "poisson")
    result <- list(u = c(0, uvec[indp]), v = calcresult[indp,  ], n = c(n.data, result[indp, 2, 1]), bins.lim = bins.lim)
  else result <- list(u = uvec[indp], v = calcresult[indp,  ], n = result[indp, 2, 1], bins.lim = bins.lim)
  result <- c(result, list(v0 = sigma2, estimator.type = estimator.type, n.data = n.data, call = call.fc))
  class(result) <- "covariogram"
  return(result)
}


"covariog.model.env" <- 
function(geodata, coords = geodata$coords, units.m = "default", obj.covariog, model.pars, nsim = 500, prob = c(0.025, 0.975), 
     messages)
{
  call.fc <- match.call()
  if(missing(messages))
    messages.screen <- ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages"))
  else messages.screen <- messages
  ##
  ## reading input
  ##
  if(any(units.m == "default")){
    if(!is.null(geodata$units.m)) units.m <- geodata$units.m
    else units.m <- rep(1, obj.covariog$n.data)
  }
  if(is.null(model.pars$beta) | !is.numeric(model.pars$beta)) {
    stop("argument beta must be provided in order to perform make covariogram envelopes")
  }
  else if(length(model.pars$beta) > 1) {
    stop("only one-dimensional beta is alllowed")
  }
  else beta <- model.pars$beta
  if(!is.null(model.pars$cov.model)) {
    cov.model <- model.pars$cov.model
  }
  else cov.model <- "exponential"
  if(!is.null(model.pars$kappa))
    kappa <- model.pars$kappa
  else kappa <- 0.5
  if(!is.null(model.pars$nugget))
    nugget <- model.pars$nugget
  else nugget <- 0
  cov.pars <- model.pars$cov.pars
  ##
  ## generating simulations from the model with parameters provided
  ##
  if(messages.screen) cat(paste("covariog.env: generating", nsim, "simulations(", obj.covariog$n.data, 
          "points). \n"))
  simula <- pois.log.grf(obj.covariog$n.data, grid = as.matrix(coords), units.m = units.m, beta = beta, cov.model = cov.model, 
          cov.pars = cov.pars, nugget = nugget, kappa = kappa, nsim = nsim, messages = FALSE)
  ##
  ## computing empirical covariograms for the simulations
  ##
  if(messages.screen) cat(paste("covariog.env: computing the empirical covariogram for the", nsim, "simulations\n"))
  simula.result <- covariog(simula, bins.lim = obj.covariog$bins.lim, units.m = units.m, estimator.type = obj.covariog$estimator.type,
            pairs.min = min(obj.covariog$n))
  ##
  ## computing envelopes
  ##
  if(messages.screen) cat("covariog.env: computing the envelopes\n")
  limits <- apply(simula.result$v, 1, quantile, prob = prob)
  if(length(prob) == 1) {
    res.env <- list(u = obj.covariog$u, v = limits)
  }
  else {
    res.env <- list(u = obj.covariog$u, v = t(limits))
  }
  res.env$call <- call.fc
  return(res.env)
}


"lines.covariomodel" <- 
function(x, max.dist = x$max.dist, ...)
{
  if (is.null(max.dist)) stop("argument max.dist needed for this object")
  my.l <- x
  if(is.null(my.l$nugget)) my.l$nugget <- 0
  C.f <- function(x, my.l){ 
    return(ifelse(x>0, 0, my.l$nugget)+cov.spatial(x, cov.model = my.l$cov.model, kappa = my.l$kappa, cov.pars = my.l$cov.pars))
  }
  curve(C.f(x,my.l=my.l), from = 0, to = max.dist, add=TRUE, ...)
  return(invisible())
}


"plot.covariogram" <- 
function(x, max.dist = max(x$u), ylim = "default", type = "b", envelope.obj = NULL, ...)
{
        u <- x$u[x$u <= max.dist]
	if(!is.matrix(x$v)) {
		v <- x$v[x$u <= max.dist]
		ymax <- max(x$v[x$u <= max.dist])
		ymin <- min(x$v[x$u <= max.dist & x$v >  - Inf])
	}
	else {
		v <- x$v[x$u <= max.dist, 1]
		ymax <- max(x$v[x$u <= max.dist,  ])
		temp <- x$v[x$u <= max.dist,  ]
		ymin <- min(temp[temp >  - Inf])
	}
	if(!is.null(envelope.obj)) {
		ymax <- max(envelope.obj$v, ymax)
		ymin <- min(envelope.obj$v[envelope.obj$v >  - Inf], ymin)
	}
	if(is.numeric(ylim))
		plot(u, v, xlim = c(0, max.dist), ylim = ylim, xlab = "Distance", ylab = "Covariance", type = type, ...)
	else plot(u, v, xlim = c(0, max.dist), ylim = c(ymin, ymax), xlab = "Distance", ylab = "Covariance", type = type, ...)
	abline(h = 0, lty = 3)
	if(is.matrix(x$v)) {
		for(k in seq(2,ncol(x$v))) {
			v <- x$v[x$u <= max.dist, k]
			lines(u, v, type = "b")
		}
	}
	if(!is.null(envelope.obj)) {
		if(ncol(as.matrix(envelope.obj$v)) == 1)
			lines(envelope.obj$u, envelope.obj$v, lty = 4)
		else apply(envelope.obj$v, 2, lines, x = envelope.obj$u, lty = 4)
	}
	return(invisible())
}


"pois.log.grf" <- 
function(n = NULL, grid, nx = round(sqrt(n)), ny = round(sqrt(n)), nsim = 1,
        xlims = c(0, 1), ylims = c(0, 1), units.m = "default", trend = "cte", beta = stop("beta parameter(s) needed"), 
        cov.model = c("exponential", "matern", "gaussian", "spherical", "cubic", "wave", "powered.exponential", "cauchy",
        "gneiting", "gneiting.matern", "pure.nugget"), cov.pars = stop("covariance parameters (sigmasq and phi) needed"), nugget = 0,
        kappa = 0.5, method = c("cholesky", "svd", "eigen", "circular.embedding"), messages, ...)
{
#
# function for simulation of a Poisson log-Gaussian random field.
#
# the syntax is similar to grf, but beta needs to be specified.
## beta   : vector of regression parameters
## trend  : by default a constant mean
## units.m   : n-vector of observation-times for data (default is 1 for all observations) 
#        
  if(missing(messages))
    messages.screen <- ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages"))
  else messages.screen <- messages
  ##
  call.fc <- match.call()
  method <- match.arg(method)
  cov.model <- match.arg(cov.model)
  work <- grf(n = n, grid = grid, nx = nx, ny = ny, nsim = nsim, xlims = xlims, ylims = ylims, cov.model = cov.model, cov.pars = 
              cov.pars, nugget = nugget, lambda = 1, kappa = kappa, method = method, messages =  messages.screen, ...)
  trend.data <- unclass(trend.spatial(trend = trend, geodata = work))
  n <- nrow(work$coords)
  if((length(beta) == 1) & trend == "cte")
    mu <- rep(beta, n)
  else mu <- trend.data %*% as.vector(beta)
  if(all(units.m == "default"))
    units.m <- rep(1, n)
  lambda <- matrix(rep(units.m * exp(mu), nsim), nrow = n, ncol = nsim) * exp(as.matrix(work$data))
  zsim <- matrix(rpois(n * nsim, lambda = lambda), nrow = n, ncol = nsim)
  if(nsim == 1)
    zsim <- as.vector(zsim)
  results <- list(coords = work$coords, data = zsim, cov.model = cov.model, nugget = nugget, cov.pars = cov.pars, kappa = kappa,
                  mu = mu, units.m = units.m, messages = work$messages, method = method)
  results$.Random.seed <- work$.Random.seed
  results$call <- call.fc
  class(results) <- "geodata"
  return(results)
}





"mcmc.bayes.pois.log" <- 
  function(data, units.m, trend, mcmc.input, messages.screen, cov.model, kappa, tausq.rel, coords, ss.sigma, df, phi.prior, phi.discrete)
{
  ##
#### This is the MCMC engine for the Bayesian analysis of a spatial Poisson log Normal model
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  if(any(mcmc.input$Htrunc=="default")) Htrunc <- 2*data + 5 
  else {
    if(is.vector(mcmc.input$Htrunc) & length(mcmc.input$Htrunc) == n) Htrunc <- mcmc.input$Htrunc
    else Htrunc <- rep(mcmc.input$Htrunc, n)
  }
  if(any(mcmc.input$S.start=="default")) {
    S <- as.vector(ifelse(data > 0, log(data), 0) - log(units.m))
  }
  else{
    if(!any(mcmc.input$S.start=="random")){
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        S <- as.vector(mcmc.input$S.start)
      }
      else  stop(" S.start must be a vector of same dimension as data ")
    }
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
  if(any(mcmc.input$phi.start=="default")) phi <- median(phi.discrete)
  else  phi <- mcmc.input$phi.start
  nmphi <-  length(phi.discrete)
  if(is.null(mcmc.input$phi.scale)) {
    if(nmphi > 1) stop("mcmc.input$phi.scale not given ")
    else phi.scale <- 0
  }
  else {
    phi.scale <- mcmc.input$phi.scale
    if(nmphi > 1 && pnorm((phi.discrete[nmphi] - phi.discrete[1])/(nmphi - 1), sd = sqrt(phi.scale)) > 0.975)
      warning("Consider making the grid in phi.discrete more dense. The algorithm may have problems moving.")
  }
  messages.C <- ifelse(messages.screen,1,0)
  ##                                                                      
##### ---------- sampling ----------- ###### 
  cov.model.number <- cor.number(cov.model)
  if(is.vector(trend)) beta.size <- 1
  else beta.size <- ncol(trend)
  n.sim <- floor(n.iter/thin)
  ## remeber this rather odd coding for telling that S.start is from the prior !!!
  if(any(mcmc.input$S.start=="random")) Sdata <- as.double(as.vector(c(rep(0, n.sim*n - 1),1)))
  else Sdata <- as.double(as.vector(c(S, rep(0, (n.sim - 1) * n))))
  result <-  .C("mcmcrun4",
                as.integer(n),
                as.double(data),
                as.double(units.m),
                as.double(as.vector(t(trend))),
                as.integer(beta.size),
                as.integer(cov.model.number),
                as.double(kappa),
                as.double(tausq.rel),
		as.double(coords[,1]),
                as.double(coords[,2]),
                as.double(S.scale),
                as.double(phi.scale),
                as.double(Htrunc),
                as.integer(n.iter),
                as.integer(thin),
                as.integer(burn.in),
                as.integer(messages.C),
                as.double(ss.sigma),
                as.integer(df),
                as.double(phi.prior),
                as.double(phi.discrete),
                as.integer(nmphi),
                Sdata = Sdata,
                phi.sample = as.double(rep(phi, n.sim)),
		acc.rate = rep(0,floor(n.iter/1000)+1), 
		acc.rate.phi = rep(0,floor(n.iter/1000)+1), DUP=FALSE, PACKAGE = "geoRglm")[c("Sdata", "phi.sample","acc.rate","acc.rate.phi" )]
  attr(result$Sdata, "dim") <- c(n, n.sim)
  if(nmphi>1) result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate,result$acc.rate.phi))
  else result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate))
  result$acc.rate.phi <- NULL
  if(burn.in==0) result$acc.rate <- result$acc.rate[-1,]
  if(nmphi>1) names(result$acc.rate) <- c("iter.numb", "Acc.rate", "Acc.rate.phi")
  else names(result$acc.rate) <- c("iter.numb", "Acc.rate")
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n"))
  return(result)
}


"mcmc.bayes.pois.boxcox" <- 
  function(data, units.m, trend, mcmc.input, messages.screen, cov.model, kappa, tausq.rel, coords, ss.sigma, df, phi.prior, phi.discrete, lambda) 
{
  ##
#### This is the MCMC engine for the Bayesian analysis of a spatial Poisson boxcox Normal model, when lambda >0
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  if(any(mcmc.input$S.start=="default")) {
    S <- as.vector(ifelse(data > 0, ((data/units.m)^lambda -1)/lambda, 0) )         
  }
  else{
    if(!any(mcmc.input$S.start=="random")){
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        S <- as.vector(mcmc.input$S.start)
      }
      else  stop(" S.start must be a vector of same dimension as data ")
    }
  }
  if(any(mcmc.input$Htrunc=="default")) Htrunc <- 2*data + 5 
  else {
    if(is.vector(mcmc.input$Htrunc) & length(mcmc.input$Htrunc) == n) Htrunc <- mcmc.input$Htrunc
    else Htrunc <- rep(mcmc.input$Htrunc, n)
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
  if(any(mcmc.input$phi.start=="default")) phi <- median(phi.discrete)
  else  phi <- mcmc.input$phi.start
  nmphi <-  length(phi.discrete)
  if(is.null(mcmc.input$phi.scale)) {
    if(nmphi > 1) stop("mcmc.input$phi.scale not given ")
    else phi.scale <- 0
  }
  else {
    phi.scale <- mcmc.input$phi.scale
    if(nmphi > 1 && pnorm((phi.discrete[nmphi] - phi.discrete[1])/(nmphi - 1), sd = sqrt(phi.scale)) > 0.975)
      warning("Consider making the grid in phi.discrete more dense. The algorithm may have problems moving. ")
  }
  messages.C <- ifelse(messages.screen,1,0)
  ##                                                                      
##### ---------- sampling ----------- ###### 
  cov.model.number <- cor.number(cov.model)
  if(is.vector(trend)) beta.size <- 1
  else beta.size <- ncol(trend)
  n.sim <- floor(n.iter/thin)
  ## remeber this rather odd coding for telling that S.start is from the prior !!!
  if(any(mcmc.input$S.start=="random")) Sdata <- as.double(as.vector(c(rep(0, n.sim*n - 1),1)))
  else Sdata <- as.double(as.vector(c(S, rep(0, (n.sim - 1) * n))))
  result <-  .C("mcmcrun4boxcox",
                as.integer(n),
                as.double(data),
                as.double(units.m),              
                as.double(as.vector(t(trend))),
                as.integer(beta.size),
                as.integer(cov.model.number),
                as.double(kappa),
                as.double(tausq.rel),
		as.double(coords[,1]),
                as.double(coords[,2]),
                as.double(S.scale),
                as.double(phi.scale),
                as.double(Htrunc),
                as.integer(n.iter),
                as.integer(thin),
                as.integer(burn.in),
                as.integer(messages.C),
                as.double(ss.sigma),
                as.integer(df),
                as.double(phi.prior),
                as.double(phi.discrete),
                as.integer(nmphi),
                as.double(lambda),
                Sdata = Sdata,
                phi.sample = as.double(rep(phi, n.sim)), 
		acc.rate = rep(0,floor(n.iter/1000)+1), 
		acc.rate.phi = rep(0,floor(n.iter/1000)+1), DUP=FALSE, PACKAGE = "geoRglm")[c("Sdata", "phi.sample","acc.rate","acc.rate.phi" )]   
  attr(result$Sdata, "dim") <- c(n, n.sim)
  if(nmphi>1) result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate,result$acc.rate.phi))
  else result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate))
  result$acc.rate.phi <- NULL
  if(burn.in==0) result$acc.rate <- result$acc.rate[-1,]
  if(nmphi>1) names(result$acc.rate) <- c("iter.numb", "Acc.rate", "Acc.rate.phi")
  else names(result$acc.rate) <- c("iter.numb", "Acc.rate")
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n"))
  return(result)
}


"mcmc.bayes.conj.pois.log" <- 
  function(data, units.m, meanS, ttvbetatt, mcmc.input, messages.screen, cov.model, kappa, tausq.rel, coords, ss.sigma, df, phi.prior, phi.discrete)
{
  ##
  ## This is the MCMC engine for the Bayesian analysis (with normal prior for beta) of a spatial Poisson logit Normal model
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  if(any(mcmc.input$S.start == "default")) {
    S <- as.vector(ifelse(data > 0, log(data), 0) - log(units.m)) - meanS
  }
  else{
    if(!any(mcmc.input$S.start=="random")){
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        S <- as.vector(mcmc.input$S.start)
      }
      else  stop(" S.start must be a vector of same dimension as data ")
    }
  }
  if(any(mcmc.input$Htrunc=="default")) Htrunc <- 2*data + 5 
  else {
    if(is.vector(mcmc.input$Htrunc) & length(mcmc.input$Htrunc) == n) Htrunc <- mcmc.input$Htrunc
    else Htrunc <- rep(mcmc.input$Htrunc, n)
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
  if(any(mcmc.input$phi.start=="default")) phi <- median(phi.discrete)
  else  phi <- mcmc.input$phi.start
  nmphi <-  length(phi.discrete)
  if(is.null(mcmc.input$phi.scale)) {
    if(nmphi > 1) stop("mcmc.input$phi.scale not given ")
    else phi.scale <- 0
  }
  else {
    phi.scale <- mcmc.input$phi.scale
    if(nmphi > 1 && pnorm((phi.discrete[nmphi] - phi.discrete[1])/(nmphi - 1), sd = sqrt(phi.scale)) > 0.975)
      warning("Consider making the grid in phi.discrete more dense. The algorithm may have problems moving. ")
  }
  messages.C <- ifelse(messages.screen,1,0)
  ##                                                                      
  ## ---------- sampling ----------- ###### 
  cov.model.number <- cor.number(cov.model)
  if(is.null(ttvbetatt)) ttvbetatt <- matrix(0,beta.size,beta.size)
  n.sim <- floor(n.iter/thin)
  ## remeber this rather odd coding for telling that S.start is from the prior !!!
  if(any(mcmc.input$S.start=="random")) Sdata <- as.double(as.vector(c(rep(0, n.sim*n - 1),1)))
  else Sdata <- as.double(as.vector(c(S, rep(0, (n.sim - 1) * n))))
  result <-  .C("mcmcrun5",
                as.integer(n),
                as.double(data),
                as.double(units.m),
                as.double(as.vector(meanS)),
                as.double(as.vector(ttvbetatt)),
                as.integer(cov.model.number),
                as.double(kappa),
                as.double(tausq.rel),
                as.double(coords[,1]),
                as.double(coords[,2]),
                as.double(S.scale),
                as.double(phi.scale),
                as.double(Htrunc),
                as.integer(n.iter),
                as.integer(thin),
                as.integer(burn.in),
                as.integer(messages.C),
                as.double(ss.sigma),
                as.integer(df),
                as.double(phi.prior),
                as.double(phi.discrete),
                as.integer(nmphi),
                Sdata = Sdata,
                phi.sample = as.double(rep(phi, n.sim)),
		acc.rate = rep(0,floor(n.iter/1000)+1), 
		acc.rate.phi = rep(0,floor(n.iter/1000)+1), DUP=FALSE, PACKAGE = "geoRglm")[c("Sdata", "phi.sample","acc.rate","acc.rate.phi" )]
  attr(result$Sdata, "dim") <- c(n, n.sim) 
  if(nmphi>1) result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate,result$acc.rate.phi))
  else result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate))
  result$acc.rate.phi <- NULL
  if(burn.in==0) result$acc.rate <- result$acc.rate[-1,]
  if(nmphi>1) names(result$acc.rate) <- c("iter.numb", "Acc.rate", "Acc.rate.phi")
  else names(result$acc.rate) <- c("iter.numb", "Acc.rate")
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n"))
  return(result)
}

"mcmc.bayes.conj.pois.boxcox" <- 
  function(data, units.m, meanS, ttvbetatt, mcmc.input, messages.screen, cov.model, kappa, tausq.rel, coords, ss.sigma, df, phi.prior,
           phi.discrete, lambda)
{
  ##
  ## This is the MCMC engine for the Bayesian analysis (with normal prior for beta) of a spatial Poisson logit Normal model
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  if(any(mcmc.input$S.start == "default")) {
    S <- as.vector(ifelse(data > 0, ((data/units.m)^lambda -1)/lambda, 0) ) - meanS
  }
  else{
    if(!any(mcmc.input$S.start=="random")){
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        S <- as.vector(mcmc.input$S.start)
      }
      else  stop(" S.start must be a vector of same dimension as data ")
    }
  }
  if(any(mcmc.input$Htrunc=="default")) Htrunc <- 2*data + 5 
  else {
    if(is.vector(mcmc.input$Htrunc) & length(mcmc.input$Htrunc) == n) Htrunc <- mcmc.input$Htrunc
    else Htrunc <- rep(mcmc.input$Htrunc, n)
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
  if(any(mcmc.input$phi.start=="default")) phi <- median(phi.discrete)
  else  phi <- mcmc.input$phi.start
  nmphi <-  length(phi.discrete)
  if(is.null(mcmc.input$phi.scale)) {
    if(nmphi > 1) stop("mcmc.input$phi.scale not given ")
    else phi.scale <- 0
  }
  else {
    phi.scale <- mcmc.input$phi.scale
    if(nmphi > 1 && pnorm((phi.discrete[nmphi] - phi.discrete[1])/(nmphi - 1), sd = sqrt(phi.scale)) > 0.975)
      warning("Consider making the grid in phi.discrete more dense. The algorithm may have problems moving. ")
  }
  messages.C <- ifelse(messages.screen,1,0)
  ##                                         
  ## ---------- sampling ----------- ###### 
  cov.model.number <- cor.number(cov.model)
  if(is.null(ttvbetatt)) ttvbetatt <- matrix(0,beta.size,beta.size)
  n.sim <- floor(n.iter/thin)
  ## remeber this rather odd coding for telling that S.start is from the prior !!!
  if(any(mcmc.input$S.start=="random")) Sdata <- as.double(as.vector(c(rep(0, n.sim*n - 1),1)))
  else Sdata <- as.double(as.vector(c(S, rep(0, (n.sim - 1) * n))))
  result <-  .C("mcmcrun5boxcox",
                as.integer(n),
                as.double(data),
                as.double(units.m),
                as.double(as.vector(meanS)),
                as.double(as.vector(ttvbetatt)),
                as.integer(cov.model.number),
                as.double(kappa),
                as.double(tausq.rel),
		as.double(coords[,1]),
                as.double(coords[,2]),
                as.double(S.scale),
                as.double(phi.scale),
                as.double(Htrunc),
                as.integer(n.iter),
                as.integer(thin),
                as.integer(burn.in),
                as.integer(messages.C),
                as.double(ss.sigma),
                as.integer(df),
                as.double(phi.prior),
                as.double(phi.discrete),
                as.integer(nmphi),
                as.double(lambda),
                Sdata = Sdata,
                phi.sample = as.double(rep(phi, n.sim)), 
		acc.rate = rep(0,floor(n.iter/1000)+1), 
		acc.rate.phi = rep(0,floor(n.iter/1000)+1), DUP=FALSE, PACKAGE = "geoRglm")[c("Sdata", "phi.sample","acc.rate","acc.rate.phi" )]
  attr(result$Sdata, "dim") <- c(n, n.sim)
  if(nmphi>1) result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate,result$acc.rate.phi))
  else result$acc.rate <- as.data.frame(cbind(burn.in + seq(0,floor(n.iter/1000))*1000, result$acc.rate))
  result$acc.rate.phi <- NULL
  if(burn.in==0) result$acc.rate <- result$acc.rate[-1,]
  if(nmphi>1) names(result$acc.rate) <- c("iter.numb", "Acc.rate", "Acc.rate.phi")
  else names(result$acc.rate) <- c("iter.numb", "Acc.rate")
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n"))
  return(result)
}

"pred.aux" <- function(S, coords, locations, model, prior, output, phi.posterior, link)
{
  n.sim <- ncol(S)
  ni <- nrow(locations)
  do.prediction <- ifelse(all(locations == "no"), FALSE, TRUE)
  beta.size <- ncol(unclass(trend.spatial(trend=model$trend.d, geodata = list(coords=coords))))
  lambda <- model$lambda
  ##
  temp.post <- list()
  temp.post$beta.mean <- array(NA, dim = c(beta.size, n.sim))
  temp.post$beta.var <- array(NA, dim = c(beta.size, beta.size, n.sim))
  temp.post$S2 <- rep(0, n.sim)
  if(do.prediction) {
    temp.pred <- list()
    temp.pred$mean <- array(NA, dim = c(ni, n.sim))
    temp.pred$var <- array(NA, dim = c(ni, n.sim))
    if(output$sim.predict) {
      num.pred <- 1
      pred.simulations <- array(NA, dim = c(ni, n.sim))
    }
    else {
      num.pred <- 0
      pred.simulations <- " no simulations from the predictive distribution "
    }
  }
  else num.pred <- 0
  model.temp <- model
  model.temp$lambda <- 1
  output.temp <- list(n.posterior = 0, n.predictive = num.pred, messages.screen = FALSE)
  prior.temp <- prior
  prior.temp$phi.prior <- "fixed"
  prior.temp$phi.discrete <- NULL
  prior.temp$tausq.rel.prior <- "fixed"
  prior.temp$priors.info <- NULL
  if(phi.posterior$phi.prior == "fixed" || length(phi.posterior$phi.discrete) == 1) {
    if(phi.posterior$phi.prior == "fixed") prior.temp$phi <- phi.posterior$phi
    else prior.temp$phi <- phi.posterior$phi.discrete
    temp.result <- krige.bayes.extnd(data = S, coords = coords, locations = locations,
                                     model = model.temp, prior = prior.temp, output = output.temp)
    temp.post$beta.mean <- temp.result$posterior$beta$pars$mean
    temp.post$beta.var <- temp.result$posterior$beta$pars$var
    temp.post$S2 <- temp.result$posterior$sigmasq$pars$S2
    if(do.prediction) {
      temp.pred$mean <- temp.result$predictive$mean
      temp.pred$var <- temp.result$predictive$variance
      if(output$sim.predict){
        if(link=="logit") pred.simulations <- plogis(temp.result$predictive$simulations)
        else{
          if(lambda==0) pred.simulations <- exp(temp.result$predictive$simulations)
          else pred.simulations <- BC.inv(temp.result$predictive$simulations, lambda)
        }
      }
    }
  }
  else {
    phi.discrete <- phi.posterior$phi.discrete
    len.phi.discrete <- length(phi.discrete)      
    step.phi.discrete <- phi.discrete[2] - phi.discrete[1]
    phi.table <- rep(0,len.phi.discrete)
    for(i in seq(length=len.phi.discrete)){
      phi.table[i] <- sum(ifelse(abs(phi.posterior$sample-phi.discrete[i])<0.5*step.phi.discrete,1,0))
    }
    phi.sample.unique <- phi.discrete[phi.table>0]
    phi.table <- phi.table[phi.table>0]
    len.phi.un <- length(phi.sample.unique)
    indic.phi <- array(rep(0, len.phi.un * max(phi.table)), dim = c(len.phi.un, max(phi.table)))
    for(i in seq(length=len.phi.un)){
      temp.num <- which(abs(phi.posterior$sample-phi.sample.unique[i])<0.5*step.phi.discrete)
      indic.phi[i, seq(along=temp.num)] <- temp.num
    }
    for(i in seq(length=len.phi.un)){
      id.phi.i <- indic.phi[i, seq(length=phi.table[i])]
      prior.temp$phi <- phi.sample.unique[i]
      if(phi.table[i]==1)
        temp.result <- krige.bayes(data = S[, id.phi.i], coords = coords, locations = locations, model = model.temp, prior = prior.temp, output = output.temp)
      else temp.result <- krige.bayes.extnd(data = S[, id.phi.i], coords = coords, locations = locations, 
                                            model = model.temp, prior = prior.temp, output = output.temp)
      temp.post$beta.mean[, id.phi.i] <- temp.result$posterior$beta$pars$mean         
      temp.post$beta.var[,  , id.phi.i] <- temp.result$posterior$beta$pars$var
      temp.post$S2[id.phi.i] <- temp.result$posterior$sigmasq$pars$S2       
      if(do.prediction) {
        temp.pred$mean[, id.phi.i] <- temp.result$predictive$mean
        temp.pred$var[, id.phi.i] <- temp.result$predictive$variance
        if(output$sim.predict) {
          if(link=="logit") pred.simulations[ , id.phi.i] <- plogis(temp.result$predictive$simulations)
          else{
            if(lambda==0) pred.simulations[ , id.phi.i] <- exp(temp.result$predictive$simulations)
            else pred.simulations[ , id.phi.i] <- BC.inv(temp.result$predictive$simulations, lambda)
          }
        }
      }
    }
  }
  remove("temp.result")
  if(do.prediction) return(list(temp.post=temp.post,temp.pred=temp.pred,pred.simulations=pred.simulations))
  else return(list(temp.post=temp.post))
}


"pred.quan.aux" <- 
  function(temp.pred, loc.coincide, df.model, ni, quantile.estimator)
{
  temp.med <- apply(temp.pred$mean, 1, median)
  temp.unc <- sqrt(apply(temp.pred$mean, 1, var) + apply(temp.pred$var, 1, median))
  not.accurate <- (!loc.coincide)
  diffe <- pmixed(temp.med, temp.pred,df.model)-0.5
  temp.med.new <- temp.med[not.accurate]+0.1*(temp.med[not.accurate]+0.1) # to get started
  inv.sl <- rep(0,ni)
  parms.temp <- list()
  while(any(not.accurate)){
    parms.temp$mean<-temp.pred$mean[not.accurate,,drop=FALSE]
    parms.temp$var<-temp.pred$var[not.accurate,,drop=FALSE]
    diffe.new <- pmixed(temp.med.new, parms.temp,df.model)-0.5
    inv.sl[not.accurate] <- (temp.med.new-temp.med[not.accurate])/(diffe.new-diffe[not.accurate])
    temp.med[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), temp.med.new,temp.med[not.accurate])
    diffe[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), diffe.new, diffe[not.accurate])
    not.accurate[not.accurate] <- ifelse(abs(diffe[not.accurate])>0.0005, TRUE, FALSE)
    temp.med.new <- temp.med[not.accurate] - diffe[not.accurate]*inv.sl[not.accurate]
  }
  temp.upper <- qnorm(rep(0.975, ni), mean = temp.med, sd = temp.unc)
  not.accurate <- (!loc.coincide)
  diffe <- pmixed(temp.upper, temp.pred,df.model)-0.975
  temp.upper.new <- temp.upper[not.accurate]+0.5*(temp.upper[not.accurate]+0.5) # to get started
  inv.sl <- rep(0,ni)      
  while(any(not.accurate)){
    parms.temp$mean<-temp.pred$mean[not.accurate,,drop=FALSE]
    parms.temp$var<-temp.pred$var[not.accurate,,drop=FALSE]
    diffe.new <- pmixed(temp.upper.new, parms.temp,df.model)-0.975
    inv.sl[not.accurate] <- (temp.upper.new-temp.upper[not.accurate])/(diffe.new-diffe[not.accurate])
    temp.upper[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), temp.upper.new,temp.upper[not.accurate])
    diffe[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), diffe.new, diffe[not.accurate])
    not.accurate[not.accurate] <- ifelse(abs(diffe[not.accurate])>0.0005, TRUE, FALSE)
    temp.upper.new <- temp.upper[not.accurate] - diffe[not.accurate]*inv.sl[not.accurate]
  }      
  temp.lower <- qnorm(rep(0.025, ni), mean = temp.med, sd = temp.unc)
  not.accurate <- (!loc.coincide)
  diffe <- pmixed(temp.lower, temp.pred,df.model)-0.025
  temp.lower.new <- temp.lower[not.accurate]+0.5*(temp.lower[not.accurate]+0.5) # to get started
  inv.sl <- rep(0,ni)
  while(any(not.accurate)){
    parms.temp$mean<-temp.pred$mean[not.accurate,,drop=FALSE]
    parms.temp$var<-temp.pred$var[not.accurate,,drop=FALSE]
    diffe.new <- pmixed(temp.lower.new,parms.temp,df.model)-0.025
    inv.sl[not.accurate] <- (temp.lower.new-temp.lower[not.accurate])/(diffe.new-diffe[not.accurate])
    temp.lower[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), temp.lower.new,temp.lower[not.accurate])
    diffe[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), diffe.new, diffe[not.accurate])
    not.accurate[not.accurate] <- ifelse(abs(diffe[not.accurate])>0.0005, TRUE, FALSE)
    temp.lower.new <- temp.lower[not.accurate] - diffe[not.accurate]*inv.sl[not.accurate]
  }
  if(any(loc.coincide)){
    temp.med[loc.coincide] <- apply(temp.pred$mean[loc.coincide,,drop=FALSE], 1, median)
    temp.upper[loc.coincide] <- apply(temp.pred$mean[loc.coincide,,drop=FALSE], 1, quantile, probs = 0.975)
    temp.lower[loc.coincide] <- apply(temp.pred$mean[loc.coincide,,drop=FALSE], 1, quantile, probs = 0.025) 
  }
  ## calculating quantiles
  if(is.logical(quantile.estimator) && (quantile.estimator)){
    temp.quan <- as.data.frame(cbind(temp.lower, temp.med, temp.upper))
  }
  if(is.numeric(quantile.estimator)){
    nmq <- length(quantile.estimator)
    if(nmq > 1) {
      temp.quan <- matrix(NA, ni, nmq)
      dig <- rep(3, nmq)
      for(i in seq(length=nmq)) {
        while(quantile.estimator[i] != round(quantile.estimator[i], digits = dig[i])) dig[i] <-dig[i] + 1
        temp.quan[, i] <- qnorm(rep(quantile.estimator[i], ni), mean = temp.med, sd = temp.unc)
        if(any(loc.coincide)) temp.quan[loc.coincide, i] <- temp.med[loc.coincide]
        not.accurate <- (!loc.coincide)
        diffe <- pmixed(temp.quan[,i], temp.pred,df.model)-quantile.estimator[i]
        numb <- 0.1+abs(quantile.estimator[i]-0.5)
        temp.quan.new <- temp.quan[not.accurate,i]+numb*(temp.quan[not.accurate,i]+numb) # to get started
        inv.sl <- rep(0,ni)
        while(any(not.accurate)) {
          parms.temp$mean <-temp.pred$mean[not.accurate,,drop=FALSE]
          parms.temp$var <-temp.pred$var[not.accurate,,drop=FALSE]
          diffe.new <- pmixed(temp.quan.new,parms.temp,df.model)-quantile.estimator[i]
          inv.sl[not.accurate] <- (temp.quan.new-temp.quan[not.accurate, i])/(diffe.new-diffe[not.accurate])
          temp.quan[not.accurate, i] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), temp.quan.new,temp.quan[not.accurate, i])
          diffe[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), diffe.new, diffe[not.accurate])
          not.accurate[not.accurate] <- ifelse(abs(diffe[not.accurate])>0.0005, TRUE, FALSE)
          temp.quan.new <- temp.quan[not.accurate, i] - diffe[not.accurate]*inv.sl[not.accurate]
        }            
        if(any(loc.coincide)){
          temp.quan[loc.coincide,i] <- apply(temp.pred$mean[loc.coincide,,drop=FALSE], 1, quantile, probs = quantile.estimator[i])
        }
      }
      temp.quan <- as.data.frame(temp.quan)
    }
    else {
      dig <- 3
      while(quantile.estimator != round(quantile.estimator,digits = dig)) dig <- dig + 1
      temp.quan <- qnorm(rep(quantile.estimator,ni), mean = temp.med, sd = temp.unc)
      not.accurate <- (!loc.coincide)
      diffe <- pmixed(temp.quan, temp.pred,df.model)-quantile.estimator
      numb <- 0.1+abs(quantile.estimator-0.5)
      temp.quan.new <- temp.quan[not.accurate]+numb*(temp.quan[not.accurate]+numb) # to get started
      inv.sl <- rep(0,ni)
      while(any(not.accurate)) {
        parms.temp$mean <-temp.pred$mean[not.accurate,,drop=FALSE]
        parms.temp$var <-temp.pred$var[not.accurate,,drop=FALSE]
        diffe.new <- pmixed(temp.quan.new,parms.temp,df.model)-quantile.estimator
        inv.sl[not.accurate] <- (temp.quan.new-temp.quan[not.accurate])/(diffe.new-diffe[not.accurate])
        temp.quan[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), temp.quan.new,temp.quan[not.accurate])
        diffe[not.accurate] <- ifelse(abs(diffe[not.accurate]) > abs(diffe.new), diffe.new, diffe[not.accurate])
        not.accurate[not.accurate] <- ifelse(abs(diffe[not.accurate])>0.0005, TRUE, FALSE)
        temp.quan.new <- temp.quan[not.accurate] - diffe[not.accurate]*inv.sl[not.accurate]
      }
      if(any(loc.coincide)){
        temp.quan[loc.coincide] <- apply(temp.pred$mean[loc.coincide,,drop=FALSE], 1, quantile, probs = quantile.estimator)
      }
      temp.quan <- as.vector(temp.quan)
    }
  }
  if(is.logical(quantile.estimator) && (quantile.estimator)){
    qname <- rep(0, 3)
    qname[1] <- paste("q0.025", sep = "")
    qname[2] <- paste("q0.5", sep = "")
    qname[3] <- paste("q0.975", sep = "")
    names(temp.quan) <- qname
  }
  if(is.numeric(quantile.estimator) && nmq > 1) {
    qname <- rep(0, length(quantile.estimator))
    for(i in seq(along=quantile.estimator))
      qname[i] <- paste("q", 100 * quantile.estimator[i], sep = "")
    names(temp.quan) <- qname
  }
  return(list(median=temp.med, upper=temp.upper, lower=temp.lower, quantiles = temp.quan))
}


"pois.krige.bayes" <- 
  function(geodata, coords = geodata$coords, data = geodata$data, units.m = "default", locations = "no", 
           model, prior, mcmc.input, output)
{
###########
  if(missing(geodata))
    geodata <- list(coords=coords, data=data, units.m=units.m)
  call.fc <- match.call()
  seed <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
  do.prediction <- ifelse(all(locations == "no"), FALSE, TRUE)
  ##
  ## Checking data configuration
  ##
  if(is.vector(coords)) {
    coords <- cbind(coords, 0)
    warning("vector of coordinates: one spatial dimension assumed")
  }
  coords <- as.matrix(coords)
  if(nrow(coords) != length(data))
    stop("number of data is different of number of data locations (coordinates)")
  n <- length(data)
  if(any(units.m == "default")){
    if(!is.null(geodata$units.m)) units.m <- geodata$units.m
    else units.m <- rep(1, n)
  }
  ####
  ## reading model input
  ##
  if(missing(model)) model <- model.glm.control()
  else model <- model.glm.check.aux(model, fct = "pois.krige.bayes")
  cov.model <- model$cov.model
  kappa <- model$kappa
  tausq.rel <- prior$tausq.rel
  lambda <- model$lambda
  if(lambda < 0) stop ("lambda < 0 is not allowed")
  ## reading prior input
  ##
  if(missing(prior)) stop("pois.krige.bayes: argument prior must be given")
  else prior <- prior.glm.check.aux(prior, fct = "pois.krige.bayes")
  beta.prior <- prior$beta.prior
  beta <- prior$beta
  beta.var <- prior$beta.var.std
  sigmasq.prior <- prior$sigmasq.prior
  if(sigmasq.prior == "fixed") sigmasq <- prior$sigmasq
  else{
    df.sigmasq <- prior$df.sigmasq
    S2.prior <- prior$sigmasq
  }
  phi.prior <- prior$phi.prior 
  phi <- prior$phi
  if(phi.prior != "fixed") phi.discrete <- prior$phi.discrete
  else phi.discrete <- phi
  ##
  ## reading output options
  ##
  if(missing(output)) output <- output.glm.control()
  else output <- output.glm.check.aux(output, fct = "pois.krige.bayes")
  quantile.estimator <- output$quantile.estimator
  probability.estimator <- output$probability.estimator  
  inference <- output$inference
  messages.screen <- output$messages.screen
  ## check == here
  data.dist <- as.vector(dist(coords))
  if(round(1000000000000. * min(data.dist)) == 0) stop("Two coords are identical; not allowed.")
  ##
  trend.d <- model$trend.d
  if(messages.screen) {
    cat(switch(as.character(trend.d)[1],
                 "cte" = "pois.krige.bayes: model with mean being constant",
                 "1st" = "pois.krige.bayes: model with mean given by a 1st order polynomial on the coordinates",
                 "2nd" = "pois.krige.bayes: model with mean given by a 2nd order polynomial on the coordinates",
                 "pois.krige.bayes: model with mean defined by covariates provided by the user"))
    cat("\n")
  }
  trend.data <- unclass(trend.spatial(trend=trend.d, geodata = geodata))
  dimnames(coords) <- list(NULL, NULL)
  dimnames(trend.data) <- list(NULL, NULL)
  beta.size <- ncol(trend.data)
  if(nrow(trend.data) != n) stop("length of trend is different from the length of the data")
  if(beta.size > 1)
    beta.names <- paste("beta", (0:(beta.size-1)), sep="")
  else beta.names <- "beta"
  if(beta.prior == "normal" |  beta.prior == "fixed"){
    if(beta.size != length(beta))
      stop("pois.krige.bayes: size of beta incompatible with the trend model (covariates)")
  }
  aniso.pars <- model$aniso.par
  if(!is.null(aniso.pars)) coords.transf <- coords.aniso(coords = coords, aniso.pars = aniso.pars)
  else coords.transf <- coords
  ##
  ## checking prediction locations
  ##
  if((inference) & (do.prediction)){
    ## Checking the consistency between coords, locations, and trends
    trend.l <- model$trend.l
    if(is.vector(locations)) {
      if(length(locations) == 2) {
        locations <- t(as.matrix(locations))
        warning("only one location to be predicted (in two-dimensional space) \n")
      }
      else locations <- as.matrix(cbind(locations, 0))
    }
    else locations <- as.matrix(locations)
    ni <- nrow(locations)
    ## Checking for 1D prediction 
    if(length(unique(locations[,1])) == 1 | length(unique(locations[,2])) == 1)
      krige1d <- TRUE
    else krige1d <- FALSE
    ##
    if(is.null(trend.l)) stop("trend.l needed for prediction")
    if(inherits(trend.d, "formula") | inherits(trend.l, "formula")){
      if((!inherits(trend.d, "formula")) | (!inherits(trend.l, "formula")))
        stop("trend.d and trend.l must have similar specification\n")
    }
    else{
      if((class(trend.d)=="trend.spatial") & (class(trend.l)=="trend.spatial")){
        if(ncol(trend.d) != ncol(trend.l))
          stop("trend.d and trend.l do not have the same number of columns")
      }
      else if(trend.d != trend.l) stop("trend.l is different from trend.d")
    }
    if(nrow(unclass(trend.spatial(trend=trend.l, geodata = list(coords = locations)))) != ni) 
      stop("pois.krige.bayes: number of points to be estimated is different of the number of trend locations")
    kb.results <- list(posterior = list(), predictive = list())
  }
  else {
    if(do.prediction & messages.screen) cat(paste("need to specify inference=TRUE to make predictions \n"))
    kb.results <- list(posterior = list(), predictive = paste("prediction not performed"))
    do.prediction <- FALSE
  }
  ##
  ## ##### preparing for MCMC -------------------------------------------------------
  ##
  if(missing(mcmc.input)) stop("pois.krige.bayes: argument mcmc.input must be given")
  mcmc.input <- mcmc.check.aux(mcmc.input, fct="pois.krige.bayes")
  ##
  if(beta.prior == "fixed" | beta.prior == "normal") mean.d <- as.vector(trend.data%*%beta)
  else mean.d <- rep(0,n)
  if(sigmasq.prior != "fixed"){
    if(beta.prior == "flat") df.model <- n - beta.size + df.sigmasq
    else df.model <- n + df.sigmasq
  }
  else df.model <- Inf
  if(beta.prior == "normal"){
    if(beta.size > 1) ttvbetatt <- trend.data%*%beta.var%*%t(trend.data)
    else ttvbetatt <- crossprod(t(trend.data))*beta.var
  }  
  else ttvbetatt <- 0
  if(sigmasq.prior == "fixed") {     ### implies that phi is fixed !
    invcov <- varcov.spatial(coords = coords, cov.model = cov.model, kappa = kappa, nugget = tausq.rel*sigmasq,
                               cov.pars = c(sigmasq,phi), inv = TRUE, func.inv = "cholesky",
                               try.another.decomposition = FALSE)$inverse
    if(beta.prior != "fixed"){
      ivtt <- invcov%*%trend.data
      if(beta.prior == "normal") invcov <- invcov-ivtt%*%solve.geoR(crossprod(trend.data, ivtt) + solve(beta.var), t(ivtt))
      else invcov <- invcov-ivtt%*%solve.geoR(crossprod(trend.data, ivtt), t(ivtt))
    }
  }
  if((phi.prior == "fixed") & (sigmasq.prior != "fixed")){
    phi.prior.prob <- 1
    phi.discrete <- phi
  }
  else phi.prior.prob <-  prior$priors.info$phi$probs   
  ##
############----------PART 2 ------------##############################
############-----------MCMC -------------##############################
  ##
  if(sigmasq.prior == "fixed"){ 
    if(lambda == 0){ 
      gauss.post <- mcmc.pois.log(data = data, units.m = units.m, meanS = mean.d, invcov=invcov, mcmc.input = mcmc.input, messages.screen=messages.screen)
    }
    else{
      gauss.post <- mcmc.pois.boxcox(data=data, units.m=units.m, meanS=mean.d, invcov=invcov, mcmc.input=mcmc.input, messages.screen=messages.screen, lambda=lambda)
    }
  } 
  else {
    kb.results$posterior$phi <- list()
    ## take care re-using gauss.post !
    if(beta.prior == "flat"){
      if(lambda == 0){ 
        gauss.post <- mcmc.bayes.pois.log(data=data, units.m=units.m, trend=trend.data, mcmc.input=mcmc.input, messages.screen=messages.screen, cov.model=cov.model, 
                                        kappa=kappa, tausq.rel = tausq.rel, coords=coords.transf, 
                                        ss.sigma = df.sigmasq*S2.prior, df = df.model, phi.prior = phi.prior.prob,
                                        phi.discrete = phi.discrete)
      }
      else{
        gauss.post <- mcmc.bayes.pois.boxcox(data=data, units.m=units.m, trend=trend.data, mcmc.input=mcmc.input, messages.screen=messages.screen, cov.model=cov.model, 
                                           kappa=kappa, tausq.rel = tausq.rel, coords=coords.transf,
                                           ss.sigma = df.sigmasq*S2.prior, df = df.model,  phi.prior = phi.prior.prob,
                                             phi.discrete = phi.discrete, lambda = lambda)
      }
    }
    else{
      if(lambda == 0){ 
        gauss.post <- mcmc.bayes.conj.pois.log(data=data, units.m=units.m, meanS = mean.d, ttvbetatt = ttvbetatt, mcmc.input=mcmc.input,
                                               messages.screen=messages.screen, cov.model=cov.model,  kappa=kappa, tausq.rel = tausq.rel,
                                               coords=coords.transf,  ss.sigma = df.sigmasq*S2.prior, df = df.model,
                                               phi.prior = phi.prior.prob, phi.discrete = phi.discrete,)
      }
      else{ 
        gauss.post <- mcmc.bayes.conj.pois.boxcox(data=data, units.m=units.m, meanS = mean.d, ttvbetatt = ttvbetatt,
                                                  mcmc.input=mcmc.input, messages.screen=messages.screen, cov.model=cov.model,  kappa=kappa, tausq.rel = tausq.rel,
                                                  coords=coords.transf,  ss.sigma = df.sigmasq*S2.prior, df = df.model,
                                                  phi.prior = phi.prior.prob, phi.discrete = phi.discrete, lambda = lambda)
      }
    }
    kb.results$posterior$phi$sample <- gauss.post$phi.sample
  }
  kb.results$posterior$acc.rate  <- gauss.post$acc.rate
  gauss.post <- gauss.post$Sdata
  ##           
##############-------------PART 3----------######################
##############------------prediction-------######################
  ## 
  n.sim <- ncol(gauss.post)
  if(inference){
    if(phi.prior=="fixed") phi.posterior <- list(phi.prior=phi.prior, phi=phi)
    else  phi.posterior <- list(phi.prior=phi.prior, phi.discrete=phi.discrete, sample=kb.results$posterior$phi$sample)
    predict.temp <- pred.aux(S=gauss.post, coords=coords, locations=locations, model=model, prior=prior, output=output, phi.posterior=phi.posterior, link="boxcox")
    temp.post <- predict.temp$temp.post
    if(do.prediction){
      temp.pred <- predict.temp$temp.pred
      kb.results$predictive$simulations <- predict.temp$pred.simulations
    }
    if(do.prediction) {
      ##
      d0mat <- loccoords(coords, locations)
      loc.coincide <- (colSums(d0mat < 1e-10) == 1)
      ##
      ## ------ median, quantiles and uncertainty 
      ##
      if((is.logical(quantile.estimator) && (quantile.estimator)) || (is.numeric(quantile.estimator))){
        predi.q <- pred.quan.aux(temp.pred, loc.coincide, df.model, ni, quantile.estimator)
        kb.results$predictive$median <- BC.inv(predi.q$median,lambda)
        kb.results$predictive$uncertainty <- (BC.inv(predi.q$upper,lambda) - BC.inv(predi.q$lower,lambda))/4      
        if(is.data.frame(predi.q$quantiles)){
          names.q <- names(predi.q$quantiles)
          kb.results$predictive$quantiles <- as.data.frame(BC.inv(as.matrix(predi.q$quantiles),lambda))
          names(kb.results$predictive$quantiles) <- names.q
        }
        else kb.results$predictive$quantiles <- BC.inv(predi.q$quantiles,lambda)
      }
      ##
      ## ------ probability estimators
      ##
      if(!is.null(probability.estimator)) {
        if(lambda == 0) transf.probab <- ifelse(probability.estimator > 0, log(probability.estimator), -1e+17)
          else transf.probab <- ifelse(probability.estimator > 0, (probability.estimator^lambda-1)/lambda, -1e+17)
        len.p <- length(probability.estimator)
        if(len.p==1){
          kb.results$predictive$probability <- round(pmixed(transf.probab, temp.pred, df.model), digits = 3)
        }
        else{
          kb.results$predictive$probability <- matrix(NA, ni,len.p)
          for(ii in seq(length=len.p)){
            kb.results$predictive$probability[,ii] <- round(pmixed(transf.probab[ii], temp.pred, df.model), digits = 3)
          }
        }
      }
      remove("temp.pred")
      ## 
      if(messages.screen) cat("pois.krige.bayes: Prediction performed \n")
    }
    else {
      kb.results$predictive <- "no locations to perform prediction were provided"
      if(messages.screen) cat(paste("Only Bayesian estimation of model parameters "))
    }
    ##
    ##----- calculating posterior summaries ----------------##
    ##
    if(beta.prior == "fixed") kb.results$posterior$beta <- paste("provided by user: ", beta) 
    else {
      kb.results$posterior$beta <- list()
      kb.results$posterior$beta$mean <- rowMeans(temp.post$beta.mean)
      names(kb.results$posterior$beta$mean) <- beta.names
      kb.results$posterior$beta$var <- rowMeans(temp.post$beta.var, dims=2) + var(t(temp.post$beta.mean))
      dimnames(kb.results$posterior$beta$var) <- list(beta.names,beta.names)
    }
    if(sigmasq.prior == "fixed") kb.results$posterior$sigmasq <- paste("provided by user: ", sigmasq) 
    else{
      kb.results$posterior$sigmasq <- list()
      kb.results$posterior$sigmasq$mean <- mean(temp.post$S2)*df.model/(df.model-2)
      kb.results$posterior$sigmasq$var <- (mean(temp.post$S2)*2/(df.model-4) + var(temp.post$S2))*df.model^2/(df.model-2)^2
    }
    if(phi.prior == "fixed") kb.results$posterior$phi <- paste("provided by user: ", phi) 
    else{
      kb.results$posterior$phi$mean <- mean(kb.results$posterior$phi$sample)
      kb.results$posterior$phi$var <- var(kb.results$posterior$phi$sample)
    }   
    ##
    ## Simulations from the posterior of parameters.
    ##
    if(output$sim.posterior){
      if(beta.size == 1) {
        if(sigmasq.prior == "fixed") {
          if(beta.prior != "fixed")
            kb.results$posterior$beta$sample <- rnorm(n.sim) * as.vector(sqrt(temp.post$beta.var)) + as.vector(temp.post$beta.mean)
        }
        else{
          kb.results$posterior$sigmasq$sample <- rinvchisq(n.sim, df.model, temp.post$S2)
          if(beta.prior != "fixed"){
            cond.beta.sd <- sqrt((as.vector(temp.post$beta.var) * kb.results$posterior$sigmasq$sample)/temp.post$S2)
            kb.results$posterior$beta$sample <- rnorm(n.sim) * cond.beta.sd + as.vector(temp.post$beta.mean)
          }
        }
      }
      else {
        if(sigmasq.prior == "fixed") {
          if(beta.prior != "fixed")
            kb.results$posterior$beta$sample <- array(apply(temp.post$beta.var,3,multgauss),dim=c(beta.size, n.sim))+temp.post$beta.mean
        }
        else {
          kb.results$posterior$sigmasq$sample <- rinvchisq(n.sim, df.model, temp.post$S2)
          if(beta.prior != "fixed"){
            if(is.R()) cond.beta.var <- temp.post$beta.var *rep(kb.results$posterior$sigmasq$sample/temp.post$S2,rep(beta.size^2,n.sim))
            else cond.beta.var <- temp.post$beta.var *rep(kb.results$posterior$sigmasq$sample/temp.post$S2,each = beta.size^2)
            kb.results$posterior$beta$sample <- array(apply(cond.beta.var,3,multgauss),dim=c(beta.size, n.sim)) + temp.post$beta.mean
          }
        }
      }
    }
    remove("temp.post")
  }
  if(output$keep.mcmc.sim) kb.results$posterior$simulations <- BC.inv(gauss.post, lambda)
  kb.results$model <- model
  kb.results$prior <- prior$priors.info
  kb.results$mcmc.input <- mcmc.input
  kb.results$.Random.seed <- seed
  kb.results$call <- call.fc
  attr(kb.results, "prediction.locations") <- call.fc$locations
  if(do.prediction) attr(kb.results, 'sp.dim') <- ifelse(krige1d, "1d", "2d")
  if(!is.null(call.fc$borders)) attr(kb.results, "borders") <- call.fc$borders
  class(kb.results) <- "glm.krige.bayes"
  return(kb.results)
}


"mcmc.aux" <- 
  function(z, data, meanS, QQ, Htrunc, S.scale, nsim, thin, QtivQ)
{
  ##
###### ------------------------ doing the mcmc-steps ----------- ############# 
  ##
  n <- length(data)
  result <- .C("mcmc1poislog",
               as.integer(n),
               z = as.double(z),
               S = as.double(rep(0, nsim * n)),
               as.double(data),
               as.double(meanS),
               as.double(as.vector(t(QQ))),
               as.double(as.vector(QtivQ)),
               as.double(rnorm(n * nsim * thin) * sqrt(S.scale)),
               as.double(runif(nsim * thin)),
               as.double(Htrunc),
               as.double(S.scale),
               as.integer(nsim),
               as.integer(thin),
               acc.rate = as.double(1), DUP=FALSE, PACKAGE = "geoRglm")[c("z", "S", "acc.rate")]
  attr(result$S, "dim") <- c(n, nsim)
  return(result)
}


"mcmc.pois.log" <- 
  function(data, units.m, meanS, invcov, mcmc.input, messages.screen)
{
  ## This is the MCMC engine for the spatial Poisson log Normal model ----
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  if(any(mcmc.input$Htrunc=="default")) Htrunc <- 2*data + 5
  else {
    if(is.vector(mcmc.input$Htrunc) & length(mcmc.input$Htrunc) == n)
      Htrunc <- mcmc.input$Htrunc
    else Htrunc <- rep(mcmc.input$Htrunc, n)
  }
  QQ <- t(chol(solve(invcov + diag(data))))
  sqrtdataQ <- sqrt(data)*QQ 
  QtivQ <- diag(n)-crossprod(sqrtdataQ)
  if(any(mcmc.input$S.start=="default")) {
    z <- as.vector(solve(QQ,ifelse(data > 0, log(data), -1.96) - meanS - log(units.m)))
  }
  else{
    if(any(mcmc.input$S.start=="random")) z <- rnorm(n)
    else{
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        else z <- as.vector(solve(QQ,mcmc.input$S.start))
      }
      else  stop(" S.start must be a vector of same dimension as data ")
    }
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
## ---------------- burn-in ----------------- ######### 
  if(burn.in > 0) {
    mcmc.output <- mcmc.aux(z, data, meanS + log(units.m), QQ, Htrunc, S.scale, 1, burn.in, QtivQ)
    if(messages.screen) cat(paste("burn-in = ", burn.in, " is finished. Acc.-rate = ", round(mcmc.output$acc.rate, digits=3), "\n"))
    acc.rate.burn.in <- c(burn.in, mcmc.output$acc.rate)
  }
  else mcmc.output <- list(z = z)
##### ---------- sampling periode ----------- ###### 
  if(n.iter <= 1000) {
    n.temp <- round(n.iter/thin)
    n.turn <- 1
  }
  else {
    n.temp <- round(1000/thin)
    n.turn <- round(n.iter/1000)
  }
  n.sim <- n.turn * n.temp
  Sdata <- matrix(NA, n, n.sim)
  acc.rate <- matrix(NA, n.turn, 2)
  for(i in seq(length=n.turn)) {
    mcmc.output <- mcmc.aux(mcmc.output$z, data, meanS + log(units.m), QQ, Htrunc, S.scale, n.temp, thin, QtivQ)
    Sdata[, seq((n.temp * (i - 1) + 1),(n.temp * i))] <- mcmc.output$S+meanS
    if(messages.screen) cat(paste("iter. numb.", i * n.temp * thin+burn.in, " : Acc.-rate = ", round(mcmc.output$acc.rate, digits=3), "\n"))
    acc.rate[i,1] <-  i * n.temp * thin
    acc.rate[i,2] <- mcmc.output$acc.rate
  }
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n"))
  if(burn.in > 0) acc.rate <- as.data.frame(rbind(acc.rate.burn.in,acc.rate))
  else acc.rate <- as.data.frame(acc.rate)
  names(acc.rate) <- c("iter.numb", "Acc.rate")
#########
  return(list(Sdata=Sdata, acc.rate=acc.rate))
}


"mcmc.boxcox.aux" <- 
  function(z, data, units.m, meanS, QQ, Htrunc, S.scale, nsim, thin, QtivQ, lambda)
{
  ##
###### ------------------------ doing the mcmc-steps ----------- ############# 
  ##
  n <- length(data)
  result <- .C("mcmc1poisboxcox",
               as.integer(n),
               z = as.double(z),
               S = as.double(rep(0, nsim * n)),
               as.double(data),
               as.double(units.m),
               as.double(meanS),
               as.double(as.vector(t(QQ))),
               as.double(as.vector(QtivQ)),
               as.double(rnorm(n * nsim * thin) * sqrt(S.scale)),
               as.double(runif(nsim * thin)),
               as.double(Htrunc),
               as.double(S.scale),
               as.integer(nsim),
               as.integer(thin),
               as.double(lambda),
               acc.rate = as.double(1), DUP=FALSE, PACKAGE = "geoRglm")[c("z", "S", "acc.rate")]
  attr(result$S, "dim") <- c(n, nsim)
  return(result)
}

"mcmc.pois.boxcox" <- 
  function(data, units.m, meanS, invcov, mcmc.input, messages.screen, lambda)
{
  ## This is the MCMC engine for the spatial Poisson - Normal model with link from the box-cox-family ----
  ##
  n <- length(data)
  S.scale <- mcmc.input$S.scale
  fisher.l <- ifelse(data>0,data^(1-2*lambda)*units.m^(2*lambda),0)
  QQ <- t(chol(solve(invcov + diag(fisher.l)))) 
  sqrtfiQ <- sqrt(fisher.l)*QQ 
  QtivQ <- diag(n)-crossprod(sqrtfiQ)
  if(any(mcmc.input$S.start=="default")) {
    S <- as.vector(ifelse(data > 0, (data/units.m)^lambda-1, -1.96)/lambda - meanS )       
    z <- as.vector(solve(QQ,S))
  }
  else{
    if(any(mcmc.input$S.start=="random")) z <- rnorm(n)
    else{
      if(is.numeric(mcmc.input$S.start)){
        if(length(mcmc.input$S.start) != n) stop("dimension of mcmc-starting-value must equal dimension of data")
        else z <- as.vector(solve(QQ,mcmc.input$S.start))
      }
      else  stop(" S.start must be a vector of same dimension as data ")
    }
  }
  if(any(mcmc.input$Htrunc=="default")) Htrunc <- 2*data + 5
  else {
    if(is.vector(mcmc.input$Htrunc) & length(mcmc.input$Htrunc) == n)
      Htrunc <- mcmc.input$Htrunc
    else Htrunc <- rep(mcmc.input$Htrunc, n)
  }
  burn.in <- mcmc.input$burn.in
  thin <- mcmc.input$thin
  n.iter <- mcmc.input$n.iter
  ## ---------------- burn-in ----------------- ######### 
  if(burn.in > 0) {
    mcmc.output <- mcmc.boxcox.aux(z, data, units.m, meanS, QQ, Htrunc, S.scale, 1, burn.in, QtivQ, lambda)
    if(messages.screen) cat(paste("burn-in = ", burn.in, " is finished. Acc.-rate = ", round(mcmc.output$acc.rate, digits=3), "\n"))
    acc.rate.burn.in <- c(burn.in, mcmc.output$acc.rate)
  }
  else mcmc.output <- list(z = z)
##### ---------- sampling periode ----------- ###### 
  if(n.iter <= 1000) {
    n.temp <- round(n.iter/thin)
    n.turn <- 1
  }
  else {
    n.temp <- round(1000/thin)
    n.turn <- round(n.iter/1000)
  }
  n.sim <- n.turn * n.temp
  Sdata <- matrix(NA, n, n.sim)
  acc.rate <- matrix(NA, n.turn, 2)
  for(i in seq(length=n.turn)){
    mcmc.output <- mcmc.boxcox.aux(mcmc.output$z, data, units.m, meanS, QQ, Htrunc, S.scale, n.temp, thin, QtivQ, lambda)
    Sdata[, seq((n.temp * (i - 1) + 1),(n.temp * i))] <- mcmc.output$S+meanS
    if(messages.screen) cat(paste("iter. numb.", i * n.temp * thin+burn.in, " : Acc.-rate = ", round(mcmc.output$acc.rate, digits=3), "\n"))
    acc.rate[i,1] <-  i * n.temp * thin
    acc.rate[i,2] <- mcmc.output$acc.rate
  }
  if(messages.screen) cat(paste("MCMC performed: n.iter. = ", n.iter, "; thinning = ", thin, "; burn.in = ", burn.in, "\n"))
  if(burn.in > 0) acc.rate <- as.data.frame(rbind(acc.rate.burn.in,acc.rate))
  else acc.rate <- as.data.frame(acc.rate)
  names(acc.rate) <- c("iter.numb", "Acc.rate")
  remove("z")
#########
  return(list(Sdata=Sdata, acc.rate=acc.rate))
}

"krige.glm.control" <-
  function (type.krige = "sk", trend.d = "cte", trend.l = "cte", obj.model = NULL, beta, cov.model, cov.pars, kappa,
            nugget, micro.scale, dist.epsilon = 1e-10, aniso.pars, lambda)
{
  if(type.krige != "ok" & type.krige != "OK" & type.krige != "o.k." & type.krige != "O.K." & type.krige != "sk" & type.krige != "SK" & type.krige != "s.k." & type.krige != "S.K.")
    stop("pois.krige: wrong option in the argument type.krige. It should be \"sk\" or \"ok\"(if ordinary or simple kriging is to be performed)")
  if(type.krige=="OK" | type.krige=="O.K." |type.krige=="o.k.")
    type.krige <- "ok"
  if(type.krige=="SK" | type.krige=="S.K." |type.krige=="s.k.")
    type.krige <- "sk"
  ##
  if(!is.null(obj.model)){
    if(missing(beta)) beta <- obj.model$beta
    if(missing(cov.model)) cov.model <- obj.model$cov.model
    if(missing(cov.pars)) cov.pars <- obj.model$cov.pars
    if(missing(kappa)) kappa <- obj.model$kappa
    if(missing(nugget)) nugget <- obj.model$nugget
    if(missing(micro.scale)) micro.scale <- nugget
    if(missing(lambda)) lambda <- obj.model$lambda
    if(missing(aniso.pars)) aniso.pars <- obj.model$aniso.pars
  }
  else{
    if(missing(beta)) beta <- NULL
    if(missing(cov.model)) cov.model <- "matern"
    if(missing(cov.pars)) stop("covariance parameters (sigmasq and phi) should be provided")
    if(missing(kappa)) kappa <- 0.5
    if(missing(nugget)) nugget <- 0
    if(missing(micro.scale)) micro.scale <- nugget
    if(missing(lambda)) lambda <- 0
    if(missing(aniso.pars)) aniso.pars <- NULL
  }
  ##
  if(type.krige == "sk")
    if(is.null(beta) | !is.numeric(beta))
      stop(" argument beta must be provided in order to perform simple kriging")
  if(micro.scale > nugget)
    stop(" micro.scale must be in the interval [0, nugget]")
  if(!is.null(aniso.pars))
    if(length(aniso.pars) != 2 | !is.numeric(aniso.pars))
      stop(" anisotropy parameters must be provided as a numeric vector with two elements: the rotation angle (in radians) and the anisotropy ratio (a number greater than 1)")
  ##
  if(inherits(trend.d, "formula") | inherits(trend.l, "formula")){
    if(!inherits(trend.d, "formula") | !inherits(trend.l, "formula"))
      stop(" trend.d and trend.l must have similar specification")
  }
  else{
    if((class(trend.d)=="trend.spatial") & (class(trend.l)=="trend.spatial")){
      if(ncol(trend.d) != ncol(trend.l))
        stop("pois.krige: trend.d and trend.l do not have the same number of columns")
    }
    else{
      if(trend.d != trend.l)
        stop(" trend.l is different from trend.d")
    }
  }
  cov.model <- match.arg(cov.model,
                         choices = c("matern", "exponential","gaussian",
                           "spherical", "circular", "cubic",
                           "wave", "power",
                           "powered.exponential", "cauchy", "gneiting",
                           "gneiting.matern", "pure.nugget"))
  if(cov.model == "power") stop("krige.glm.control: correlation function does not exist for the power variogram")
  res <- list(type.krige = type.krige,
              trend.d = trend.d, trend.l = trend.l, 
              beta = beta,
              cov.model = cov.model, 
              cov.pars = cov.pars, kappa = kappa,
              nugget = nugget,
              micro.scale = micro.scale, dist.epsilon = dist.epsilon, 
              aniso.pars = aniso.pars, lambda = lambda)
  class(res) <- "krige.geoRglm"
  return(res)
}

"krige.glm.check.aux" <-
  function(krige,fct)
{
  if(class(krige) != "krige.geoRglm"){
    if(!is.list(krige))
      stop(paste(fct,": the argument krige only takes a list or an output of the function krige.glm.control"))
    else{
      krige.names <-c("type.krige","trend.d","trend.l","obj.model","beta","cov.model",
                      "cov.pars","kappa","nugget","micro.scale","dist.epsilon","lambda","aniso.pars")
      krige <- object.match.names(krige,krige.names)
      if(is.null(krige$type.krige)) krige$type.krige <- "sk"  
      if(is.null(krige$trend.d)) krige$trend.d <-  "cte"
      if(is.null(krige$trend.l)) krige$trend.l <-  "cte"
      if(is.null(krige$cov.model)) krige$cov.model <- "matern"
      if(is.null(krige$kappa)) krige$kappa <-  0.5
      if(is.null(krige$nugget)) krige$nugget <-  0
      if(is.null(krige$micro.scale)) krige$micro.scale <- krige$nugget
      if(is.null(krige$dist.epsilon)) krige$dist.epsilon <-  1e-10
      krige <- krige.glm.control(type.krige = krige$type.krige,	
                                 trend.d = krige$trend.d, trend.l = krige$trend.l,
                                 obj.model = krige$obj.model,
                                 beta = krige$beta, cov.model = krige$cov.model,
                                 cov.pars = krige$cov.pars, kappa = krige$kappa,
                                 nugget = krige$nugget, micro.scale = krige$micro.scale,
                                 dist.epsilon = krige$dist.epsilon, 
                                 aniso.pars = krige$aniso.pars)
    }
  }
  return(krige)
}


"pois.krige" <- 
function(geodata, coords = geodata$coords, data = geodata$data, units.m = "default", locations = NULL,  borders = NULL, mcmc.input, krige, output)
{
  if(missing(geodata))
    geodata <- list(coords=coords, data=data, units.m=units.m)
  call.fc <- match.call()
  n <- length(data)
  if(any(units.m == "default")){
    if(!is.null(geodata$units.m)) units.m <- geodata$units.m
    else units.m <- rep(1, n)
  }
  if(missing(krige)) stop("must provide object krige")
  krige <- krige.glm.check.aux(krige,fct="pois.krige")
  cov.model <- krige$cov.model
  kappa <- krige$kappa
  beta <- krige$beta
  cov.pars <- krige$cov.pars
  nugget <- krige$nugget
  micro.scale <- krige$micro.scale
  aniso.pars <- krige$aniso.pars
  trend.d <- krige$trend.d
  trend.l <- krige$trend.l
  dist.epsilon <- krige$dist.epsilon
  lambda <- krige$lambda
  if(krige$type.krige == "ok") beta.prior <- "flat"
  if(krige$type.krige == "sk") beta.prior <- "deg"
  if(missing(output)) output <- output.glm.control()
  output <- output.glm.check.aux(output, fct="pois.krige")
  sim.predict <- output$sim.predict
  messages.screen <- output$messages.screen
  ##
  if(is.vector(coords)) {
    coords <- cbind(coords, 0)
    warning("vector of coordinates: one spatial dimension assumed")
  }
  coords <- as.matrix(coords)
  dimnames(coords) <- list(NULL, NULL)
  ## Checking for 1D prediction 
  if(length(unique(locations[,1])) == 1 | length(unique(locations[,2])) == 1)
    krige1d <- TRUE
  else krige1d <- FALSE
  ##
  if(is.null(locations)) {
    if(messages.screen) cat(paste("locations need to be specified for prediction; prediction not performed \n"))
  }
  else {
    if(is.null(trend.l))
      stop("trend.l needed for prediction")
  }
  trend.data <- unclass(trend.spatial(trend=trend.d, geodata = geodata))
  beta.size <- ncol(trend.data)
  if(nrow(trend.data) != n) stop("length of trend is different from the length of the data")
  if(beta.prior == "deg")
    if(beta.size != length(beta))
      stop("size of mean vector is incompatible with trend specified") 
  if(beta.size > 1)
    beta.names <- paste("beta", (0:(beta.size-1)), sep="")
  else beta.names <- "beta"
  ##
  ## preparing for MCMC 
  ##
  if(missing(mcmc.input)) stop("pois.krige: argument mcmc.input must be given")
  mcmc.input <- mcmc.check.aux(mcmc.input, fct="pois.krige")
  ##
  if(beta.prior == "deg") mean.d <-  as.vector(trend.data %*% beta)
  else mean.d <- rep(0,n)
  if(!is.null(aniso.pars)) {
    invcov <- varcov.spatial(coords = coords.aniso(coords = coords, aniso.pars = aniso.pars), cov.model = cov.model, kappa = kappa, 
                             nugget = nugget, cov.pars = cov.pars, inv = TRUE, func.inv = "cholesky",
                             try.another.decomposition = FALSE)$inverse
  }
  else {
    invcov <- varcov.spatial(coords = coords, cov.model = cov.model, kappa = kappa, nugget = nugget, cov.pars = cov.pars,
                             inv = TRUE, func.inv = "cholesky", try.another.decomposition = FALSE)$inverse
  }
  ##
########################----- MCMC ------#####################
  ##
  if(beta.prior == "flat") {
    ivtt <- invcov%*%trend.data
    invcov <- invcov-ivtt%*%solve.geoR(crossprod(trend.data, ivtt),t(ivtt))
  }
  if(lambda == 0){
    intensity <- mcmc.pois.log(data = data, units.m = units.m, meanS = mean.d, invcov=invcov, mcmc.input = mcmc.input, messages.screen=messages.screen)
    acc.rate <- intensity$acc.rate
    intensity <- exp(intensity$Sdata)
  }
  else{
    intensity <- mcmc.pois.boxcox(data=data, units.m=units.m, meanS=mean.d, invcov=invcov, mcmc.input=mcmc.input, messages.screen=messages.screen, lambda=lambda)
    acc.rate <- intensity$acc.rate
    intensity <- BC.inv(intensity$Sdata, lambda)    
  }
  ##
  ##------------------------------------------------------------
######################## ---- prediction ----- #####################
  if(!is.null(locations)) {
    if(!is.null(borders)){
      locations <- locations.inside(locations, borders)
      if(nrow(locations) == 0)
        stop(" pois.krige : there are no prediction locations inside the borders")
      if(messages.screen)
        cat(" pois.krige: results will be returned only for prediction locations inside the borders\n")
    }
    krige <- list(type.krige = krige$type.krige, beta = beta, trend.d = trend.d, trend.l = trend.l, cov.model = cov.model, 
                  cov.pars = cov.pars, kappa = kappa, nugget = nugget, micro.scale = micro.scale, dist.epsilon = dist.epsilon, 
                  aniso.pars = aniso.pars, lambda = lambda)
    kpl.result <- krige.conv.extnd(data = intensity, coords = coords, locations = locations, krige = krige,
                                   output = list(n.predictive = ifelse(sim.predict,1,0), signal = TRUE, messages = FALSE))
    remove(list = c("intensity"))
    kpl.result$krige.var <- rowMeans(kpl.result$krige.var) + apply(kpl.result$predict, 1, var) 
    if(nrow(locations) > 1) kpl.result$mcmc.error <- sqrt(asympvar(kpl.result$predict)/ncol(kpl.result$predict))
    else kpl.result$mcmc.error <- sqrt(asympvar(as.vector(kpl.result$predict), messages = FALSE)/length(as.vector(kpl.result$predict)))
    kpl.result$predict <- rowMeans(kpl.result$predict)
    if(beta.prior == "flat") {
      kpl.result$beta.est <- rowMeans(kpl.result$beta)
      names(kpl.result$beta.est) <- beta.names
    }
    kpl.result$beta <- NULL
  }
  else{
    if(beta.prior == "flat") {
      ## GLS
      beta.est <- solve.geoR(crossprod(trend.data, ivtt),t(ivtt))%*%rowMeans(log(intensity))
      kpl.result <- list(intensity=intensity, beta.est = beta.est, acc.rate=acc.rate)
    }
    else kpl.result <- list(intensity=intensity, acc.rate=acc.rate)
  }
  kpl.result$call <- call.fc
#######################################
  attr(kpl.result, "prediction.locations") <- call.fc$locations
  if(!is.null(locations)) attr(kpl.result, 'sp.dim') <- ifelse(krige1d, "1d", "2d")
  if(!is.null(call.fc$borders)) attr(kpl.result, "borders") <- call.fc$borders
  class(kpl.result) <- "kriging"
  return(kpl.result)
}

"proflik.glsm" <-
  function (mcmc.obj, obj.likfit.glsm, 
            phi.values, nugget.rel.values, messages, ...)
{
  ##
  ## Checking input
  ##
  geodata <- list(coords=mcmc.obj$coords)
  call.fc <- match.call()
  temp.list <- list()
  if(missing(messages))
    messages.screen <- ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages"))
  else messages.screen <- messages
  ##
  cov.model <- obj.likfit.glsm$cov.model
  kappa <- obj.likfit.glsm$kappa
  aniso.pars <- obj.likfit.glsm$aniso.pars
  lambda <- obj.likfit.glsm$lambda
  ##
  if(is.null(mcmc.obj$S)){
    if(is.null(mcmc.obj$mu)) stop("mcmc.obj should include either an object mu or an object S.")
    n <- temp.list$n <- nrow(mcmc.obj$mu)
    temp.list$mu <- mcmc.obj$mu
    est.boxcox <- TRUE
  }
  else{
    n <- temp.list$n <- nrow(mcmc.obj$S)
    temp.list$z <- mcmc.obj$S 
    est.boxcox <- FALSE
  }
  if(is.null(aniso.pars)) coords <- as.matrix(mcmc.obj$coords)
  else coords <- coords.aniso(coords = as.matrix(mcmc.obj$coords), aniso.pars = aniso.pars)
  temp.list$xmat <- obj.likfit.glsm$trend.matrix
  temp.list$xmat <- unclass(trend.spatial(trend=obj.likfit.glsm$trend, geodata = geodata))
  temp.list$beta.size <- dim(temp.list$xmat)[2]
  temp.list$coords <- coords
  temp.list$cov.model <- cov.model
  temp.list$kappa <- kappa
  proflik.val <- matrix(NA,length(phi.values),length(nugget.rel.values))
  if(messages.screen) cat("proflik.glsm: computing 2-D profile likelihood for the phi and relative nugget parameters \n")
  fixed.values <- list()
  ## where to find the parameter values :
  if(est.boxcox){
    fixed.values$lambda <- lambda
    ip <- list(f.tausq.rel = FALSE, f.lambda = TRUE)
  }
  else ip <- list(f.tausq.rel = FALSE)
  ##
  temp.list$log.f.sim <- mcmc.obj$log.f.sim
  temp.list$messages.screen <- getOption("verbose") ## we want a default FALSE here
  ##
  for(i in seq(along=phi.values)){
    for(j in seq(along=nugget.rel.values)){
      ini <- c(phi.values[i],nugget.rel.values[j])      
      if(est.boxcox) lik.val[i,j] <- lik.sim.boxcox(pars=ini, fp = fixed.values, ip = ip, temp.list = temp.list)
      else proflik.val[i,j] <- lik.sim(pars=ini, fp = fixed.values, ip = ip, temp.list = temp.list)
    }
  }
  result <- list()
  result$rangenugget.rel <- list(range=phi.values, nugget.rel = nugget.rel.values, proflik.rangenugget.rel = proflik.val, est.rangenugget.rel = obj.likfit.glsm$loglik)
  result$n.bi <- 1
  result$n.uni <-  0
  result$method.lik <- "ML"
  result$call <- call.fc
  class(result) <- "proflik" 
  return(result)
}

".First.lib" <- function(lib, pkg)
{
  messages.screen <- ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages"))
  if(messages.screen){
    cat("-----------------------------------------------------------\n")
    if(!require(geoR)){
      cat("\n")
      cat("Package geoR is required by geoRglm\n")
      cat("It should also be installed and loaded\n")
    }
  }
  library.dynam("geoRglm", package=pkg, lib.loc=lib)
  if(messages.screen){
    pkg.info <- packageDescription("geoRglm", lib.loc = lib, fields=c("Title","Version","Date"))
    cat(pkg.info$Title)
    cat("\n")
    cat(paste("geoRglm version ", pkg.info$Version, " (", pkg.info$Date, ") is now loaded\n", sep=""))
    cat("-----------------------------------------------------------\n")
    cat("\n")
  }
  return(invisible(0))
}
