## Authors 
## Martin Schlather, schlather@math.uni-mannheim.de
##
##
## Copyright (C) 2015 -- 2017 Martin Schlather
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License
## as published by the Free Software Foundation; either version 3
## of the License, or (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.  

# source("rf.R")
# see getNset.R for the variable .methods



RFboxcox <- function(data, boxcox, vdim=1, inverse=FALSE, ignore.na=FALSE) {
  if (missing(boxcox)) boxcox <- .Call(C_get_boxcox)
  if (any(is.na(boxcox)) && !ignore.na)
    stop("non-finte values in Box-Cox transformation")
  if (!all(is.finite(boxcox))) return(data)
  if (is.list(data)) {
    for (i in 1:length(data))
      data[[i]] <- RFboxcox(data[[i]], boxcox=boxcox, vdim=vdim,
                            inverse=inverse, ignore.na=ignore.na)
    return(data)
  }
  Data <- data + 0
  .Call(C_BoxCox_trafo, as.double(boxcox), as.double(Data), as.integer(vdim),
        as.logical(inverse));
  return(Data)
}


RFlinearpart <- function(model, x, y = NULL, z = NULL, T=NULL, grid=NULL,
                         data, params,  distances, dim, set=0, ...) {
  Reg <- MODEL_USER  
  RFoptOld <- internal.rfoptions(..., RELAX=is(model, "formula"))
  on.exit(RFoptions(LIST=RFoptOld[[1]]))
  RFopt <- RFoptOld[[2]]

  model <- list("linearpart", PrepareModel2(model, params=params, ...))

  rfInit(model=model, x=x, y=y, z=z, T=T, grid=grid,
         distances=distances, dim=dim, reg = Reg, dosimulate=FALSE)

  .Call(C_get_linearpart, Reg, as.integer(set))
}



setvector <- function(model, preceding, len, factor) {
  if (model[[1]] == SYMBOL_PLUS) {
    return(c(SYMBOL_PLUS,
             lapply(model[-1], setvector, preceding=preceding, len=len,
                    factor=factor)))
  }
  if (found <- model[[1]] == R_C) {
    if (length(model) != len + 1) stop("bug")
    model <- c(model[1], rep(0, preceding), model[-1])
    if (!missing(factor)) model <- list(SYMBOL_MULT, model, factor)
  } else if (found <- model[[1]] == R_CONST) {
    if (length(model[[1]]) != len) stop("bug")
    model[[2]] <- c(rep(0, preceding), model[[2]])
    if (!missing(factor)) model <- list(SYMBOL_MULT, model, factor)
  } else if (model[[1]] == SYMBOL_MULT) {
    for (i in 2:length(model)) {
      if (found <- model[[i]][[1]]==R_C) {
        if (length(model[[i]]) != len + 1) stop("bug")
        model[[i]] <- c(R_C, rep(0, preceding), model[[i]][-1])
        if (!missing(factor)) model[[length(model) + 1]] <- factor
        break
      }      
    }
  }
  if (!found) {
    bind <- c(list(R_C), rep(0, preceding), rep(1, len))
    splittingC <- function(model, preceding, factor) {
      const <- sapply(model[-1],
                      function(m) (is.numeric(m) && !is.na(m)) ||
                                  (m[[1]] == R_CONST && !is.na(m[[2]]))
                      )
      if (all(const)) {
        model <- c(model[1], if (preceding > 0) rep(0, preceding), model[-1])
        return(list(SYMBOL_MULT, model, if (!missing(factor)) list(factor)))
      }
      for (i in 2:length(model)) {
        vdim <- preceding + (if (i==2) 0 else GetDimension(model[[i-1]]))
        m  <- ReplaceC(model[[i]])
        L <- GetDimension(m)
        model[[i]] <- setvector(m, preceding = vdim, len = L, factor=factor)
      }
      model[[1]] <- SYMBOL_PLUS
      names(model) <- NULL
      L <- GetDimension(model[[length(model)]])
      model <- SetDimension(model, L)
    }


    if (model[[1]] == SYMBOL_MULT)
      model <- c(SYMBOL_MULT, list(bind), model[-1])
    else model <- list(SYMBOL_MULT, bind, model)
    if (!missing(factor)) model[[length(model) + 1]] <- factor
  }
  return(model)
}


