## 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.  



PrepareModel2 <- function(model, ..., params=NULL, allow.formulae=FALSE,
                          indices = NULL,
                          further.models = NULL,  arg,
                          x=NULL) {
  if (missing(model) || is.null(model)) stop("'model' must be given.")
  method <- "ml"
  if (is(model, "RF_fit")) model <- model[[method]]$model
  else if (is(model, "RFfit")) model <- model[method]
  
  transform <- NULL
  matrixandformulae <- FALSE
  if (missing(params) || is.null(params)) {
    allow.formulae <- FALSE
    m <- parseModel(model, ..., x=x)
  } else {
    Names <- names(params)
    tmpENV <- new.env(parent=.GlobalEnv)
    orig.params <- params
    formulae <- singleNA <- logical(length(params))
#    Print(formulae, params)
    for (i in 1:length(params)) {
      if (singleNA[i] <- !(formulae[i] <- is(params[[i]], "formula"))) {
        assign(Names[i], params[[i]], envir=tmpENV)
        singleNA[i] <- length(params[[i]]) == 1 && is.na(params[[i]])
        if(is.matrix(params[[i]]) && any(is.na(params[[i]])))
          matrixandformulae <- TRUE
        if (any(is.nan(params[[i]])))
          stop("only NAs are allowed to indicate parameters to be estimated")
      }
    }

    if (!is.null(further.models)) {
      if (is.character(further.models[[1]]))
        stop("Do not mix different styles of model definition.\n  See ?RFformula and ?RFformulaAdvanced for details.")
        ## Print(further.models, Names);
      FurtherModels <- "Bounds, initial values and parscales "
      for (n in names(further.models)) {
        ##Print(n, Names)
        idx <- which(n == Names)
        if (length(idx) > 0) {
#          Print(Names[idx])
          if (length(idx) != 1)
            stop("more than one formala given for '",  Names[idx], "'.")
          if (formulae[idx])
            stop(FurtherModels, "are not allowed for formulae, here '",
                 Names[idx], "'.")
          if (length(params[[n]]) != length(further.models[[n]])) {
            ## Print(n, params[[n]], sum(is.na(params[[n]])), length(further.models[[n]]), further.models[[n]])
            ## Print(get(n, envir=tmpENV))
            if (sum(is.na(params[[n]])) != length(further.models[[n]]))
              stop("Specification of '", n, "' (length) in '", arg,
                   "' does not match the model definition")
            params[[n]][is.na(params[[n]])] <- further.models[[n]]
          } else {
            #Print(params, further.models)            
            idx <- is.na(params[[n]])
            params[[n]][idx] <- further.models[[n]][idx]

            ##        Print(n, params[[n]], further.models[[n]], idx)
            
            idx <- !idx & !is.na(further.models[[n]]) &
              !is.na(further.models[[n]])

            #Print(idx)
            ok <-
              (if (arg == "lower")
               all(params[[n]][idx] >= further.models[[n]][idx])
              else if (arg == "upper")
               all(params[[n]][idx] <= further.models[[n]][idx])
              else 
               all(params[[n]][idx] == further.models[[n]][idx])
               )
            if (!ok) stop("Specification of '", n, "' in '", arg,
                          "' does not match the model definition.")
          }
        } else {
          ## Print(further.models, n, model); stop("unexpected unknown parameter name")
          assign(n, further.models[[n]], envir=tmpENV)
          params[[n]] <- further.models[[n]]
        }
      }
      if (length(params) != length(formulae))
        stop(FurtherModels, "may not define additional variables.\n  Here ",
             paste("'", names(params)[!(names(params) %in% Names)], "'",
                   sep="", collapse=","), ".")
      
      #Print(params, further.models, formulae)
    }
      
    
    ##    Print("Formulae", formulae, params)
    if (any(formulae)) {
      for (i in 1:length(params)) {
#        Print(i, params[[i]],try(eval(as.expression(params[[i]][[2]]),
                                        #                                 envir=tmpENV), silent = TRUE))
#        Print(i, formulae)
        if (formulae[i]) {
          #Print(params[[i]][[2]])
          params[[i]] <- try(eval(as.expression(params[[i]][[2]]),
                                  envir=tmpENV), silent = TRUE)
          if (is(params[[i]], "try-error")) params[[i]] <- NaN
          else if (any(idx <- is.na(params[[i]]))) params[[i]][idx] <- NaN
          else formulae[i] <- FALSE
        }
      }

                                        #
#      Print("after eval", params,allow.formulae , any(formulae), indices)
      
      if (allow.formulae && any(formulae)) {
        stopifnot(is.null(further.models))
        nas <- unlist(lapply(params,
                             function(x) {
                               if (is(x, "formula")) 1 else sum(is.na(x))
                             }
                             ))
        cumnas <- c(0, cumsum(nas))

#        Print(params, nas, cumnas);

        if (length(indices) == 0) {           
          ## prepare model to determine the indexing between R and C
           ## of the parameters that are NA
           values <- 1e306 +1e300 * (1:cumnas[length(cumnas)]) # in
               ## rf_interface.cc, includeparam, these values cannot appear
               ## as parameter values

          for (i in 1:length(params)) {
            params[[i]][is.na(params[[i]])] <-
              values[(cumnas[i] + 1) : cumnas[i + 1]]
          }
           
        } else  { ## further.models is set at least to list()
          if (matrixandformulae)
            warning("Note that some of the matrices are transposed internally.\nBe cautious in interpreting any result if matrices with NAs are used in formulae.")
          if (length(indices)<=1)
            stop("Got a wrong index set. Please contact maintainer")
          rev.idx <- order(indices)
          ftext <- character(sum(formulae)) # for each formula a text

          ff <- orig.params[formulae]
          f.exp <- vector("list", length(ff))
          for (i in 1:length(ff)) {
            f.exp[[i]] <- as.expression(ff[[i]][[2]])
            ftext[i] <- as.character(f.exp[[i]])
          }

          ## formula: 1 NaN
          ## constant: neither NA nor NaN
          ## variable: NAs
          f <- unlist(lapply(params, ## unknown referring to formulae
                             function(x) {
                               nans <- sum(is.nan(x))
                               nas <- sum(is.na(x) & !is.nan(x))
                               if (nans * nas > 0) stop("NAs and forumulae may not be given at the same time within one parameter.")
                               #Print(x, nans, nas)
                               rep(nans > 0, nas + nans)
                             }))

          where.from <- sapply(params, function(x) {
            idx <- is.na(x)
            if (all(idx)) ""
            else paste0("[c(", paste(which(idx), collapse=","), ")]")
            })[formulae]
            
     
          ##          Print(rev.idx, f, params, ftext, ff); kkkkkk
          small.idx <- rev.idx[which(!f)]
          small2big <- sort(small.idx)
          which.f <- which(f)
          formulae2big <- paste(rev.idx[which.f], collapse=", ")
          

       #   formulae | !f

          
          s <- unlist(lapply(params,
                             function(x) {
                               if (any(is.nan(x))) NULL else
                               rep(length(x) == 1, sum(is.na(x)))
                             }))
          known <- sapply(params, function(x) !any(is.na(x)))


          ## formulae may depend on each other; find the ordering in which
          ## the formulae should be evaluated
          names.formulae <- Names[formulae] ## ftext, formuale2big
          len <- length(names.formulae)
          ordering <- rep(0, len)
          zaehler <- 0
          vectors <- FALSE
          while (zaehler < len) {
            ok <- FALSE ## each round, at least one other formula should
            ##             be evaluable
            for (i in 1:len) {
              if (ordering[i] != 0) next
              #Print(vectors, names.formulae[i])
              if (!vectors && length(params[[names.formulae[i]]]) > 1) next
              ## cat(names.formulae[i], "= ") ; print(f.exp[[i]])
              if (!is(f0 <- try(eval(f.exp[[i]], envir=tmpENV), silent=TRUE),
                      "try-error")) {
                ## only scalars may be used in other formulae!
                if (!vectors) assign(names.formulae[i], f0, envir=tmpENV) 
                ordering[i] <- zaehler <- zaehler + 1
                ok <- TRUE ## one formula was evaluable
              }
            }
            if (!ok) {
              if (vectors)
                stop("Formula", if (sum(ordering == 0) > 1) "e", " for ",
                     paste("'", names.formulae[ordering == 0], "'",
                           sep="", collapse=","),
                     " not evaluable.\n  Only scalars may be used within a formula and all of them must be defined.")
              else vectors <- TRUE
            }
          }
          ordering <- order(ordering)           
          
          if ("..bca.." %in% Names)
            stop("'..bca..' is not allowed as parameter name")
          n.output <- length(small2big) + length(which.f)
                                        #          Print(small2big, which.f)

          text <- paste("function(variab) {\n",
                        paste("\t", Names[singleNA], "<- variab[",
                              rank(small.idx)[s], "]",
                              collapse="\n"),
                        "\n\t",paste(Names[known], orig.params[known], sep="<-",
                                     collapse="\n\t"),
                        "\n\t",paste(names.formulae[ordering], "<-",
                                     ftext[ordering], collapse="\n\t"),
                        "\n\t..bca.. <- double(", n.output, ")\n",
                        "\t..bca..[ c(", paste(small2big, collapse=", "),
                        ") ] <- variab", "\n",
                        "trafos <- c(", paste(names.formulae, where.from,
                                              sep="",collapse=", "), ")\n",
                        "trafos[is.na(trafos)] <- NaN \n",
                        "\t..bca..[ c(", formulae2big, ") ] <- trafos\n",
                         "\treturn(..bca..)\n}")
          ##          Print(params, text, formulae, Names[formulae]);pppp
          fctn <- eval(parse(text=text), envir=NULL)
          environment(fctn) <- baseenv()


          ##Print(f[indices]); pppp
          
          transform <- list(!f[indices], fctn)

          if (FALSE) # Bsp in RMdeclare funktioniert sonst nicht
            if (any(sort(c(small.idx, which.f)) != 1:n.output)) {
              # Print(small.idx, which.f, formulae2big, n.output, s, known, text, transform) 
            stop("error in creation of 'transform'. Please inform author of the package")
          }
              #          Print(rank(small.idx)[s], Names[singleNA], model, indices,                f, rev.idx, text, Names, singleNA, small2big, formulae2big, transform, transform[[2]](4:5), small.idx, rank(small.idx)); 
          #pppp
        } # ind
      } # if alllow.formulae
    } # if anyformulae


#    Print("parsemodel",further.models, c(list(model=model, ..x..=x), params));
     m <- do.call("parseModel", c(list(model=model, x=x), params))  
    ##
#    Print(m, c(list(model=model, x=x), params)); 
  } # params !null
        
  if (notplus <- !(m[[1]] %in% RM_PLUS)) m <- list(SYMBOL_PLUS, m)     
  if (notplus) m <- m[[2]]
  class(m) <- "RM_model"

  if (allow.formulae) {
    if (is.null(indices)) ## different from length(indices) == 0 !!
      ## latter happened if not formula and called second time from Unifydata
      ## null happens if called first from UnifyData
      return (if (any(formulae)) list(model=m, values=values) else m)
    else {
      #Print(transform, text);
      return(list(model=m, transform = transform))
    }
  }

 # Print("ende prepare", m)
   
  return(m)
}


