Kriging <- function(krige.method, x, y=NULL, z=NULL, T=NULL,
                    grid, gridtriple=FALSE,
                    model, param, given, data, trend,                 
                    pch=".") {

  krige.methodlist <- c("S", "O")
  if (is.na(krige.method.nr <- pmatch(krige.method,krige.methodlist)))
    stop(paste("kriging method not identifiable from the list.",
               "Possible values are", paste(krige.methodlist, collapse=",")))

  x <- CheckXT(x=x, y=y, z=z, T=T, grid=grid, gridtriple=gridtriple)
  y <- z <- NULL
  pm <- PrepareModel(model=model, param=param,
                     timespacedim=x$spacedim + x$Time,
                     trend=trend)

  is.matrix.data <- is.matrix(data) && (ncol(data)>1)
  data <- as.matrix(data)
  given <- as.matrix(given)
  xdim <- ncol(given)

  if (pm$timespacedim!= xdim)
    stop("dimensions of the kriged points and the given points do not match")  
  stopifnot(pm$timespacedim == ncol(x$x) + !is.null(x$T)) ## program wrong

  nd <- nrow(given)
  pos <- integer(nd)
  ## lexicographical ordering of vectors --- necessary to check
  ## whether any location is given twice, but with different value of
  ## the data
  .C("Ordering", as.double(t(given)), as.integer(nd), as.integer(xdim), pos,
     PACKAGE="RandomFields", DUP=FALSE)
  pos <- pos + 1

  given <- given[pos, , drop=FALSE]
  data <- data[pos, , drop=FALSE]
  
  ## are locations given twice with different data values?
  dup <- c(FALSE, apply(abs(given[-nd, , drop=FALSE] -
                            given[-1, , drop=FALSE]), 1, sum))
  if (any(dup <- c(FALSE, apply(abs(given[-1, , drop=FALSE] -
                                    given[-nd, , drop=FALSE]), 1, sum)==0))) {
    if (any(data[dup, ] != data[c(dup[-1], FALSE), ]))
      stop("duplicated conditioning points with different data values")
    given <- given[!dup, , drop=FALSE]
    data <- data[!dup, , drop=FALSE]  
  }
  
  tgiven <-  t(given)
  nd <- nrow(given)

  ## the next passage might be replaced by a direct call of
  ## covarianceMatrix...
  covmatrix <- diag(nd) *
    (CovarianceFct(if (pm$anisotropy) t(rep(0, pm$timespacedim)) else 0,
                   model, param) / 2)
  low.tri <- lower.tri(covmatrix, diag=FALSE)
  covmatrix[low.tri] <-
    CovarianceFct(if (pm$anisotropy)
                  matrix(.C("vectordist", as.double(given),
                            as.integer(dim(given)),
                            vd=double(xdim * nd * (nd - 1)/2), as.integer(FALSE),
                            PACKAGE="RandomFields")$vd, ncol=xdim) else 
                  as.matrix(dist(given))[low.tri],
                  model, param)
  covmatrix <- covmatrix + t(covmatrix)
  given <- NULL

  if (grid) {
    zz <- cbind(x$x, x$T)
    eval(parse(text=paste("dimension <- c(",
                 paste(paste("length(seq(",zz[1,],",", zz[2,],",",
                             zz[3,],"))"), collapse=","),
                 ")")))
    ## `x' will be used in apply within kriging
    if ( (l <- ncol(zz))==1 ) x <- matrix(seq(zz[1],zz[2],zz[3]),ncol=1)
    else {
      text <- paste("x <- expand.grid(",
                    paste("seq(zz[1,",1:l,
                          "],zz[2,",1:l,
                          "],zz[3,",1:l,"])",collapse=","),
                    ")")
      eval(parse(text=text))
    }
  } else x <- x$x

  env <- environment()
  nd.step <-  ceiling(nrow(x) / 79) ## for the user's entertainment
  assign("enu", 0, envir=env)       ## ="=
  nn <- as.integer(ncol(tgiven))
  error <- integer(1)
  .C("InitUncheckedCovFct",
     as.integer(pm$covnr),
     as.double(pm$param),
     as.integer(length(pm$param)),
     as.integer(pm$timespacedim),
     as.integer(xdim),
     as.integer(length(pm$covnr)),
     as.integer(pm$anisotropy),
     as.integer(pm$op),
     as.integer(RFparameters()$PracticalRange),
     error, PACKAGE="RandomFields", DUP=FALSE);
  if (error) stop(" Error in definition of covariance function")

  switch(krige.method.nr,
         { ## simple kriging#
           stopifnot(is.null(pm$trend))
           if (!(is.matrix(try(data <- as.matrix(solve(covmatrix,
                                                       data-pm$mean)))))) {
             res <- x * NA
           } else {
             ## time consuming "apply" -- to do : replace by c call
             res <- apply(x, 1,
                          function(z){
                            if ((enu %% nd.step) == 0) cat(pch);
                            assign("enu", enu+1, envir=env);
                            .C("UncheckedCovFct", tgiven - z, nn, cov=double(nn),
                               PACKAGE="RandomFields", DUP=FALSE)$cov %*% data}
                          ) + pm$mean
           }
         }, {
           ## ordinary kriging
           stopifnot(is.null(pm$trend))
           covmatrix <- rbind(cbind(covmatrix,1), c(rep(1,nd),0)) 
            if (!(is.matrix(try(data <- solve(covmatrix,
                                              rbind(as.matrix(data),0)))))) {
              res <- x * NA
            } else {
              ## time consuming "apply" -- to do : replace by c call
              res <-
                apply(x, 1,
                      function(z){
                        if ((enu %% nd.step) == 0) cat(pch);
                        assign("enu", enu+1, envir=env);
                        c(.C("UncheckedCovFct", tgiven-z, nn, cov=double(nn),
                             PACKAGE="RandomFields", DUP=FALSE)$cov,
                          1) %*% data
                      })
            }
         } , {
           ## universal kriging
           stopifnot(!is.null(trend))
           stop("not programmed yet")
         }
         ) # switch

  if (pch!="") cat("\n")
  ncol.data <- ncol(data)
  x <- data <- NULL
  
  if (is.matrix.data) {
    if (grid) return(array(t(res), dim=c(dimension, ncol.data)))
    else return(t(res))
  } else {
    if (grid) return(array(res, dim=dimension))
    else return(res)
  }
}