GetDimension <- function(model){
  if (model[[1]] == SYMBOL_PLUS) return(max(sapply(model[-1], GetDimension)))
  if (model[[1]] == R_C) return(length(model) - 1)
  if (model[[1]] == R_CONST) return(length(model[[2]]))
  if (model[[1]] == SYMBOL_MULT) 
    return(max(sapply(model, function(m) if (m[[1]]==R_C) length(m)-1 else 1)))
  return(1)
}


SetDimension <- function(model, L){
  if (model[[1]] == SYMBOL_PLUS) {
    return(c(SYMBOL_PLUS, lapply(model[-1], SetDimension, L=L)))
  }
  if (model[[1]] == R_C) {
    if (length(model) <= L) for (i in length(model) : L) model[[i+1]] <- 0
    names(model) <- c("", letters[1:L])
  } else if (model[[1]] == R_CONST) {
    if (length(model[[2]]) < L)
      model[[2]] <- c(model[[2]], rep(0, L - length(model[[2]])))
  } else if (model[[1]] == SYMBOL_MULT) {
    for (i in 2:length(model)) {
      if (model[[i]][[1]]==R_C) {
        if (length(model[[i]]) <= L)
          for (j in length(model[[i]]) : L) model[[i]][[j+1]] <- 0
        names(model[[i]]) <- c("", letters[1:L])
      }
    }
  }
  return(model)
}


splittingC <- function(model, preceding, factor) {
  const <- sapply(model[-1],
                  function(m) {
		    ((is.numeric(m) || is.logical(m)) && !is.na(m)) ||
                      (is.list(m) && m[[1]] == R_CONST && !is.na(m[[2]]))
		  })
  if (all(const)) {
    model <- c(model[1], if (preceding > 0) rep(0, preceding), model[-1])
    return(list(SYMBOL_MULT, model, if (!missing(factor)) list(factor)))
  }
  for (i in 2:length(model)) {
    vdim <- preceding + (if (i==2) 0 else GetDimension(model[[i-1]]))
    m <- model[[i]]
    if (is.list(m)) {
      m  <- ReplaceC(m)
      L <- GetDimension(m)
    } else if (is.numeric(m) || is.logical(m)) {
      L <- length(m)
      if (L ==1 && is.na(m)) {
	m <- list(SYMBOL_MULT, list(R_CONST, a=m)) 
      } else {
	m <- list(R_CONST, m)
      }
    } else stop(m, "not allowed")
    model[[i]] <- setvector(m, preceding = vdim, len = L, factor=factor)
  }
  model[[1]] <- SYMBOL_PLUS
  names(model) <- NULL
  L <- GetDimension(model[[length(model)]])
  model <- SetDimension(model, L)
}


SplitC <- function(model, factor) {
  model <- splittingC(model, 0, factor)
  L <- GetDimension(model)
  return(SetDimension(model, L))
}



AnyIsNA <- function(model) {
  if (is.list(model)) {
    for (i in 1:length(model)) if (AnyIsNA(model[[i]])) return(TRUE)
    return(FALSE)
  } else return((is.numeric(model) || is.logical(model)) && any(is.na(model)))
}


##RReplaceC <-
ReplaceC <- function(model) {
  if (model[[1]] == SYMBOL_PLUS) {
    for (i in 2:length(model)) model[[i]] <- ReplaceC(model[[i]])
  } else if (model[[1]] == SYMBOL_MULT) {
    if (length(model) <= 2) {
      stop("here, products must have at least 2 factors")
    }
    cs <- sapply(model, function(m) m[[1]] == R_C)
    s <- sum(cs)
    if (s > 0) {
      if (s > 1)
        stop("multiplication with '", R_C, "' may happen only once")
      cs <- which(cs)
      C <- model[[cs]]      
      if (!AnyIsNA(model[-cs])) {
        return(SplitC(C, factor=model[-cs]))
      }
    }      
  } else if (model[[1]] == R_C) {
    return(SplitC(model))
  }
  return(model)
}