PrepareModel <-  function(model, param, trend=NULL, 
                          nugget.remove=TRUE, method=NULL) {
  ## any of the users model definition (standard, nested, list) for the
  ## covariance function is transformed into a standard format, used
  ## especially in the c programs
  ##
  ## overwrites in some situation the simulation method for nugget.
  ## allows trend to be NA (or any other non finite value  -- is not checked!)
  ## trend has not been implemented yet!

  if (is(model, CLASS_CLIST))
    stop("models of class CLASS_CLIST cannot be combined with obsolete RandomFields functions")

  if (!is.null(method)) stop("to give method in suPrepareModel is obsolete")
  
  if (!is.null(trend))      
     if (!is.numeric(trend) || length(trend)!=1)
        stop("in the obsolete setting, only constant mean can used")
 
  if (is.list(model) && is.character(model[[1]]) &&
      (is.null(names(model)) || names(model)[[1]]=="")) {
     if (!missing(param) && !is.null(param))
        stop("param cannot be given in the extended definition")
     if (is.null(trend)) return(model)
     trend <- list(RM_TREND[2], mean=trend)
     if (model[[1]] %in% RM_PLUS) return(c(model, list(trend)))
     else return(list(SYMBOL_PLUS, model, trend))
  }
    
  printlevel <- RFoptions()$basic$printlevel
  STOP <- function(txt) {
    if (printlevel>=PL_ERRORS) {
      cat("model: ")
      if (!missing.model) Print(model) else cat(" missing.\n") #
      cat("param: ")
      if (!missing.param) Print(param) else cat(" missing.\n") #
      cat("trend: ")
      Print(trend) # 
     }
    stop("(in PrepareModel) ", txt, call.=FALSE)
  }
   
  transform <- function(model) {
    if (!is.list(model)) {
      STOP("some elements of the model definition are not lists")
    }
    m <- list(DOLLAR[1], var=model$v)
    lm <- length(model) - 3 # var, scale/aniso, name
    if (!is.null(model$a)) m$aniso <- model$a else m$scale <- model$scale
##    model <- c(model, if (!is.null(model$a))
##               list(aniso=model$a) else list(scale=model$s)) ## ???

    
    if (!is.na(p <- pmatch("meth", names(model), duplicates.ok=TRUE))) {
      if (printlevel>=PL_ERRORS)  Print(p, model) #
      stop("method cannot be given with the model anymore. It must be given as a parameter to the function. See 'RFoptions' and 'RFsimulate'")
     }
  
    if (!is.null(model$me))
      stop("'mean' seems to be given within the inner model definitions"); 
    if (!is.character(model$m)) {
       stop("'model' was not given extacly once each odd number of list entries or additional unused list elements are given.")
    }
    m1 <- list(model$m)
    if (!is.null(model$k)) {
      lm <- lm - 1
      if (length(model$k) != 0)
        for (i in 1:length(model$k)) {
          eval(parse(text=paste("m1$k", i, " <- model$k[", i, "]", sep="")))
      }
    }
    if (lm != 0) {
      if (printlevel>=PL_ERRORS) Print(lm, model) #
      stop("some parameters do not fit")
    }
    m <- c(m, list(m1))

    return(m)
    
  } # end transform

  op.list <- c(SYMBOL_PLUS, SYMBOL_MULT)  ## if others use complex list definition !
  missing.model <- missing(model)
  missing.param <- missing(param) || is.null(param)

  if (missing.param && is.null(model$param)) { ## full model
    if (RFoptions()$internal$warn_oldstyle)
      warning("the sequential list format is depreciated.")
    if (missing.model || (length(model)==0)) model <- list()
    else if (!is.list(model))
      STOP("if param is missing, model must be a list of lists (or a list in the extended notation)")
    if (is.null(trend) + is.null(model$mean) + is.null(model$trend)<2)
      STOP("trend/mean is given twice")
    if (!is.null(model$mean)) trend <- model$mean else
    if (!is.null(model$trend)) trend <- model$trend else trend <- NULL

    model$trend <- model$mean <- NULL
    ## the definition might be given at a deeper level as element
    ## $model of the list:
    if (is.list(model$model)) {
      if (!is.list(model$model[[1]]))
        STOP("if param is missing, the model$model must be a list of lists")
      model <- model$model
    }
    if (length(model)==0) { ## deterministic      
      return(if (is.null(trend)) NULL else list(RM_TREND[2], mean=trend))
    }
    if (length(model) %% 2 !=1) STOP("list for model definition should be odd")
    if (length(model)==1)
      return(if (is.null(trend) ||
                 is.numeric(trend) && length(trend)==1 && !is.na(trend)&&trend==0)
             transform(model[[1]]) 
             else list(SYMBOL_PLUS, transform(model[[1]]),
                       list(RM_TREND[2], mean=trend)));

    op <- pmatch(c(model[seq(2, length(model), 2)], recursive=TRUE),
                 op.list, duplicates.ok=TRUE) - 1
    if (!all(is.finite(op))) STOP("operators are not all allowed; see the extended list definition for extensions")
    model <- model[seq(1, length(model), 2)]

    plus <- which(op==0)
    if (length(plus) == 0) {
      m <- list("*", lapply(model, transform))
    } else {
      plus <- c(0, plus, length(op)+1)
      m <- list(SYMBOL_PLUS)
      for (i in 1:(length(plus) - 1)) {
        m[[i+1]] <-
          if (plus[i] + 1 == plus[i+1]) transform(model[[plus[i] + 1]])
          else list(SYMBOL_MULT,
                    lapply(model[(plus[i] + 1) : plus[i+1]], transform))
      }
    }
   model <- m
  } else { ## standard definition or nested model
    if (missing.param) { ## a simple list of the model and the
      ##                    parameters is also possible
      if (is.null(param <- model$p)) STOP("is.null(model$param)")
      stopifnot(is.null(trend) || is.null(model$trend))
      if (is.null(trend)) trend <- model$trend
      if (!is.null(model$mean)) {
        if (!is.null(trend)) STOP("mean and trend given twice")
        trend <- model$mean
      }
      model <- model$model
    }
    stopifnot(is.character(model), length(model)==1)
    if (is.matrix(param)) { ## nested
      if (nrow(param) == 1)
        return(PrepareModel(model=model, param=c(param[1], 0, param[-1]),
                            trend=trend))
      name <- model
      model <- list(SYMBOL_PLUS)#, method=method)
      for (i in 1:nrow(param)) {
        model <- c(model,
                   if (is.na(param[i, 2]) || param[i, 2] != 0)
                   list(list(DOLLAR[1], var=param[i, 1], scale=param[i, 2],
                             if (ncol(param) >2) list(name, k=param[i,-1:-2])
                             else list(name)))
                   else list(list(DOLLAR[1], var=param[i,1],
                                  list(RM_NUGGET[2]))))
      }
    } else if (is.vector(param)) {  ## standard, simple way
      ## falls trend gegeben, dann ist param um 1 Komponente gekuerzt
      if (is.null(trend)) {
        trend <- param[1]
        param <- param[-1]
      } else message("It is assumed that no mean is given so that the first component of param is the variance")
      if (model == RM_NUGGET[2]) {
        model <- transform(list(model=model, var=sum(param[1:2]), scale=1)) 
      } else {
        if  (length(param) > 3)
          model <- transform(list(model=model, var=param[1], scale=param[3],
                                  k=param[-1:-3]))
        else 
          model <- transform(list(model=model, var=param[1], scale=param[3]))
        if (is.na(param[2]) || param[2] != 0 || !nugget.remove) {# nugget
          model <- list(SYMBOL_PLUS,
                        model,
                        transform(list(model=RM_NUGGET[2], var=param[2], scale=1)))
        }
        ## if (!is.null(method)) model <- c(model, method=method) ## doppelt
      }
    } else stop("unknown format")  # end nested/standard definition
  }

  return(if (is.null(trend) ||
             is.numeric(trend) && length(trend)==1 &&  !is.na(trend) &&trend==0)
         return(model)
         else if (model[[1]] %in% RM_PLUS)
                 c(model, list(list(RM_TREND[2], mean=trend)))
         else list(SYMBOL_PLUS, model, list(RM_TREND[2], mean=trend)))
}