CondSimu <- function(krige.method, x, y=NULL, z=NULL, T=NULL,
                     grid, gridtriple=FALSE,
                     model, param, method=NULL,
                     given, data, trend,
                     n=1, register=0, 
                     err.model=NULL, err.param=NULL, err.method=NULL,
                     err.register=1, 
                     tol=1E-5, pch=".",
                     na.rm=FALSE
                     ) {
  op.list <- c("+","*")
  
  if (is.character(method) && (!is.na(pmatch(method, c("S","O")))))
    stop("Sorry. The parameters of the function `CondSimu' as been redefined. Use `krige.method' instead of `method'. See help(CondSimu) for details")

  x  <- CheckXT(x, y, z, T, grid, gridtriple)
  pm <- PrepareModel(model, param, x$spacedim + x$Time, trend, method=method)
  y <- z <- NULL   
  total <- x$total

  krige.mean <- 0 ## pm$mean + err$mean ??
  krige.trend <- pm$trend
  if (!is.null(krige.trend)) stop("not programmed yet")

  if (!is.null(err.model)) {
    err <- PrepareModel(err.model, err.param,x$spacedim+x$Time,method=err.method)
    if (xor(pm$anisotropy, err$anisotropy))
      stop("both, the data model and the error model must be either istropic or anisotropic")

    if (xor(is.null(pm$method), is.null(err$method)))
      stop("the simulation method must be defined for the data model and the error model or for none of them")
    if (((length(err$cov)>1 &&
          err$cov!=.C("GetModelNr", "nugget", 1, nr=integer(1),
            PACKAGE="RandomFields")$nr
          ) || err$timespacedim>1)
        && (err$timespacedim != pm$timespacedim))
      stop("time-space dimensions do not match") 
    if (!is.null(err$trend)) stop("trend for the error term not allowed")
    
    krige <- convert.to.readable(list(covnr=c(pm$covnr, err$covnr),
                                     param=c(pm$param, err$param),
                                     anisotropy = pm$anisotropy,
                                     op = c(pm$op, pmatch("+", op.list) - 1,
                                       err$op),
                                     mean = krige.mean,  
                                     trend = krige.trend,
                                     method = c(pm$method, err$method),
                                     timespacedim = pm$timespacedim
                                     )
                                )
  } else krige <-convert.to.readable(list(covnr=pm$covnr,
                                     param=pm$param,
                                     anisotropy = pm$anisotropy,
                                     op = pm$op,
                                     mean = krige.mean,  
                                     trend = krige.trend, 
                                     method = pm$method,
                                     timespacedim = pm$timespacedim
                                     )
                                )

  simu.grid <- grid
  given <- as.matrix(given)
  if (nrow(given)!=length(data)) {
    cat("dimension of 'given' does not match 'data'")
    return(NA)
  }
  if (na.rm && any(data.idx <- is.na(data))) {
    data <- data[data.idx]
    given <- given[data.idx, ]
  }
  
  txt <- "kriging in space time dimensions>3 where not all the point ly on a grid is not possible yet"
  ## if 4 dimensional then the last coordinates should ly on a grid

  ## now check whether and if so, which of the given points belong to the
  ## points where conditional simulation takes place
  if (grid) {
    zz <-  cbind(x$x, x$T)
    ind <- 1 + (t(given) - zz[1,]) / zz[3,]
    index <-  round(ind)
    endpoints <- 1 + floor((zz[2,]-zz[1,])/zz[3,])
    outside.grid <- apply((abs(ind-index)>tol) | (index<1) |
                          (index>endpoints), 2, any)   
    if (any(outside.grid)) {
      ## at least some data points are not on the grid
      ## simulate as if there is no grid
      simu.grid <- FALSE
      l <- ncol(zz)
      if (l>3) stop(txt)
      if (l==1) xx <- matrix(seq(x$x[1],x$x[2],x$x[3]),nrow=1)
      else  eval(parse(text=paste("xx <-  t(as.matrix(expand.grid(",
                         paste("seq(zz[1,",1:l,
                               "],zz[2,",1:l,
                               "],zz[3,",1:l,"])",collapse=","),
                         ")))")))
      eval(parse(text=paste("ll <- c(",
                   paste("length(seq(zz[1,",1:l,"],zz[2,",1:l,"],zz[3,",1:l,
                         "]))",
                         collapse=","),
                   ")")))

      new.index <- rep(0,ncol(index))
      ## data points that are on the grid, must be registered,
      ## so that they can be used as conditioning points of the grid
      if (!all(outside.grid)) {
        new.index[!outside.grid] <- 1 +
          apply((index[,!outside.grid,drop=FALSE]-1) *
                cumprod(c(1,ll[-length(ll)])), 2, sum)
      }
      index <- new.index
      new.index <- NULL
    } else {
      ## data points are all lying on the grid     
      z <- GaussRF(x=x$x, T=x$T,  grid=TRUE, model=model, param=param,
                   method=method, n=n, register=register,
                   gridtriple=TRUE)
      ## for all the other cases of simulation see, below
      index <- t(index)
    }
  } else {
    xx <- t(as.matrix(x$x)) ## not a grid
    
    if (!is.null(x$T)) {
      if (ncol(xx) > 2) stop(txt)
      T <- seq(x$T[1], x$T[2], x$T[3])
      ## multiple the xx structure by length of T,
      ## then add first component of T to first part ... last component of T
      ## to last part
      xx <- rbind(matrix(xx, nrow=nrow(xx), ncol=ncol(xx) * length(T)),
                  as.vector(t(matrix(T, nrow=length(T),ncol(xx)))))
    }
    
    ## the next step can be pretty time consuming!!!
    ## to.do: programme it in C
    ##
    ## identification of the points that are given twice, as points to
    ## be simulated and as data points (up to a tolerance distance !)
    ## this is important in case of nugget effect, since otherwise
    ## the user will be surprised not to get the value of the data at
    ## that point
    one2ncol.xx <- 1:ncol(xx)
    index <- apply(as.matrix(given),1,function(z){
      i <- one2ncol.xx[apply(abs(xx - z),2,sum) < tol]
      if (length(i)==0) return(0)
      if (length(i)==1) return(i)
      return(NA)
    })
  }  
 
  if (!simu.grid) {
    ## otherwise the simulation has already been performed (see above)
    tol <- tol * nrow(xx)
    if (any(is.na(index)))
      stop("identification of the given data points is not unique - `tol' too large?")
    if (any(notfound <- (index==0))) {
      index[notfound] <- (ncol(xx)+1):(ncol(xx)+sum(notfound))
    }
    xx <- rbind(t(xx), given[notfound,,drop=FALSE])
    z <- GaussRF(x=xx, grid=FALSE, model=model, param=param,
                 method=method, n=n, register=register)
    xx <- NULL
  }
  if (is.null(z)) stop("random fields simulation failed")
  
  if (n==1) {
    ## z values at the `given' locations
    zgiven <- z[index]
    z <- z[1:total]
  } else {
    ## this is a bit more complicated since index is either a vector or
    ## a matrix of dimension dim(z)-1
    zgiven <- apply(z, length(dim(z)), function(m) m[index])
    z <- as.vector(apply(z,length(dim(z)),function(m) m[1:total]))
  }
  
  if (!is.null(err.model)) {
     error <- GaussRF(given, grid=FALSE, model=err.model, param=err.param,
                     method=err.method, n=n, register=err.register)  
     if (is.null(error)) stop("error field simulation failed")
     zgiven <- zgiven + as.vector(error)
     error <- NULL
   }
  
  z +  Kriging(krige.method=krige.method,
              x=x$x, grid=grid,
              model=krige,
              given=given,data=as.vector(data)-zgiven,
              gridtriple=TRUE, pch=pch)
}