initRFlikelihood <- function(model, x, y = NULL, z = NULL, T=NULL, grid=NULL,
                             data, params, distances, dim, likelihood,
                             estimate_variance = NA,
                              Reg, ignore.trend=FALSE, ...) {

  if (!missing(likelihood)) ## composite likelihood etc
    stop("argument 'likelihood' is a future feature, not programmed yet")


  model <- PrepareModel2(model, params=params, ...)
  model <- ReplaceC(model); ## multivariates c() aufdroeseln
  
  model <- list("loglikelihood", model, data = data,
                estimate_variance=estimate_variance,
                betas_separate = FALSE, ignore_trend=ignore.trend)

  rfInit(model=model, x=x, y=y, z=z, T=T, grid=grid,
                 distances=distances, dim=dim, reg = Reg, dosimulate=FALSE)

}


RFlikelihood <- function(model, x, y = NULL, z = NULL, T=NULL, grid=NULL,
                         data, params, distances, dim, likelihood,
                         estimate_variance = NA,
                         ...) {
  relax <- is(model, "formula")
  RFoptOld <-
    if (missing(likelihood)) internal.rfoptions(..., RELAX=relax)
    else internal.rfoptions(likelihood=likelihood, ..., RELAX=relax)
  on.exit(RFoptions(LIST=RFoptOld[[1]]))
  RFopt <- RFoptOld[[2]]
  Reg <- RFopt$register$likelihood_register

  initRFlikelihood(model=model, x=x, y = y, z = z, T=T, grid=grid,
                   data=data, params=params,
                   distances=distances, dim=dim, likelihood=likelihood,
                   estimate_variance = estimate_variance,
                   Reg = Reg, ...)  

  likeli <- .Call(C_EvaluateModel, double(0),  Reg)
  info <- .Call(C_get_likeliinfo, Reg)
  globalvariance <- info$estimate_variance
  where <- 1 + globalvariance
  param <- likeli[-1:-where]
  if (length(param) > 0) names(param) <- info$betanames

  return(list(loglikelihood = likeli[1], likelihood = exp(likeli[1]),
              global.variance = if (globalvariance) likeli[where] else NULL,
              parameters = param
              )
         )
}


rfInit <- function(model, x, y = NULL, z = NULL, T=NULL, grid=FALSE,
                   distances, dim, reg, dosimulate=TRUE, old.seed=NULL) {
 
  stopifnot(xor(missing(x), #|| length(x)==0,
                missing(distances) || length(distances)==0))

  RFopt <- RFoptions() 
  if (!is.na(RFopt$basic$seed)) {
    allequal <- all.equal(old.seed, RFopt$basic$seed)
    allequal <- is.logical(allequal) && allequal
    if (dosimulate && RFopt$basic$printlevel >= PL_IMPORTANT &&
        (is.null(old.seed) || (!is.na(old.seed) && allequal)
         )
        ) {
      if (RFopt$internal$warn_seed) {
         RFoptions(internal.warn_seed = FALSE)
         txt <- "\nSet 'RFoptions(seed=NA)' to make the seed arbitrary."
      } else txt <- ""
      message("NOTE: simulation is performed with fixed random seed ",
               RFopt$basic$seed, ".", txt)
    }
    set.seed(RFopt$basic$seed)
  }
  ##  if (missing(x) || length(x) == 0) stop("'x' not given")

  new <- C_UnifyXT(x, y, z, T, grid=grid, distances=distances, dim=dim,
                 y.ok=!dosimulate)

##  Print(C_Init, as.integer(reg), model, new, NAOK=TRUE) 
  
  vdim <- .Call(C_Init, as.integer(reg), model, new, NAOK=TRUE) # ok
  
  if (is.null(old.seed)) return(vdim) else return(!is.na(RFopt$basic$seed))
}



