#
#   Décembre  2007
#
#    modif de Michel Goulard pour faire du voisinage glissant
#
`binom.krige1` <-
function (geodata, coords = geodata$coords, data = geodata$data, 
    units.m = "default", locations = NULL, borders, mcmc.input, 
    krige, output,nwin="full",trend0=T) 
{
    if (missing(geodata)) 
        geodata <- list(coords = coords, data = data, units.m = units.m)
    if (missing(borders)) 
        borders <- geodata$borders
    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 (any(units.m <= 0)) 
        stop("units.m must be postive")
    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)
    if (is.null(locations)) {
        if (messages.screen) 
            cat(paste("locations need to be specified for prediction; prediction not performed \n"))
    }
    else {
        locations <- .check.locations(locations)
        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
    }
    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"
    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
    }
    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
    if (!is.null(locations)) {
      if(trend0){
        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")
      }
      else{
        krige <- list(type.krige = krige$type.krige, 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.aux1(data = res.mcmc$Sdata, coords = coords, 
            locations = locations, borders = borders, krige = krige, 
            output = list(n.predictive = ifelse(sim.predict, 
                1, 0), signal = TRUE, messages = TRUE),nwin=nwin)
        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)
	kpl.result$nvar=rowMeans(kpl.result$nvar)+ 
            apply(kpl.result$predict, 1, var)
	 kpl.result$npred = rowMeans(kpl.result$npred)
        if (beta.prior == "flat") {
#            kpl.result$beta.est <- rowMeans(kpl.result$beta.est)
#            kpl.result$mean.est <- rowMeans(kpl.result$mean.est)
            names(kpl.result$beta.est) <- beta.names
        }
        kpl.result$beta <- NULL
    }
    else {
        if (beta.prior == "flat") {
            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.aux1` <-
function (data, coords, locations, borders, krige, output,nwin="full") 
{
#
#   Décembre 2007 ajout voisinage glissant
#       et rajout récup du prédicteur  avant transfo.
#    le prédicteur de la proba est npred et sa variance de pred nvar
#
    krige$lambda <- 1
    krige$link <- NULL
    kc.result <- krige.conv.extnd1(data = data, coords = coords, 
        locations = locations, borders = borders, krige = krige, 
        output = output,nwin=nwin)
    kc.result$message <- NULL
    ivlogit2 <- ifelse(kc.result$predict < 700, exp(kc.result$predict) * 
        (-expm1(kc.result$predict))/(1 + exp(kc.result$predict))^3, 
        0)
    kc.result$npred <- plogis(kc.result$predict) + 0.5 * ivlogit2 * 
        kc.result$krige.var
    kc.result$nvar <- 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.conv.extnd1` <-
function (geodata, coords = geodata$coords, data = geodata$data, 
    locations, borders, krige, output,nwin="full") 
{
    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)
    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
    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
    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)
    locations <- .check.locations(locations)
    if (!is.null(borders)) {
        nloc0 <- nrow(locations)
        ind.loc0 <- .geoR_inout(locations, borders)
        locations <- locations[ind.loc0, , drop = FALSE]
        if (nrow(locations) == 0) 
            stop("\n .krige.conv.extnd : there are no prediction locations inside the borders")
        if (messages.screen) 
            cat("krige.conv: results will be returned only for prediction locations inside the borders\n")
    }
    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)))
    if (!is.null(borders) && nrow(trend.l) == nloc0) 
        trend.l <- trend.l[ind.loc0, , drop = FALSE]
    if (nrow(trend.l) != nrow(locations)) 
        stop("locations and trend.l have incompatible sizes")
    ni <- nrow(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)
    }
    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
    }
    tausq <- nugget
    if (is.vector(cov.pars)) {
        sigmasq <- cov.pars[1]
        phi <- cov.pars[2]
    }
    else stop("covariance parameter should be given as a vector")
    cpars <- c(1, phi)
    tausq.rel <- tausq/sigmasq
    nug.factor <- ifelse(signal, krige$micro.scale/sigmasq, tausq.rel)
    kc.result <- list()
    