seq2grid <- function(x, name, grid, warn_ambiguous, gridtolerance) {
  xx <- matrix(nrow=3, ncol=length(x))
  step0 <- rep(FALSE, length(x))
  gridnotgiven <- missing(grid) || length(grid) == 0
  
  for (i in 1:length(x)) {
    if (length(x[[i]]) == 1) {
      xx[,i] <- c(x[[i]], 0, 1)
      next
    }
    step <- diff(x[[i]])
    if (step[1] == 0.0) {
      
      ok <- step0[i] <- all(step == 0.0)      
    } else {
      ok <- max(abs(step / step[1] - 1.0)) <= gridtolerance
    }

    if (!ok) {
      if (gridnotgiven) return(FALSE)
      if (!TRUE)
        Print(i, x[[i]][1:min(100, length(x[[i]]))], #
              step[1:min(100,length(step))],
              range(diff(step[1:min(100,length(step))])))
      stop("Different grid distances detected, but the grid must ",
           "have equal distances in each direction -- if sure that ",
           "it is a grid, increase the value of 'gridtolerance' which equals ",
           gridtolerance,".\n")
    }

    xx[,i] <- c(x[[i]][1], step[1], if (step0[i]) 1 else length(x[[i]]))
  }

  if (FALSE && gridnotgiven && warn_ambiguous && length(x) > 1) {
    RFoptions(internal.warn_ambiguous = FALSE)
    message("Ambiguous interpretation of coordinates. Better give 'grid=TRUE' explicitly. (This message appears only once per session.)")
  }

  if (any(step0)) {
    if (all(step0)) {
      if (gridnotgiven) return(FALSE)
      else stop("Within a grid, the coordinates must be distinguishable")
    } else {
      if (gridnotgiven && warn_ambiguous) {
        RFoptions(internal.warn_ambiguous = FALSE)
        warning("Interpretation as degenerated grid. Better give 'grid' explicitely. (This warning appears only once per session.)")
      }
    }
  }

  return(xx)
}