rfdistr <- function(model, x, q, p, n, params, dim=1, ...) {
  ## note: * if x is a matrix, the points are expected to be given row-wise
  ##       * if not anisotropic Covariance expects distances as values of x
  ##
  ## here, in contrast to Covariance, nonstatCovMatrix needs only x

  RFoptOld <- internal.rfoptions(..., RELAX=is(model, "formula"))
  on.exit(RFoptions(LIST=RFoptOld[[1]]))
  RFopt <- RFoptOld[[2]]

  if (!missing(n) && n>10 && RFopt$internal$examples_reduced) {
    message("number of simulations reduced")
    n <- 10
  }
    
  model<- list("Distr", PrepareModel2(model, params=params, ...),
               dim=as.integer(dim));
  if (!missing(x)) {
    model$x <- if (is.matrix(x)) t(x) else x
  }
  if (!missing(q)) {
    model$q <- if (is.matrix(q)) t(q) else q
  }
  if (!missing(p)) {
    model$p <- if (is.matrix(p)) t(p) else p
  }
  if (!missing(n)) {
    model$n <- n
  }
  
  old.seed <- if (exists(".Random.seed")) .Random.seed else NULL
  if (rfInit(model=model, x=matrix(0, ncol=dim, nrow=1),
         y=NULL, z=NULL, T=NULL, grid=FALSE, reg = MODEL_DISTR,
         dosimulate=FALSE, old.seed=RFoptOld[[1]]$basic$seed) &&
      !is.null(old.seed))
    on.exit(.Random.seed <<- old.seed, add = TRUE)


  res <-  .Call(C_EvaluateModel, double(0), as.integer(MODEL_DISTR))

  if (RFoptOld[[2]]$general$returncall) attr(res, "call") <-
    as.character(deparse(match.call(call=sys.call(sys.parent()))))
  attr(res, "coord_system") <- c(orig=RFoptOld[[2]]$coords$coord_system,
                                 model=RFoptOld[[2]]$coords$new_coord_system)
  return(res)
}

RFdistr <- function(model, x, q, p, n, params, dim=1, ...) {
   rfdistr(model=model, x=x, q=q, p=p, n=n, params=params, dim=dim, ...)
}
RFddistr <- function(model, x, params, dim=1, ...) {
  if (hasArg("q") || hasArg("p") || hasArg("n")) stop("unknown argument(s)");
  rfdistr(model=model, x=x, params=params, dim=dim, ...)
}
RFpdistr <- function(model, q, params, dim=1, ...) {
  if (hasArg("x") || hasArg("p") || hasArg("n")) stop("unknown argument(s)");
  rfdistr(model=model, q=q, params=params, dim=dim, ...)
}
RFqdistr <- function(model, p, params, dim=1, ...){
  if (hasArg("x") || hasArg("q") || hasArg("n")) stop("unknown argument(s)");
  rfdistr(model=model, p=p, params=params, dim=dim, ...)
}
RFrdistr <- function(model, n, params, dim=1, ...) {
  if (hasArg("x") || hasArg("q") || hasArg("p")) stop("unknown argument(s)");
  rfdistr(model=model, n=n, params=params, dim=dim, ...)
}




rfeval <- function(model, x, y = NULL, z = NULL, T=NULL, grid=NULL,
                  params, distances, dim, ..., 
                  ##                  dim = ifelse(is.matrix(x), ncol(x), 1),
                  fctcall=c("Cov", "CovMatrix", "Variogram",
                    "Pseudovariogram", "Fctn"), reg=MODEL_USER) {

  ## note: * if x is a matrix, the points are expected to be given row-wise
  ##       * if not anisotropic Covariance expects distances as values of x
  ##
  ## here, in contrast to Covariance, nonstatCovMatrix needs only x



  RFoptOld <-internal.rfoptions(#xyz_notation=2*(length(y)!=0 && !is.matrix(y)),
				..., RELAX=is(model, "formula"))
  on.exit(RFoptions(LIST=RFoptOld[[1]]))
  
  fctcall <- match.arg(fctcall)

  if (fctcall != "CovMatrix" && !missing(distances) && !is.null(distances)) {
    if(missing(dim) || length(dim) != 1) {
      warning("'dim' not given or not of length 1, hence set to 1");
      dim <- 1
    }
    if (length(y) != 0 || length(z) != 0 || length(T) != 0 ||
        (!missing(grid) && length(grid) != 0))
      stop("if distances are given 'y', 'z', 'T', 'grid' may not be given")
    x <- (if (is.matrix(distances)) distances else
          cbind(distances, matrix(0, nrow=length(distances), ncol = dim - 1)))
    distances <- NULL
    dim <- NULL
    grid <- FALSE
  }
  
  p <- list(fctcall, PrepareModel2(model, params=params, ...))
  if (exists(".Random.seed")) {
    old.seed <- .Random.seed 
    on.exit(.Random.seed <<- old.seed, add = TRUE)
  }
  
  rfInit(model=p, x=x, y=y, z=z, T=T, grid=grid,
         distances=distances, dim=dim, reg = reg, dosimulate=FALSE,
         old.seed=RFoptOld[[1]]$basic$seed)
  res <- .Call(C_EvaluateModel, double(0), as.integer(reg))
  
  if (RFoptOld[[2]]$general$returncall) attr(res, "call") <-
    as.character(deparse(match.call(call=sys.call(sys.parent()))))
  attr(res, "coord_system") <- .Call(C_GetCoordSystem, reg,
              RFoptOld[[2]]$coords$coord_system,
              RFoptOld[[2]]$coords$new_coord_system)
   return(res)
}