#
#   modif voisinage unique
#
    if(nwin=="full"){
    out.message <- "krige.conv.extnd: Kriging performed using global neighbourhood"
    if (messages.screen) 
           cat(paste(out.message, "\n"))
    Vcov <- varcov.spatial(coords = coords, cov.model = cov.model, 
        kappa = kappa, nugget = tausq.rel, cov.pars = cpars)$varcov
    ivtt <- solve(Vcov, trend.data)
    ttivtt <- crossprod(ivtt, trend.data)
    if (beta.prior == "flat") 
        beta.flat <- solve(ttivtt, crossprod(ivtt, data))
    remove("ivtt")
    temp <- NULL
    v0 <- loccoords(coords = coords, locations = locations)
    if (n.predictive > 0) {
        loc.coincide <- apply(v0, 2, function(x, min.dist) {
            any(x < min.dist)
        }, min.dist = krige$dist.epsilon)
        if (any(loc.coincide)) 
            loc.coincide <- (1:ni)[loc.coincide]
        else loc.coincide <- NULL
        if (!is.null(loc.coincide)) {
            temp.f <- function(x, data, dist.eps) {
                return(data[x < dist.eps, , drop = FALSE])
            }
            data.coincide <- apply(v0[, loc.coincide, drop = FALSE], 
                2, temp.f, data = data, dist.eps = krige$dist.epsilon)
        }
        else data.coincide <- NULL
    }
    ind.v0 <- which(v0 < krige$dist.epsilon)
    v0 <- cov.spatial(obj = v0, cov.model = cov.model, kappa = kappa, 
        cov.pars = cpars)
    v0[ind.v0] <- 1 + nug.factor
    ivv0 <- solve(Vcov, v0)
    tv0ivv0 <- colSums(v0 * ivv0)
    b <- t(trend.l - crossprod(ivv0, trend.data))
    tv0ivdata <- crossprod(ivv0, data)
    if (n.predictive == 0) {
        remove(list = "v0")
        gc(verbose = FALSE)
    }
    if (beta.prior == "deg") {
        kc.result$predict <- array(tv0ivdata + rep(drop(crossprod(b, 
            beta)), n.datasets), dim = c(ni, n.datasets))
        if (n.predictive == 0) 
            remove(list = c("b", "tv0ivdata"))
        else remove("tv0ivdata")
        kc.result$krige.var <- sigmasq * drop(1 + nug.factor - 
            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)) + 
            crossprod(b, beta.flat)
        remove(list = c("tv0ivdata"))
        bitb=colSums(b * solve(ttivtt, b))
        kc.result$krige.var <- sigmasq * drop(1 + nug.factor - 
            tv0ivv0 + bitb)
        kc.result$beta.est <- apply(beta.flat,1,mean)
        remove("beta.flat")
    }
}
#
#   voisinage glissant nwin est numérique
#
else{# nwin numérique
        out.message <- "krige.conv.extnd: Kriging performed using local neighbourhood"
        if (messages.screen) cat(paste(out.message, "\n"))
	kc.result <- list(predict = matrix(0,ni,n.datasets), krige.var = rep(0,ni), var=rep(0,ni), mean.est=rep(0,ni),beta.est=matrix(0,ni,ncol(trend.l)))
	v0tot <- loccoords(coords = coords, locations = locations)
	for (i in 1:ni) {
	      if (messages.screen){
        	if(i%%100==0) cat(paste("krige.conv: kriging location: ", i, "out of", ni, "\n"))
              }
      ##
      ## choix de la fenêtre
      ## 
	dm0 <- v0tot[,i]
	coordswin <- coords[order(dm0)[1:nwin], ]
	trend.dwin <- trend.data[order(dm0)[1:nwin],]
	trend.li<- trend.l[i,]
	trend.li <- matrix(trend.li,nrow=1)
	datawin <- data[order(dm0)[1:nwin],]
	Vcov <- varcov.spatial(coords = coordswin, cov.model = cov.model, 
        kappa = kappa, nugget = tausq.rel, cov.pars = cpars)$varcov
	ivtt <- solve(Vcov, trend.dwin)
    	ttivtt <- crossprod(ivtt, trend.dwin)
        if (beta.prior == "flat") 
          beta.flat <- solve(ttivtt, crossprod(ivtt, datawin))
        remove("ivtt")
        temp <- NULL
	dm0=dm0[order(dm0)[1:nwin]]
    	v0 <- cov.spatial(obj = dm0, cov.model = cov.model, kappa = kappa, cov.pars = cpars)
    	ind.v0 <- which(dm0 < krige$dist.epsilon)
	v0[ind.v0] <- 1 + nug.factor
	ivv0 <- solve(Vcov, v0)
    	tv0ivv0 <- sum(v0 * ivv0)
	tv0ivx=crossprod(ivv0, trend.dwin)
     	b <- trend.li - tv0ivx
	b=as.double(b)
    	tv0ivdata <- crossprod(ivv0, datawin)
    if (beta.prior == "deg") {
        kc.result$predict[i,] <- tv0ivdata + rep(crossprod(b,beta),n.datasets)
        kc.result$krige.var[i] <- sigmasq * (1 + nug.factor - tv0ivv0)
        kc.result$beta.est <- paste(".krige.conv.extnd: Simple kriging (beta provided by user)\n")
        kc.result$mean.est <- paste(".krige.conv.extnd: Simple kriging (beta provided by user)\n")
    }
    if (beta.prior == "flat") {
        kc.result$predict[i,] <- tv0ivdata + crossprod(b,beta.flat)
        bitb=sum(b * solve(ttivtt, b))
        kc.result$krige.var <- sigmasq * drop(1 + nug.factor - 
            tv0ivv0 + bitb)
        kc.result$mean.est[i] <- mean(crossprod(matrix(trend.li,ncol=1),beta.flat))
	kc.result$beta.est[i,]=apply(beta.flat,1,mean)
        remove("beta.flat")

    }
}
}
#
# fin modif nwin
#    
    kc.result$krige.var[kc.result$krige.var < 1e-12] <- 0
    if (any(kc.result$krige.var < 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"))
    if (n.predictive > 0) {
        if (messages.screen) 
            cat(".krige.conv.extnd: sampling from the predictive distribution (conditional simulations)\n")
        Dval <- 1 + nug.factor
        if (beta.prior == "deg") 
            vbetai <- matrix(0, ncol = beta.size, nrow = beta.size)
        else vbetai <- matrix(.solve.geoR(ttivtt), ncol = beta.size, 
            nrow = beta.size)
        coincide.cond <- (((round(1e+12 * 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) {
            invcov <- varcov.spatial(coords = coords, cov.model = cov.model, 
                kappa = kappa, nugget = tausq.rel, cov.pars = cpars, 
                inv = TRUE, only.inv.lower.diag = TRUE)
            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"))
    }
    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) {
            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 = out.message, call = cl))
    return(kc.result)
}
binom.krige.simul1=function(data, coords, units.m,mcmc.input, cov.model,
	kappa,beta,cov.pars,nugget,aniso.pars,trend.data,
	beta.prior,messages.screen,nwin)
{
    n=nrow(coords)	
    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
    }
    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)
}