UnifyXT <- function(x, y=NULL, z=NULL, T=NULL, grid, distances=NULL,
                    dim=NULL, # == spatialdim!
                    length.data,
                    y.ok = FALSE, 
                    printlevel = RFopt$basic$printlevel,
                    allow_duplicated =
                      RFopt$internal$allow_duplicated_locations ){

  
  ## do not pass anything on "..." ! --- only used for internal calls
  ## when lists are re-passed

  ## converts the given coordinates into standard formats
  ## (one for arbitrarily given locations and one for grid points)
  #print("UnifyXT in convert.R")#Berreth
  
  RFopt <- RFoptions()
  
  if (!missing(x)) {
    if (is(x, "UnifyXT")) return(x)  
    if (is.list(x) && !is.data.frame(x))  {
      if (!is.list(x[[1]]))
        return(do.call("UnifyXT", c(x, list(printlevel=printlevel))))
      L <- list()
      for (i in 1:length(x)) {
        L[[i]] <-
          if (is(x[[i]], "UnifyXT")) x[[i]]
          else do.call("UnifyXT", c(x[[i]], list(printlevel=printlevel)))
      }
      if (length(x) > 1) {    
        if (!all(diff(sapply(L, function(x) x$has.time.comp)) == 0) ||
          !all(diff(sapply(L, function(x) x$spatialdim)) == 0))
          stop("all sets must have the same dimension")
        if (!all(diff(sapply(L, function(x) x$dist.given)) == 0))
        stop("either all the sets must be based on distances or none")
      }
      class(L) <- "UnifyXT"
      return(L)
    }
  }

  
  curunits <- RFopt$coords$coordunits
  newunits <-  RFopt$coords$new_coordunits
  coord_system <-  RFopt$coords$coord_system
  new_coord_system <-  RFopt$coords$new_coord_system
  ex.red <- RFopt$internal$examples_reduced
  
  if (!missing(distances) && !is.null(distances)) { ## length==0 OK!
    stopifnot(is.matrix(distances) || is(distances, "dist") ||
                                       is.vector(distances),
              (!missing(dim) && !is.null(dim)),
              (missing(grid) || length(grid) == 0),
              missing(x) || is.null(x),
              length(y)==0,
              length(z)==0,
              length(T)==0)
    
    if (coord_system != new_coord_system && new_coord_system != "keep")      
      stop("coordinate systems differ")
    
    if (is.list(distances)) {
      L <- list()
      for (i in 1:length(distances))
        L[[i]] <- do.call("UnifyXT", list(distances=distances[[i]], dim=dim,
                                          printlevel=printlevel))
       class(L) <- "UnifyXT"
      return(L)
    }
        
    if (is(distances, "dist")) {
      x <- as.vector(distances)
      len <- length(distances)
    } else if (is.matrix(distances) || is.vector(distances)) {
      if (is.matrix(distances)) {        
        len <- nrow(distances)
        if (is.null(dim)) dim = ncol(distances)
        else if (dim != ncol(distances))
          stop("matrix of distances does not fit the given dimension")
      } else {
        len <- length(distances)
        if (is.null(dim))
          stop("dim is not given although 'distances' are used")
      }
      x <- distances
    } else {
      stop("'distances' not of required format.")
    }

    if (ex.red && len > ex.red^2 / 2) {
      LEN <- as.integer(ex.red)
      len <- as.integer(LEN * (LEN - 1) / 2)
      x <- if (is.matrix(x)) x[1:len ,] else x[1:len]
    } else {
      LEN <- as.integer(1e-9 + 0.5 * (1 + sqrt(1 + 8 * len)))
      if (LEN * (LEN-1) / 2 != len) LEN <- NaN
    }

    ## keep exactly the sequence up to 'distances'
    if (storage.mode(x) != "double") storage.mode(x) <- "double"
    L <- list(x = as.matrix(x), #0
              y = double(0),   #1
                T= double(0),  #2
              grid = FALSE, #3
              spatialdim=as.integer(dim),#4
              has.time.comp=FALSE, #5
              dist.given = TRUE, #6
              restotal = LEN, ## number of points
              l = LEN, ## ?? physical length??
              coordunits = curunits,
              new_coordunits = newunits
              )
    class(L) <- "UnifyXT"
    return(L)
  }

  stopifnot(!missing(x))
  if (is(x, "RFsp") || isSpObj(x)) {
    return(UnifyXT(x=coordinates(x), y=y, z=z, T=T, grid=grid,
                   distances=distances, dim=dim, length.data=length.data,
                   y.ok=y.ok, printlevel=printlevel))
  }    
  if (is.raster(x)) x <- as(x, 'GridTopology')
 
  if ((missing(grid) || length(grid) == 0) && !missing(length.data)) {
    new <-  try(UnifyXT(x=x, y=y, z=z, T=T, grid=TRUE, distances=distances,
                        dim=if (!missing(dim)) dim,
                        length.data = length.data, y.ok =y.ok,
                        printlevel = printlevel), silent=TRUE)
    if (grid <- !is(new, "try-error")) {
      ratio <- length.data / new$restotal

      if (grid <- ratio == as.integer(ratio)) {
        if (printlevel>=PL_IMPORTANT && new$spatialdim > 1)
          message("Grid detected. If it is not a grid, set grid=FALSE.\n")
      }
    }
    return(if (grid) new else {
      UnifyXT(x, y, z, T, grid=FALSE, distances,
              if (!missing(distances) && length(distances) > 0) dim=1,
              length.data = length.data,
              printlevel = printlevel) }
           )
  } # if (missing(grid) && !missing(length.data))


  gridtriple <- FALSE

  if (is.GridTopology <- is(x, "GridTopology")){
    x <- rbind(x@cellcentre.offset,
               x@cellsize,
               x@cells.dim)
    if ((missing(grid) || length(grid) == 0)) grid <- TRUE else stopifnot(grid)
    gridtriple <- TRUE
  }
  ##else {
  ##  is.GridTopology <- FALSE
  ##}

  
  if (is.data.frame(x)) {
    if (ncol(x)==1) x <- as.vector(x) else x <- as.matrix(x)
  }
  
  stopifnot(length(x) != 0)
#  stopifnot(all(unlist(lapply(as.list(x), FUN=function(li) is.numeric(li))))) ## wann benoetigt???

  stopifnot(is.numeric(x))# um RFsimulte(model, data) statt data=data abzufangen
  
  
#  stopifnot(all(is.finite(x)), all(is.finite(y)), all(is.finite(z))) ; s.u. unlist
 
     
  if (is.matrix(x)) {
    if (!is.numeric(x)) stop("x is not numeric.")
    if (length(z)!=0) stop("If x is a matrix, then z may not be given")
    if (length(y)!=0) {
      if (!y.ok) stop("If x is a matrix, then y may not be given")
      if (length(T)!=0)
        stop("If x is a matrix and y is given, then T may not be given")
      if (!is.matrix(y) || ncol(y) != ncol(x) ||
          nrow(x)==3 && nrow(y)!=3 && ((missing(grid) || length(grid) == 0) ||
                                grid))
        stop("y does not match x (it must be a matrix)")
    }

    if (coord_system == COORD_SYS_NAMES[coord_auto + 1] && ncol(x) >= 2
        && ncol(x) <= 3 && !is.null(n <- dimnames(x)[[2]])) {
      if (any(idx <- earth_coordinate_names(n))) {
        if (length(idx) == 2 && !all(idx == 1:2))
          stop("earth coordinates not in order longitude/latitude")
        cur <- curunits[1]
        newunits <- RFopt$coords$new_coordunits
        curunits <- RFopt$coords$coordunits
        curunits[1:2] <- COORD_NAMES_EARTH[1:2]
        if (newunits[1] == "") newunits[1] <-  UNITS_NAMES[units_km + 1]
        newunits[2:3] <- newunits[1]                
        if (RFopt$internal$warn_coordinates) {
          message("\n\nNOTE: current units are ",
                  if (cur=="") "not given and" else paste("'", cur, "', but"),
                  " earth coordinates detected:\n",
                  "earth coordinates will be transformed into units of '",
                  newunits[1],
                  "'.\nIn particular, the values of all scale parameters of ",
                  "any model defined\nin R^3 (currently all models!) are ",
                  "understood in units of '", newunits[1],
                  "'.\nChange options 'coord_system' and/or 'units' if ",
                  "necessary.\nNOTE: angles in R.cos, R.sin, R.tan, RMangle ",
                  "are now expected in DEGREE and\n",
                  "R.acos, R.asin, R.atan, R.atan2 return results in degree.\n",
                 "(This message appears in long form only once per session.)\n")
        } else message(" earth coordinates detected that will be transformed ",
                     " into units of ",  newunits[1], " .")
        coord_system <- COORD_SYS_NAMES[earth + 1]
        RFoptions(coords.coord_system = coord_system,
                  coords.coordunits = curunits,
                  coords.new_coordunits = newunits,
                  internal.warn_coordinates=FALSE)
       } else {
         RFoptions(coords.coord_system =  COORD_SYS_NAMES[cartesian + 1])
      }
    }    
   
    spatialdim <- ncol(x)
    len <- nrow(x)
    if (spatialdim==1 && len != 3 && (missing(grid) || length(grid) == 0)) {
      if (length(x) <= 2) grid <- TRUE
      else {
        dx <- diff(x)
        grid <- max(abs(diff(dx))) < dx[1] * RFopt$general$gridtolerance
      }
    } # else {

    if ((missing(grid) || length(grid) == 0) &&
        any(apply(x, 2, function(z) (length(z) <= 2) || max(abs(diff(diff(z))))
                  > RFopt$general$gridtolerance))) {
      grid <- FALSE
    }

    if ((missing(grid) || length(grid) == 0) || !is.logical(grid)) {
      grid <- TRUE
      if (spatialdim > 1 && RFopt$internal$warn_ambiguous) {
        RFoptions(internal.warn_ambiguous = FALSE)
        warning("Ambiguous interpretation of the coordinates. Better give the logical parameter 'grid=TRUE' explicitely. (This warning appears only once per session.)")
      }
    }

    if (grid && !is.GridTopology) {
      if (gridtriple <- len==3) {
        if (printlevel >= PL_SUBIMPORTANT && RFopt$internal$warn_oldstyle) {
#          print(x); dsddfdsf
          message("x was interpreted as a gridtriple; the new gridtriple notation is:\n  1st row of x is interpreted as starting values of sequences,\n  2nd row as step,\n  3rd row as number of points (i.e. length),\n  in each of the ", ncol(x), " directions.")
        } 
      } else len <- rep(len, times=spatialdim)   # Alex 8.10.2011
    }

    if (grid && !gridtriple) {
      ## list with columns as list elements -- easier way to
      ## do it??
      x <- lapply(apply(x, 2, list), function(r) r[[1]])
      if (length(y) != 0) y <- lapply(apply(y, 2, list), function(r) r[[1]])
    }
  } else { ## x, y, z given separately
    if (length(y)==0 && length(z)!=0) stop("y is not given, but z")
    xyzT <- list(x=if (!missing(x)) x, y=y, z=z, T=T)
    for (i in 1:4) {
      if (!is.null(xyzT[[i]]) && !is.numeric(xyzT[[i]])) {
        if (printlevel>PL_IMPORTANT) 
          message(names(xyzT)[i],
                  " not being numeric it is converted to numeric")
        assign(names(xyzT)[i], as.numeric(xyzT[[i]]))
      }
    }
    remove(xyzT)
    spatialdim <- 1 + (length(y)!=0) + (length(z)!=0)
    if (spatialdim==1 && ((missing(grid) || length(grid) == 0) || !grid)) {
      ## ueberschreibt Einstellung des Nutzers im Falle d=1
      if (length(x) <= 2) newgrid <- TRUE
      else {
        dx <- diff(x)
        newgrid <- max(abs(diff(dx))) < dx[1] * RFopt$general$gridtolerance
      }
      if ((missing(grid) || length(grid) == 0)) grid <- newgrid
      else if (xor(newgrid, grid) && RFopt$internal$warn_on_grid) {
        RFoptions(internal.warn_on_grid = FALSE)
        message("coordinates", if (grid) " do not",
                " seem to be on a grid, but grid = ", grid)
      }
    }
    len <- c(length(x), length(y), length(z))[1:spatialdim]
    
    if (!(missing(grid) || length(grid) == 0) && !grid) { ## sicher nicht grid, ansonsten ausprobieren
      if (any(diff(len) != 0)) stop("some of x, y, z differ in length")
      x <- cbind(x, y, z)
      ## make a matrix out of the list
      len <- len[1]
    } else {
      if ((missing(grid) || length(grid) == 0) && any(len != len[1]))
        grid <- TRUE
      x <- list(x, y, z)[1:spatialdim]
    }
    y <- z <- NULL ## wichtig dass y = NULL ist, da unten die Abfrage
  }  ## end of x, y, z given separately 
  
  if (!all(is.finite(unlist(x)))) {
    stop("coordinates are not all finite")
  }

  if ((missing(grid) || length(grid) == 0) || grid) {
    if (gridtriple) {
      if (len != 3)
        stop("In case of simulating a grid with option gridtriple, exactly 3 numbers are needed for each direction")
      lr <- x[3,] # apply(x, 2, function(r) length(seq(r[1], r[2], r[3])))
      ##x[2,] <- x[1,] + (lr - 0.999) * x[3,] ## since own algorithm recalculates
      ##                               the sequence, this makes sure that
      ##                               I will certainly get the result of seq
      ##                               altough numerical errors may occurs
      restotal <- prod(x[3, ])
      if (length(y)!=0 && !all(y[3,] == x[3,]))
        stop("the grids of x and y do not match ")        
    } else {     
      xx <- seq2grid(x, "x",  grid,
                     RFopt$internal$warn_ambiguous, RFopt$general$gridtolerance)
     if (length(y)!=0) {
        yy <- seq2grid(y, "y", grid,
                       RFopt$internal$warn_ambiguous,
                       RFopt$general$gridtolerance)
        if (xor(is.logical(xx), is.logical(yy)) ||
            (!is.logical(xx) && !all(yy[3,] == xx[3,])))
          stop("the grids for x and y do not match")      
      }
      if (missing(grid) || length(grid) == 0) grid <- !is.logical(xx)       
      if (grid) {
        x <- xx
        if (length(y) != 0) y <- yy
        restotal <- prod(len)
        len <- 3
      } else {
        x <- sapply(x, function(z) z)
        if (length(y) != 0) y <- sapply(y, function(z) z)
      }
    }
    if (grid && any(x[3, ] <= 0))
      stop(paste("step must be postive. Got as steps",
                 paste(x[3,], collapse=",")))
    ##if (len == 1) stop("Use grid=FALSE if only a single point is simulated")
  }
 
  if (!grid) {
    restotal <- nrow(x)
    if (length(y)==0 && !allow_duplicated) {
      if (restotal < 200 && any(as.double(dist(x)) == 0)) {
        d <- as.matrix(dist(x))
        diag(d) <- 1
        idx <-  which(as.matrix(d) ==0)
        if (printlevel>PL_ERRORS)
          Print(x, dim(d), idx , cbind( 1 + ((idx-1)%% nrow(d)), #
                                       1 + as.integer((idx - 1)  / nrow(d))) ) 
        stop("Some locations are not distinguishable. If this is intentional, see ?RPnugget and ?RFoptions 'allow_duplicated_locations' to deal with repeated measurements.")
      }
      ## fuer hoehere Werte con total ist ueberpruefung nicht mehr praktikabel
    }
  }
 
  if (coord_system == "earth") {
    # if (ncol(x) > 4) stop("earth coordinates have maximal 3 components")
    opt <- RFoptions()$coords ## muss nochmals neu sein
    global.units <- opt$new_coordunits[1]
    if (global.units[1] == "") global.units <- "km"
   
    Raumdim <- ncol(x) #if (grid) ncol(x) else
    new_is_cartesian <- new_coord_system %in% CARTESIAN_SYS_NAMES
    if (new_is_cartesian) {
      if (sum(idx <- is.na(opt$zenit))) {
        zenit <- (if (grid) x[1, 1:2] + x[2, 1:2] * (x[3, 1:2] - 1)
                  else if (opt$zenit[!idx] == 1) colMeans(x[, 1:2])
                  else if (opt$zenit[!idx] == Inf) colMeans(apply(x[, 1:2], 2,
                                                                  range))
                  else stop("unknown value of zenit"))
         RFoptions(zenit = zenit)
      }

      code <- switch(new_coord_system,
                     "cartesian" = CARTESIAN_COORD,
                     "gnomonic" = GNOMONIC_PROJ,
                     "orthographic" = ORTHOGRAPHIC_PROJ,
                     stop("unknown projection method")
                     )
      message("New coordinate system: ", new_coord_system, ".\n")
      x <- RFfctn(RMtrafo(new=code), x, grid=grid, 
                   coords.new_coordunits=global.units,
                   coords.new_coord_system = "keep")
      
       if (length(y) != 0)         
         y <- RFfctn(RMtrafo(new=code), y, grid=grid, 
                   coords.new_coordunits=global.units,
                   coords.new_coord_system = "keep")
     
      if (new_coord_system == "cartesian") {
        Raumdim <- max(3, Raumdim)
        spatialdim <- Raumdim
      }
      dim(x) <- c(length(x) /Raumdim, Raumdim)
      #x <- t(x)

      ## never try to set the following lines outside the 'if (new_coord_system'
      ## as in case of ..="keep" none of the following lines should be set
      RFoptions(coords.coord_system = 
                if (new_is_cartesian) "cartesian" else new_coord_system)
      grid <- FALSE
    } else if (!(new_coord_system %in% c("keep", "sphere", "earth"))) {
      warning("unknown new coordinate system")
    }
  }

  if (has.time.comp <- length(T)!=0) {
    Ttriple <- length(T) == 3;
    if (length(T) <= 2) Tgrid <- TRUE
      else {
        dT <- diff(T)
        Tgrid <- max(abs(diff(dT))) < dT[1] * RFopt$general$gridtolerance
      }
    if (is.na(RFopt$general$Ttriple)) {
      if (Ttriple && Tgrid)
        stop("ambiguous definition of 'T'. Set RFoptions(Ttriple=TRUE) or ",
             "RFoptions(Ttriple=FALSE)")
      if (!Ttriple && !Tgrid) stop("'T' does not have a valid format")
    } else if (RFopt$general$Ttriple) {
      if (!Ttriple)
        stop("'T' is not given in triple format 'c(start, step, length)'")
      Tgrid <- FALSE
    } else {
      if (!Tgrid) stop("'T' does not define a grid")
      Ttriple <- FALSE
    }
    if (Tgrid)
      T <- as.vector(seq2grid(list(T), "T", Tgrid,
                              RFopt$internal$warn_ambiguous,
                              RFopt$general$gridtolerance))
    restotal <- restotal * T[3]
  }

   if (!missing(dim) && !is.null(dim) && spatialdim != dim) {
    stop("'dim' should be given only when 'distances' are given. Here, 'dim' contradicts the given coordinates.")
  }

  if (ex.red) {
    if (grid) {
      x[3, ] <- pmin(x[3, ], ex.red)
      if (length(y) > 0) y[3, ] <- pmin(y[3, ], ex.red)
      restotal <- as.integer(prod(x[3, ]))
    } else {
      len <- restotal <- as.integer(min(nrow(x), ex.red^spatialdim))
      x <- x[1:len, , drop=FALSE]
      if (length(y) > 0) y <- y[1:len, , drop=FALSE]
    }
    
    if (has.time.comp) {
      T[3] <- min(T[3], 3)
      restotal <- as.integer(restotal * T[3])
    }
  }
 
   
  ## keep exactly the sequence up to 'grid'
  if (length(x) > 0) {
    if (storage.mode(x) != "double") storage.mode(x) <- "double"
  } else x <- double(0)
  if (length(y) > 0) {
    if (storage.mode(y) != "double") storage.mode(y) <- "double"
  } else y <- double(0)

  L <- list(x=x, #0
            y=y, #1
            T=as.double(T), #2
            grid=as.logical(grid), #3
            spatialdim=as.integer(spatialdim), #4
            has.time.comp=has.time.comp, #5
            dist.given=FALSE, #6
            restotal=as.integer(restotal), ## 7, nr of locations
            l=as.integer(len),             ## 8, physical "length/rows" of input
            coordunits = curunits,  #9
            new_coordunits = newunits) #10
    class(L) <- "UnifyXT"
    
  return(L)  
}