RFcov <- function(model, x, y = NULL, z = NULL, T=NULL, grid,
		  params, distances, dim, ...,
		  data, bin=NULL, phi=NULL, theta = NULL,
		  deltaT = NULL, vdim=NULL) {
  if (hasArg("data")) {
    if (hasArg("model")) {
      stop("If both model and data are given, a least square fit will be performed soon")
    } else {
      rfempirical(x=x, y=y, z=z, T=T, data=data, grid=grid, bin=bin,
			   phi=phi, theta=theta, deltaT=deltaT, vdim=vdim,
			   method=METHOD_COVARIANCE, ...)
    }
  } else rfeval(model=model, x=x, y=y, z=z, T=T, grid=grid, params=params,
                distances=distances, dim=dim, ..., fctcall="Cov",
                reg=MODEL_COV)
}


RFcovmatrix <- function(model, x, y = NULL, z = NULL, T=NULL, grid,
                        params, distances, dim,  ...) {  
  rfeval(model=model, x=x, y=y, z=z, T=T, grid=grid, 
         distances=distances, dim=dim, params=params, ..., fctcall="CovMatrix",
         reg=MODEL_COVMATRIX)
}


extractVal <- function(L, val) {
  ans <- NULL
  for (i in 1:length(L)) {
    p <- L[[i]]
    if (is.list(p)) ans <- c(ans, extractVal(L, val))
    else if (is.numeric(p)) {
      for (j in 1:length(p)) {
        q <- pmatch(p[j], val)
        if (!is.na(q)) ans <- c(ans, q)
      }
    } ## else string -- does not matter 
  }
  return (ans)
}


RFvariogram <- function (model, x, y=NULL, z = NULL, T=NULL, grid,
			 params, distances, dim, ...,
			 data, bin=NULL, phi=NULL, theta = NULL,
			 deltaT = NULL, vdim=NULL){
  if (hasArg("data")) {
    if (hasArg("model")) {
      stop("If both model and data are given, a least square fit will be performed soon")
    } else {
      rfempirical(x=x, y=y, z=z, T=T, data=data, grid=grid, bin=bin,
			   phi=phi, theta=theta, deltaT=deltaT, vdim=vdim,
			   method=METHOD_VARIOGRAM, ...)
    }
  } else rfeval(model=model, x=x, y=y, z=z, T=T, grid=grid, params=params,
                distances=distances, dim=dim, ..., fctcall="Variogram",
                reg=MODEL_VARIOGRAM)
}

RFpseudovariogram <- function(model, x, y=NULL,  z = NULL, T=NULL, grid,
			      params, distances, dim, ...,
			      data, bin=NULL, phi=NULL, theta = NULL,
			      deltaT = NULL, vdim=NULL){
   if (hasArg("data")) {
    if (hasArg("model")) {
      stop("If model and data are given, a least square fit will be performed soon")
    } else {
      rfempirical(x=x, y=y, z=z, T=T, data=data, grid=grid, bin=bin,
			   phi=phi, theta=theta, deltaT=deltaT, vdim=vdim,
			   method=METHOD_PSEUDO, ...)
    }
   } else rfeval(model=model, x=x, y=y, z=z, T=T, grid=grid, params=params,
                distances=distances, dim=dim, ..., fctcall="Pseudovariogram",
                reg=MODEL_PSEUDO)
}