trafo.to.C_UnifyXT <- function(new) {
  if (is.list(new[[1]])) {
    for(i in 1:length(new)) {
      if (length(new[[i]]$x)>0 && !new[[i]]$grid) new[[i]]$x = t(new[[i]]$x)
      if (length(new[[i]]$y)>0 && !new[[i]]$grid) new[[i]]$y = t(new[[i]]$y)
    }
  } else {
    if (length(new$x)>0 && !new$grid) new$x = t(new$x)
    if (length(new$y)>0 && !new$grid) new$y = t(new$y)
  }
  new
}


C_UnifyXT <- function(...){
  neu <- UnifyXT(...)
  return(trafo.to.C_UnifyXT(neu))
}
   
    

RFearth2cartesian <- function(coord, units=NULL, system = "cartesian",
                              grid=FALSE) {
  if (is.character(system)) system <- pmatch(system, ISO_NAMES) - 1
  stopifnot(system %in%
            c(CARTESIAN_COORD, GNOMONIC_PROJ, ORTHOGRAPHIC_PROJ))
  if (is.null(units)) {
    global.units <- RFoptions()$coords$new_coordunits[1]
    units <- if (global.units[1] == "") "km" else global.units 
  }
  if (!is.matrix(coord)) coord <- t(coord)
  res <- RFfctn(RMtrafo(new=system), coord, grid=grid,
                coords.new_coord_system = "keep",
                coords.new_coordunits=units,
                coords.coord_system="earth")
  dimnames(res) <- list(NULL, c("X", "Y", "Z", "T")[1:ncol(res)])
  return(res)
}

RFearth2dist <- function(coord, units=NULL, system="cartesian",
                         grid=FALSE, ...) {
  if (is.character(system)) system <- pmatch(system, ISO_NAMES) - 1
  stopifnot(system %in%
            c(CARTESIAN_COORD, GNOMONIC_PROJ, ORTHOGRAPHIC_PROJ))
  if (is.null(units)) {
    global.units <- RFoptions()$coords$new_coordunits[1]
    units <- if (global.units[1] == "") "km" else global.units 
  }
  if (!is.matrix(coord)) coord <- t(coord)
  z <- RFfctn(RMtrafo(new=system), coord, grid=grid,
                coords.new_coord_system = "keep",
                coords.new_coordunits=units,
                coords.coord_system="earth")
  return(dist(z, ...))
}




## used by RFratiotest, fitgauss, Crossvalidation, likelihood-ratio,  RFempir
UnifyData <- function(model, x, y=NULL, z=NULL, T=NULL,  grid, data,
                      distances=NULL, RFopt, mindist_pts=2,
                      dim=NULL, allowFirstCols=TRUE, vdim = NULL,
                      params, further.models=NULL, ...) {
  ##  if (missing(x)) Print(data, T) else Print(data, T, x)
##  if (!missing(model)) print(model)
  
  RFoptOld <- internal.rfoptions(internal.examples_reduced=FALSE)
  RFopt <- RFoptOld[[2]]
  PL <- as.integer(RFopt$basic$printlevel)
  
  if (missing(dim)) dim <- NULL
  if (missing(grid)) grid <- NULL
  
  dist.given <- !missing(distances) && length(distances)>0
##  if (dist.given) {   printf("geht nicht");  bug; }

  prepRecall <- matrix.indep.of.x.assumed <- FALSE  
  rangex <- neu <- gridlist <- RFsp.info <- mindist <- data.col <- NULL
  if (missing(data)) stop("missing data")
  missing.x <- missing(x)  
  
  if (isSpObj(data)) data <- sp2RF(data)
  if (isRFsp <- is(data, "RFsp") || (is.list(data) && is(data[[1]], "RFsp"))){
    if ( (!missing.x && length(x)!=0) || length(y)!=0   || length(z) != 0 ||
	length(T) !=  0 || dist.given || length(dim)!=0 || length(grid) != 0)
      stop("data object already contains information about the locations. So, none of 'x' 'y', 'z', 'T', 'distance', 'dim', 'grid' should be given.")
    if (!is.list(data)) data <- list(data)
    sets <- length(data)
    x <- RFsp.info <- vector("list", sets)
    
    if (!is.null(data[[1]]@.RFparams)) {
      if (length(vdim) > 0) stopifnot( vdim == data[[1]]@.RFparams$vdim)
      else vdim <- data[[1]]@.RFparams$vdim
    }
    
    ##   dimdata <- NULL
    for (i in 1:length(data)) {
      xi <- list()
      xi$grid <- isGridded(data[[i]])
      compareGridBooleans(grid, xi$grid)
      columns <- data.columns(data=data[[i]], model=model, RFopt=RFopt,
			      vdim=vdim)

##      Print(columns, vdim)
      
      sel <- columns$is.data
      if (!is.null(sel)) {	
	data[[i]] <- data[[i]][sel]
	if (!is.null(names(sel))) colnames(data[[i]]@data) <- names(sel)
##	stopifnot(!is.logical(sel))
##	data[[i]]@.RFparams$vdim <- length(sel) ### 24.12.17
      }
  
      dimensions <-
        if (xi$grid) data[[i]]@grid@cells.dim else nrow(data[[i]]@data)
      if (data[[i]]@.RFparams$vdim > 1) {
        dimensions <- c(dimensions, data[[i]]@.RFparams$vdim)      
        if (RFopt$general$vdim_close_together)
          dimensions <- dimensions[c(length(dimensions),
                                     1:(length(dimensions)-1))]
      }
      L <- nrow(data[[i]]@data) * ncol(data[[i]]@data)
#      Print(data[[i]]@data, dim(data[[i]]@data), L, dimensions, prod(dimensions))
      repet <- L / prod(dimensions)
      if (repet != as.integer(repet)) stop("number of calculated repetitions does not match the length of the data.")
      if (repet > 1) dimensions <- c(dimensions, repet)
      ##    dimdata <- rbind(dimdata, c(dimensions, data[[i]]@.RFparams$n))

      RFsp.info[[i]] <- list(data.params = data[[i]]@.RFparams,
                             dimensions = dimensions,
			     gridTopology=if (xi$grid) data[[i]]@grid else NULL,
			     coords = if (!xi$grid) data[[i]]@coords else NULL)
      
 

##      Print(data[[i]])
      
      tmp <- RFspDataFrame2conventional(data[[i]])
      xi$x <- tmp$x
      if (!is.null(tmp$T)) xi$T <- tmp$T
      data[[i]] <- as.matrix(tmp$data)   
      x[[i]] <- xi
    }   
    varnames <- names(columns$is.data)    
    coordnames <- names(columns$is.x)

  #  Print("A", coordnames)

 #   idx <- if (RFopt$general$vdim_close_together) 1 else length(dimensions)
#    if (all(dimdata[, idx] == 1))
 #     dimdata <- dimdata[, -idx, drop=FALSE]
 #   if (all(dimdata[, ncol(dimdata)] == 1)) # repet
 #     dimdata <- dimdata[, -ncol(dimdata), drop=FALSE]
    
  } else { # !isRFsp
    ## dimdata wird spaeter bestimmt

    if (is.data.frame(data) || !is.list(data)) data <- list(data)
    sets <- length(data)
    for (i in 1:sets) {
      if (is.data.frame(data[[i]]) || is.vector(data[[i]]))
	data[[i]] <- as.matrix(data[[i]])
    }
    
#Print(data[[1]], model, xdim=dim,
#			    force=allowFirstCols, halt=!allowFirstCols)

    columns <- data.columns(data[[1]], model=model, RFopt = RFopt, xdim=dim,
			    force=allowFirstCols, halt=!allowFirstCols,
			    vdim=vdim)
    varnames <- names(columns$is.data)
    coordnames <- names(columns$is.x)
#    Print("B", coordnames)
  
    if (dist.given) {
      stopifnot(missing(x) || length(x)==0, length(y)==0, length(z)==0)
      if (!is.list(distances)) distances <- list(distances)
      if (length(distances) != length(data))
          stop("number of sets of distances does not match number of sets of data")
      for (i in 1:length(distances)) {
        if (any(is.na(data)))
          stop("missing data are not allowed if distances are used.")
      }

      stopifnot(missing(T) || length(T)==0)
      if (is.matrix(distances[[1]])) {

        dimensions <- sapply(distances, nrow)
        spatialdim <- tsdim <- xdimOZ <- dimensions[1]
        if (length(dim) > 0 && dim != spatialdim)
          stop("unclear specification of the distances: either the distances is given as a vector or distance vectors should given, where the number of rows matches the spatial dimension")
        lcc <- sapply(distances, function(x) 0.5 * (1 + sqrt(1 + 8 * ncol(x))) )
        if (!all(diff(dimensions) == 0))
          stop("sets of distances show different dimensions")
        range_distSq <- function(M) range(apply(M, 2, function(z) sum(z^2)))
        rangex <- sqrt(range(sapply(distances, range_distSq)))
      } else {        
        xdimOZ <- 1L
        spatialdim <- tsdim <- as.integer(dim)      
        lcc <- sapply(distances, function(x) if (is.matrix(x)) -1
                                        else 0.5 * (1 + sqrt(1 + 8* length(x))))
        rangex <- range(sapply(distances, range))
      }
#      Print(mindist, rangex, RFopt$nugget$tol)      
      mindist <- min(rangex)
      if (is.na(mindist)) mindist <- 1 ## nur 1 pkt gegeben, arbitraerer Wert
      if (mindist <= RFopt$nugget$tol) {
        if (!RFopt$general$allowdistanceZero)
          stop("distance with value 0 identified -- use allowdistanceZero=T?")
        mindist <- 1e-15 * (RFopt$nugget$tol == 0) + 2 * RFopt$nugget$tol

        for (i in 1:length(distances))
          if (is.vector(distances[[i]]))
            distances[[i]][distances[[i]] == 0] <- mindist
          else distances[[i]][1, apply(distances[[i]], 2,
                                       function(z) sum(z^2))] <- mindist
      }

      if (any(as.integer(lcc) != lcc))
	stop("number of distances not of form k(k-1)/2")
      neu <- UnifyXT(distances=distances, dim = spatialdim) 
      coordunits <- RFopt$coords$coordunits
      has.time.comp <- FALSE      
    } else { ## distances not given
      
       
      if (missing(x)) { ## dec 2012: matrix.indep.of.x.assumed        
	x <- list()
	if (length(columns$is.x)==0) {
	  matrix.indep.of.x.assumed <- TRUE
	  if (PL > 0)
	    message("Columns representing coordinates not found. So no spatial context is assumed.")
	  for (i in 1:sets) {
	    x[[i]] <- 1:nrow(as.matrix(data[[i]]))
	    storage.mode(x[[i]]) <- "numeric"
	  }
	} else {
	  for (i in 1:sets) {
	    xx <- data[[i]][, columns$is.x, drop=FALSE]
	    storage.mode(xx) <- "numeric"
	    colnames(xx) <- coordnames
	    x[[i]] <- list(x=xx, grid=FALSE)
	  }
        } ## end mising
        
        for (i in 1:sets) {
          data[[i]] <- data[[i]][ , columns$is.data, drop=FALSE]
	  colnames(data[[i]]) <- varnames
         storage.mode(data[[i]]) <- "numeric"
        }
        
      } ## missing xgiven; KEIN ELSE, auch wenn nachfolgend z.T. gedoppelt wird

    
      if (is.data.frame(x)) x <- as.matrix(x)   
      if (is.list(x)) {
        if (length(y)!=0 || length(z)!=0 || length(T)!=0)
          stop("if x is a list, then 'y', 'z', 'T' may not be given")       

       if (!is.list(x[[1]])) {
         if (length(data) == 1) x <- list(x)
         else stop("number of sets of 'x' and 'data' differ")
       }
       
        
       } else {
        x <- list(x=x)
        if (length(y)!=0) {
          stopifnot(!is.list(y))
          x$y <- y
        }
        if (length(z)!=0) {
          stopifnot(!is.list(z))
          x$z <- z
        }
        if (length(T)!=0) {
          stopifnot(!is.list(T))
          x$T <- T
        }
        if (!is.null(grid))
          x$grid <- grid
        if (!is.list(data)) data <- list(as.matrix(data))
        x <- list(x)          
      }
      ## {}
    } # ! distance
    sets <- length(data)

 ##   next * two * linse * also * previous * vesrion
  ##  dimdata <- matrix(nrow=sets, ncol=length(base::dim(data[[1]])))
  ##  for (i in 1:sets) dimdata[i, ] <- base::dim(data[[i]])
    
  } # !isRFsp

  attr(data, "Unified") <- TRUE
 
  if (!dist.given) { ##   x coordinates, not distances

    neu <- UnifyXT(x=x, printlevel=min(PL, PL_IMPORTANT)) #, y=y, z=z, T=T, grid=grid, distances=distances,
    ## dim=dim, length # , length.data=length(data[[i]]), printlevel = 0

    if (!is.list(neu[[1]])) neu <- list(neu)

    coordunits<- neu[[1]]$coordunits
    spatialdim <- as.integer(neu[[1]]$spatialdim)
    has.time.comp <- neu[[1]]$has.time.comp
    tsdim <- as.integer(spatialdim + has.time.comp)

    getrange <- function(x)
      if (x$grid) rbind(x$x[1, ], x$x[1, ] + x$x[2, ] * (x$x[3, ] - 1))
      else apply(x$x, 2, range)
    rangex <- sapply(neu, getrange)

    ## falls mehrere datasets:
    if (ncol(x[[1]]$x) > 1 || is.null(x[[1]]$dist.given) || !x[[1]]$dist.given){
      rangex <- t(rangex) 
      base::dim(rangex) <- c(length(rangex) / spatialdim, spatialdim)
    }
    rangex <- apply(rangex, 2, range)

    getmindistSq <- function(x) {
      if (x$grid) sum(x$x[2,]^2)
      else if (nrow(x$x) < 2) NA
      else if (nrow(x$x) <= mindist_pts) min(dist(x$x))
      else min(dist(x$x[sample(nrow(x$x), mindist_pts), ]))
    }

    if (has.time.comp && any(sapply(neu, function(x) x$T[2]) <= RFopt$nugget$tol))
      stop("step of time component smaller than nugget tolerance 'tol'")
    
    if (any(sapply(neu, function(x) x$grid && any(x$x[2, ]<=RFopt$nugget$tol))))
      stop("step of some spatial component smaller than nugget tolerance 'tol'")

    zaehler <- 0

    repeat {
      mindist <- sqrt(min(sapply(neu, getmindistSq)))      
      if (is.na(mindist)) mindist <- 1 ## nur 1 pkt gegeben, arbitraerer Wert
      if (mindist <= RFopt$nugget$tol) {
        if (!RFopt$general$allowdistanceZero)
          stop("Distance with value 0 identified -- use allowdistanceZero=T?")
        if ((zaehler <- zaehler + 1) > 10)
          stop("unable to scatter point pattern")
        for (i in 1:length(neu)) if (!neu[[i]]$grid)
          neu[[i]]$x <- neu[[i]]$x + rnorm(length(neu[[i]]$x), 0,
                                           10 * RFopt$nugget$tol)
      } else break;
    }
  
    xdimOZ <- ncol(neu[[1]]$x)
  }
 
 
  ## geht x[[1]]$x immer gut ??
  ##  Print(missing(x), neu)

#  Print("CD", coordnames, neu[[1]])
  if (is.null(coordnames))
    coordnames <- SystemCoordNames(locinfo=neu[[1]], RFopt=RFopt) 

#     Print("C", coordnames)
 
  restotal <- sapply(neu, function(x) x$restotal)
  ldata <- sapply(data, length)

  ## unklar wo dies gebraucht wird. Macht schwierigkeiten 27.2.19. Warum
  ## wird onedim.x gebraucht? mit false mal ausgeschaltet
  onedim.x <- FALSE && all(sapply(data, function(x) is.vector(x) || ncol(x) == 1))

  
  if (missing(model)) {
    if (length(further.models) > 0)
      stop("'model' is not given, but 'further.models'")
    if (length(vdim) == 0) vdim <- if (onedim.x) 1 else 1# warum war da else NA?
  } else {
    model.vdim <- integer(1)
    neu_C <- trafo.to.C_UnifyXT(neu)

    # internal <- RFoptions(GETOPTIONS="internal", declare_PL = 0)
    newmodel <- PrepareModel2(model=model, ..., params=params,
                              allow.formulae = TRUE, x=neu_C)

    if (is(newmodel, "RM_model")) {
      values <- double(0)
      pos <- integer(0)
    } else {
      values <- newmodel$values
      newmodel <- newmodel$model
    }

    pos <- .Call(C_GetNAPositions, MODEL_AUX, list("Cov", newmodel), neu_C,
                 as.double(values),
                 ## spatialdim, has.time.comp,  xdimOZ = xdimOZ,
                 FALSE, model.vdim,  ##           set by side effect !
                 PL - 5L
                 )    

    ##    Print(pos, newmodel, values, pos,further.models, model.vdim);# kkk
   
    if (prepRecall <- length(values) > 0 || length(further.models) > 0) {
      ans <- PrepareModel2(model=model, ..., params=params,
                           allow.formulae = TRUE, indices = pos, x=neu_C)
      newmodel <- ans$model
    }

    if (!missing(params) && length(further.models)>0) {
      #RFoptions(GETOPTIONS="internal", declare_PL = PL)
      for (i in 1:length(further.models)) {
        further.models[[i]] <-
          PrepareModel2(model=model, ..., params=params, allow.formulae = FALSE,
                        further.models = further.models[[i]],
                        arg = names(further.models)[[i]],
                        x=neu_C)
      }
    }

##    Print(vdim, onedim.x, model.vdim)
    
    if (length(vdim) == 0) vdim <- if (onedim.x) 1 else model.vdim
    else if (vdim != model.vdim)
      stop("given multivariate dimension differs from the dimensions expected by the model")  
  }
  
    
  repetitions <- as.integer(ldata / (restotal * vdim))
                                        #
##  Print(data, ldata, repetitions, restotal, vdim, neu, dist.given)
  if (any(repetitions)==0) stop("no or not sufficiently many data are given")
  if (!is.na(vdim) && any(ldata != repetitions * restotal * vdim))
    stop("mismatch of data dimensions")

  vrep <- repetitions * vdim

  ##Print(ldata, repetitions, restotal, vdim, vrep, data)

  for (i in 1:length(data)) {
    base::dim(data[[i]]) <- c(ldata[i] / vrep[i], vrep[i])
  }
 
  ## Print(data, length(data), repetitions, vdim)
  ## data <- lapply(data, dim(data) <- c(length(data) / (repetitions * vdim) , repetitions * vdim))
  
  
  return(list(
    ## coord = expandiertes neu # #
    model = if (missing(model)) NULL else newmodel,
    orig.model = if (missing(model)) NULL else model,
    further.models = further.models,
    transform = if (prepRecall) ans$transform,
 ##   dimdata=dimdata, isRFsp = isRFsp,  # 3.2.0 : not given anymore
## oben auch noch loeschen!!
 ##  len = len,    # 3.2.0 : not given anymore
    data=data,
    RFsp.info = RFsp.info,
    coord = neu,
    dist.given=dist.given,
    spatialdim=spatialdim,
    tsdim=tsdim,
    rangex = as.matrix(rangex),
    coordunits=coordunits,
    has.time.comp = has.time.comp,
    matrix.indep.of.x.assumed = matrix.indep.of.x.assumed,
    mindist = mindist,
    xdimOZ = xdimOZ,
    vdim = vdim,
    coordnames=coordnames,
    varnames=varnames,
    data.col = columns,
    repetitions = repetitions
  ))
}