RFmadogram <- function(model, x, y=NULL,  z = NULL, T=NULL, grid, params,
			      distances, dim, ...,
			      data, bin=NULL, phi=NULL, theta = NULL,
			      deltaT = NULL, vdim=NULL){
   if (hasArg("data")) {
    if (hasArg("model")) {
      stop("If model and data are given, a least square fit will be performed soon")
    } else {
      rfempirical(x=x, y=y, z=z, T=T, data=data, grid=grid, bin=bin,
			   phi=phi, theta=theta, deltaT=deltaT, vdim=vdim,
			   method=METHOD_MADOGRAM, ...)
    }
  }  else stop("theoretical values for madograms cannot be calculated yet")
}

RFpseudomadogram <- function(model, x, y=NULL,  z = NULL, T=NULL, grid,
			      params, distances, dim, ...,
			      data, bin=NULL, phi=NULL, theta = NULL,
			      deltaT = NULL, vdim=NULL){
   if (hasArg("data")) {
    if (hasArg("model")) {
      stop("If model and data are given, a least square fit will be performed soon")
    } else {
      rfempirical(x=x, y=y, z=z, T=T, data=data, grid=grid, bin=bin,
		  phi=phi, theta=theta, deltaT=deltaT, vdim=vdim,
		  method=METHOD_PSEUDOMADOGRAM, ...)
    }
  }  else stop("theoretical values for madograms cannot be calculated yet")
}


RFfctn <- function(model, x, y=NULL,  z = NULL, T=NULL, grid,
                   params, distances, dim, ...) {
  rfeval(model=model, x=x, y=y, z=z, T=T, grid=grid, params=params,
                distances=distances, dim=dim, ..., fctcall="Fctn",
                reg=MODEL_FCTN )
}

RFcalc <- function(model, params, ...) {
  if (is.numeric(model)) return(model)
  rfeval(model=model, x=0, params=params, ...,
         coord_system="cartesian", new_coord_system="keep", spConform = FALSE,
         fctcall="Fctn",reg=MODEL_CALC)
}


######################################################################
######################################################################


rfDoSimulate <- function(n = 1, reg, spConform) {
  stopifnot(length(n) == 1, n>0, is.finite(n))
  RFopt <- RFoptions()
  if (missing(spConform)) spConform <- RFopt$general$spConform

  if (RFopt$gauss$paired && (n %% 2 != 0))
    stop("if paired, then n must be an even number")

  info <- RFgetModelInfo(RFopt$registers$register, level=3)

  vdim <- info$vdim
  total <- info$loc$totpts
  if (is.null(total) || total <= 0)
    stop("register ", RFopt$registers$register, " does not look initialized")

  result <- .Call(C_EvaluateModel, as.double(n), as.integer(reg)) #userdefined,
  if (!spConform) return(result)
  
  prep <- prepare4RFspDataFrame(info=info, RFopt=RFopt) 
  ## attributes(result)$varnames <- extractVarNames(model)  #prep$names$varnames
  
  res2 <- conventional2RFspDataFrame(result,
                                     coords=prep$coords,
                                     gridTopology=prep$gridTopology,
                                     n=n,
                                     vdim=info$vdim,
                                     T=info$loc$T,
                                     vdim_close_together
                                     =RFopt$general$vdim_close_together)
  return(res2)
}


    
RFsimulate <- function (model, x, y = NULL, z = NULL, T = NULL, grid=NULL,
                        distances, dim, data, given = NULL, err.model,
                        params, err.params, n = 1, ...) {

#  print("test"); 
 # Print(model, x, y,z,T, grid, given)
  
  if (!missing(model) && is(model, "RFfit"))
    stop("To continue with the output of 'RFfit' use 'predict' or give the components separately")
  mc <- as.character(deparse(match.call()))

### preparations #########################################################  
  if (!missing(distances) && length(distances)  > 0)
    RFoptOld <- internal.rfoptions(xyz_notation=length(y)!=0,
                                   expected_number_simu=n, ..., 
                                   general.spConform = FALSE,
                                   RELAX=!missing(model) && is(model,"formula"))
  else 
    RFoptOld <- internal.rfoptions(xyz_notation=length(y)!=0,
                                   expected_number_simu=n, ...,
                                   RELAX=!missing(model) && is(model,"formula"))
  
  on.exit(RFoptions(LIST=RFoptOld[[1]]))
  RFopt <- RFoptOld[[2]]

  if (n>2 && RFopt$internal$examples_reduced) {
    message("number of simulations reduced")
    n <- 2
  }
  
  cond.simu <- !missing(data) && !is.null(data) 
  reg <- RFopt$registers$register

  ### simulate from stored register ########################################
  mcall <- as.list(match.call(expand.dots=FALSE))
  if (length(mcall)==1 ||
      length(mcall)==2 && !is.null(mcall$n) ||
      length(mcall)==3 && !is.null(mcall$n) && "..." %in% names(mcall)) {
    if (cond.simu) {
      stop("repeated performance of conditional simulation not programmed yet")
    } else {
      ## userdefined <- GetParameterModelUser(PrepareModel2(model, params=params, ...))
      res <- rfDoSimulate(n=n, reg=reg, spConform=RFopt$general$spConform
                          #userdefined=userdefined
                          )
      if (RFopt$general$returncall) attr(res, "call") <- mc
      attr(res, "coord_system") <- .Call(C_GetCoordSystem, reg,
                                     RFopt$coords$coord_system,
                                     RFopt$coords$new_coord_system)
      return(res)
    }
  }
  
  ### preparations #########################################################
  stopifnot(!missing(model) && !is.null(model))

  model.orig <- model

  model <- PrepareModel2(model, params=params, ...)
  err.model <-
    if (missing(err.model)) NULL
    else PrepareModel2(err.model, params=err.params, ...)


  ### conditional simulation ###############################################
  if (cond.simu) {
    if (isSpObj(data)) data <- sp2RF(data)
    stopifnot(missing(distances) || is.null(distances))
    res <- switch(GetProcessType(model),
                  RPgauss = 
                  rfCondGauss(model=model.orig, x=x, y=y, z=z, T=T,
                              grid=grid, n=n, data=data, given=given,
                              err.model=err.model,
                              ## next line to make sure that this part
                              ## matches with predictGauss
                              predict_register = MODEL_PREDICT,
                              ...),
                  stop(GetProcessType(model),
                       ": conditional simulation of the process not programmed yet")
                  )
  } else { ## unconditional simulation ####
    if(!is.null(err.model))
      warning("error model is unused in unconditional simulation")

    if (exists(".Random.seed")) {
      old.seed <- .Random.seed 
      on.exit(.Random.seed <<- old.seed, add = TRUE)
    }
    rfInit(model=list("Simulate",
                         ##setseed=eval(parse(text="quote({set.seed(seed=seed); print(.Random.seed)})")),
                     setseed=eval(parse(text="quote(set.seed(seed=seed))")),
                  env=.GlobalEnv, model), x=x, y=y, z=z, T=T,
               grid=grid, distances=distances, dim=dim, reg=reg,
               old.seed=RFoptOld[[1]]$basic$seed)
    if (n < 1) return(NULL)
    res <- rfDoSimulate(n=n, reg=reg, spConform=FALSE)
  } # end of uncond simu


  ## output: RFsp   #################################
  if ((!cond.simu || (!missing(x) && length(x) != 0)) ## not imputing
      && RFopt$general$spConform) {
    info <- RFgetModelInfo(if (cond.simu) MODEL_PREDICT else reg, level=3)
    if (length(res) > 1e7) {
      message("Too big data set (more than 1e7 entries) to allow for 'spConform=TRUE'. So the data are returned as if 'spConform=FALSE'")
      return(res)
    }
    
    prep <- prepare4RFspDataFrame(info, RFopt, x, y, z, T, grid,
				  coordnames = attributes(res)$coordnames)
    attributes(res)$varnames <- attributes(res)$varnames
    
    
    res <- conventional2RFspDataFrame(data=res, coords=prep$coords,
                                      gridTopology=prep$gridTopology,
                                      n=n,
                                      vdim=info$vdim,
                                      T=info$loc$T,
                                      vdim_close_together=
                                      RFopt$general$vdim_close_together)
    if (is.raster(x)) {
      res <- raster::raster(res)
      raster::projection(res) <- raster::projection(x)
    }
  }
  
  if (RFopt$general$returncall) attr(res, "call") <- mc
  attr(res, "coord_system") <- .Call(C_GetCoordSystem,
                                     as.integer(reg),
                                     RFopt$coords$coord_system,
                                     RFopt$coords$new_coord_system)
  return(res)
}

					# RFreplace <- function(model, by) { }


RFoptions <- function(...) RandomFieldsUtils::RFoptions(...)