From cd12b795883d3e1095421626d4b97878ba5669ca Mon Sep 17 00:00:00 2001
From: Raimon Tolosana-Delgado <tolosa53@fwg206.ad.fz-rossendorf.de>
Date: Thu, 15 Jul 2021 17:52:43 +0200
Subject: [PATCH] gmGeostats v0.10-8 submission

---
 DESCRIPTION                   | 15 +++++---
 NAMESPACE                     | 11 +++---
 R/Anamorphosis.R              |  6 +--
 R/abstractClasses.R           | 14 +++----
 R/accuracy.R                  |  7 ++--
 R/closeup.R                   |  0
 R/compositionsCompatibility.R | 60 ++++++++++++++++--------------
 R/genDiag.R                   |  3 +-
 R/geostats.R                  | 12 ++----
 R/gmSimulation.R              |  4 +-
 R/gmSpatialMethodParameters.R | 10 ++---
 R/gmSpatialModel.R            | 39 +++++++++++--------
 R/gmValidationStrategy.R      | 70 +++++++++++++----------------------
 R/grids.R                     |  2 +-
 R/gstatCompatibility.R        |  2 +-
 R/preparations.R              | 47 +++++++++++++++++++++++
 R/variograms.R                | 21 ++++-------
 src/gmGeostats.c              |  6 +--
 18 files changed, 180 insertions(+), 149 deletions(-)
 create mode 100644 R/closeup.R
 create mode 100644 R/preparations.R

diff --git a/DESCRIPTION b/DESCRIPTION
index ff80217..3bf20c5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: gmGeostats
 Version: 0.10-8
-Date: 2020-10-05
+Date: 2021-07-15
 Title: Geostatistics for Compositional Analysis
 Authors@R: c(person(given = "Raimon", 
         family = "Tolosana-Delgado", 
@@ -32,7 +32,8 @@ Suggests:
     knitr,
     rmarkdown (>= 2.3),
     magrittr,
-    readxl
+    readxl,
+    RandomFields
 Imports: 
     methods, 
     gstat, 
@@ -60,21 +61,23 @@ VignetteBuilder: knitr
 LazyData: true
 Collate: 
     'Anamorphosis.R'
-    'compositionsCompatibility.R'
+    'preparations.R'
     'gstatCompatibility.R'
-    'variograms.R'
     'gmAnisotropy.R'
-    'gmValidationStrategy.R'
+    'compositionsCompatibility.R'
+    'variograms.R'
+    'gmSpatialMethodParameters.R'
     'abstractClasses.R'
     'accuracy.R'
+    'closeup.R'
     'data.R'
     'exploratools.R'
     'genDiag.R'
     'geostats.R'
     'gmDataFrameStack.R'
     'gmSimulation.R'
-    'gmSpatialMethodParameters.R'
     'gmSpatialModel.R'
+    'gmValidationStrategy.R'
     'grids.R'
     'mask.R'
     'spSupport.R'
diff --git a/NAMESPACE b/NAMESPACE
index 9874ed5..feecde3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -99,6 +99,8 @@ S3method(plot,gmEVario)
 S3method(plot,logratioVariogramAnisotropy)
 S3method(plot,swarmPlot)
 S3method(precision,accuracy)
+S3method(predict,LMCAnisCompo)
+S3method(predict,genDiag)
 S3method(predict,gmCgram)
 S3method(print,mask)
 S3method(setMask,DataFrameStack)
@@ -124,6 +126,8 @@ S3method(swath,rcomp)
 S3method(unmask,DataFrameStack)
 S3method(unmask,SpatialPixels)
 S3method(unmask,data.frame)
+S3method(validate,LeaveOneOut)
+S3method(validate,NfoldCrossValidation)
 S3method(variogramModelPlot,gmEVario)
 S3method(variogramModelPlot,gstatVariogram)
 S3method(variogramModelPlot,logratioVariogram)
@@ -139,6 +143,7 @@ export(LMCAnisCompo)
 export(LeaveOneOut)
 export(Maf)
 export(NfoldCrossValidation)
+export(Predict)
 export(RJD)
 export(SequentialSimulation)
 export(TurningBands)
@@ -191,8 +196,6 @@ export(noStackDim)
 export(pairsmap)
 export(precision)
 export(predict)
-export(predict.LMCAnisCompo)
-export(predict.genDiag)
 export(pwlrmap)
 export(setCgram)
 export(setGridOrder)
@@ -211,7 +214,6 @@ export(swarmPlot)
 export(swath)
 export(unmask)
 export(validate)
-export(validate_gmSpatialModel)
 export(variogram)
 export(variogramModelPlot)
 export(variogram_gmSpatialModel)
@@ -231,11 +233,10 @@ exportClasses(gmSpatialModel)
 exportClasses(gmTrainingImage)
 exportClasses(gmUnconditionalSpatialModel)
 exportClasses(gmValidationStrategy)
+exportMethods(Predict)
 exportMethods(as.gstat)
 exportMethods(logratioVariogram)
-exportMethods(predict)
 exportMethods(stackDim)
-exportMethods(validate)
 exportMethods(variogram)
 import(RColorBrewer)
 import(compositions)
diff --git a/R/Anamorphosis.R b/R/Anamorphosis.R
index 488cd44..61ea31a 100644
--- a/R/Anamorphosis.R
+++ b/R/Anamorphosis.R
@@ -153,8 +153,7 @@ anaForward <- function(x,Y,sigma0,sigma1=1+sigma0,steps=30,plt=FALSE,sphere=TRUE
      wY=checkDouble(weights,ncol(Y)),
      steps=checkInt(steps,1),
      sigma0=checkDouble(sigma0,1),
-     sigma1=checkDouble(sigma1,1),
-     PACKAGE = "gmGeostats"
+     sigma1=checkDouble(sigma1,1)#,     PACKAGE = "gmGeostats"
      )
   x<- t(structure(erg$x, dim=dim(x)))
   colnames(x) = paste("flow", 1:ncol(x), sep="")
@@ -237,8 +236,7 @@ anaBackward <- function(x,Y,sigma0,sigma1=1+sigma0,steps=30,plt=FALSE,sphere=TRU
      wY=checkDouble(weights,ncol(Y)),
      steps=checkInt(steps,1),
      sigma0=checkDouble(sigma0,1),
-     sigma1=checkDouble(sigma1,1),
-     PACKAGE = "gmGeostats"
+     sigma1=checkDouble(sigma1,1) #,     PACKAGE = "gmGeostats"
      )
   x<- t(structure(erg$x, dim=dim(x)))
   st(x,inv=TRUE)
diff --git a/R/abstractClasses.R b/R/abstractClasses.R
index 51c9155..cbec8af 100644
--- a/R/abstractClasses.R
+++ b/R/abstractClasses.R
@@ -4,14 +4,6 @@
 # S3 -> S4 classes
 # cat("creating spatial method parameter classes\n")
 
-#' @include gmAnisotropy.R
-#' @include gmValidationStrategy.R 
-#' @include variograms.R
-setOldClass("gmKrigingNeighbourhood")
-setOldClass("gmDirectSamplingParameters")
-setOldClass("gmTurningBands")
-setOldClass("gmSequentialSimulation")
-setOldClass("gmCholeskyDecomposition")
 
 
 # abstract classes
@@ -20,16 +12,19 @@ setOldClass("gmCholeskyDecomposition")
 #' @title Neighbourhood description
 #' @description abstract class, containing any specification of a spatial neighbourhood
 #' @export
+#' @include gmSpatialMethodParameters.R
+#' @include preparations.R
 setClassUnion(name="gmNeighbourhoodSpecification", members=c("gmKrigingNeighbourhood","NULL"))
 
 #' @title Validation strategy description
 #' @description abstract class, containing any specification of a validation strategy for spatial models
 #' @export
 setClassUnion(name="gmValidationStrategy", 
-              members=c("NULL",
+              members=c("NULL", 
                         "LeaveOneOut", 
                         "NfoldCrossValidation"))
 
+
 #' @title parameters for Multiple-Point Statistics methods
 #' @description abstract class, containing any parameter specification of a spatial multipoint algorithm 
 #' @export
@@ -97,6 +92,7 @@ setClassUnion(name="gmTrainingImage",
 #' @description abstract class, containing any specification of an unconditional
 #' spatial model 
 #' @export
+#' @include variograms.R
 setClassUnion(name="gmUnconditionalSpatialModel", 
               members=c("NULL",
                         "gmGaussianModel",
diff --git a/R/accuracy.R b/R/accuracy.R
index 742265a..4a2e6cf 100644
--- a/R/accuracy.R
+++ b/R/accuracy.R
@@ -142,7 +142,7 @@ accuracy.data.frame <- function(x, observed=x$observed,
 accuracy.DataFrameStack <- function(x, observed, 
                                     ivars = intersect(colnames(observed), dimnames(x)[[noStackDim(x)]]),
                                     prob = seq(from=0, to=1, by=0.05),
-                                    method = ifelse(length(vars)==1, "simulation", "Mahalanobis"),
+                                    method = ifelse(length(ivars)==1, "simulation", "Mahalanobis"),
                                     outMahalanobis=FALSE, ...){
   methods = c("simulation", "Mahalanobis","mahalanobis", "flow")
   mm = methods[pmatch(method, methods)]
@@ -329,13 +329,13 @@ xvErrorMeasures <- function(x,...) UseMethod("xvErrorMeasures", x)
 #' @param krigVar a vector containing the kriging variances
 #' @param observed a vector containing the true values
 #' @param output  which output do you want? a vector of one or several of  c("ME","MSE","MSDR","Mahalanobis")
+#' @param ... extra arguments for generic functionality
 #'
-#' @return
 #' 
 #' @export
 #' @method xvErrorMeasures default
 #' @family accuracy functions
-xvErrorMeasures.default <- function(x, krigVar, observed, output="MSDR1"){
+xvErrorMeasures.default <- function(x, krigVar, observed, output="MSDR1", ...){
   if(length(output)>1)
     return(sapply(output, function(o) xvErrorMeasures(x, krigVar, observed, output=o)))
   
@@ -359,6 +359,7 @@ xvErrorMeasures.default <- function(x, krigVar, observed, output="MSDR1"){
 #' @param observed a vector (if univariate) or a matrix/dataset of true values
 #' @param output which output do you want? a vector of one or several of  c("ME","MSE","MSDR","MSDR1","MSDR2","Mahalanobis")
 #' @param univariate logical control, typically you should not touch it
+#' @param ... extra arguments for generic functionality
 #'
 #' @return If just some of c("ME","MSE","MSDR","MSDR1","MSDR2") are requested, the output is a named
 #' vector with the desired quantities. If only "Mahalanobis" is requested, the output is a vector 
diff --git a/R/closeup.R b/R/closeup.R
new file mode 100644
index 0000000..e69de29
diff --git a/R/compositionsCompatibility.R b/R/compositionsCompatibility.R
index 10344e1..86e9b81 100644
--- a/R/compositionsCompatibility.R
+++ b/R/compositionsCompatibility.R
@@ -33,7 +33,7 @@ logratioVariogram <- function(data, ...) UseMethod("logratioVariogram", data)
 #' @method logratioVariogram default
 #' @export
 #' @param loc if `data` is a composition (or if `comp` is provided), spatial coordinates of the sampling locations 
-#' @param comp an alias for `data`, provided for back-compatibility with [compositions::logratioVariogram()] 
+#' @param comp an alias for `data`, provided for back-compatibility with compositions::logratioVariogram 
 logratioVariogram.default <- function(data, loc, ..., comp=data){
   res = try(compositions::logratioVariogram(comp=acomp(comp), loc=loc, ...), silent=TRUE)
   if(class(res)!="try-error") return(res)
@@ -66,7 +66,7 @@ logratioVariogram_gmSpatialModel <- function(data, ..., azimuth=0, azimuth.tol=1
   logratioVariogram.default(comp, loc=coords, ...)    
 }
 
-
+#' @describeIn logratioVariogram S4 generic for gmSpatialModel objects 
 setGeneric("logratioVariogram", logratioVariogram)
 
 
@@ -706,13 +706,14 @@ LMCAnisCompo = function(Z, models=c("nugget", "sph", "sph"),
 #'
 #' @return an array of dimension (nr of lags, D, D) with D the number of variables in the model.
 #' @export
-#'
+#' @include gmAnisotropy.R
+#' @method predict LMCAnisCompo
 #' @examples
 #' data("jura", package="gstat")
 #' Zc = compositions::acomp(jura.pred[,7:9])
 #' lrmd = LMCAnisCompo(Zc, models=c("nugget", "sph"), azimuths=c(0,45))
 #' predict(lrmd, outer(0:5, c(0,1)))
-predict.LMCAnisCompo <- function(object, newdata, ...){
+predict.LMCAnisCompo = function(object, newdata, ...){
   K = ncol(object)
   if(is.data.frame(newdata)) newdata = as.matrix(newdata)
   if(is.null(dim(newdata))) newdata = cbind(newdata, 0)
@@ -735,6 +736,8 @@ predict.LMCAnisCompo <- function(object, newdata, ...){
 
 
 
+
+
 #' Recast a model to the variogram model of package "compositions"
 #' 
 #' Recast a variogram model specified in any of the models of "gstat" or "gmGeostats" in
@@ -1298,30 +1301,31 @@ as.gmCgram.LMCAnisCompo = function(m, V=attr(m,"contrasts"),
   return(output)
 }
 
-
-# @describeIn as.variogramModel
-as.variogramModel.CompLinModCoReg <- function(m, V=NULL, ...){
-  stop("not yet available")
-  # extract the substructures from m-variogram
-  rs = gsi.extractCompLMCstructures(m) # elements: "models", "ranges", "sills"
-  # construct the vgm template
-    # 1.- set the nugget (if needed) 
-  if(any(rs$models=="nugget")){
-    vg0 = vgm(psill=1, model="Nug")
-  }else{
-    vg0 = NULL
-  }
-   # 2.- add each structure
-  if(sum(rs$models!="nugget")>0){
-    for(i in which(rs$models!="nugget")){
-      vg0 = vgm(add.to = vg0, model=rs$models[i], range = rs$ranges[i], psill=1)
-    }
-  }
-  # cast the sills to the requested logratio coordinates
-  
-  # expand the vgm template to the new coordinates
-  
-} 
+# 
+# # @describeIn as.variogramModel
+#  EXISTS in gstatCompatibility.R
+# as.variogramModel.CompLinModCoReg <- function(m, V=NULL, ...){
+#   stop("not yet available")
+#   # extract the substructures from m-variogram
+#   rs = gsi.extractCompLMCstructures(m) # elements: "models", "ranges", "sills"
+#   # construct the vgm template
+#     # 1.- set the nugget (if needed) 
+#   if(any(rs$models=="nugget")){
+#     vg0 = vgm(psill=1, model="Nug")
+#   }else{
+#     vg0 = NULL
+#   }
+#    # 2.- add each structure
+#   if(sum(rs$models!="nugget")>0){
+#     for(i in which(rs$models!="nugget")){
+#       vg0 = vgm(add.to = vg0, model=rs$models[i], range = rs$ranges[i], psill=1)
+#     }
+#   }
+#   # cast the sills to the requested logratio coordinates
+#   
+#   # expand the vgm template to the new coordinates
+#   
+# } 
 
 
 #' extract information about the original data, if available
diff --git a/R/genDiag.R b/R/genDiag.R
index dd65d38..dd2b516 100644
--- a/R/genDiag.R
+++ b/R/genDiag.R
@@ -430,7 +430,9 @@ RJD.rcomp <- RJD.acomp
 #'
 #' @return A data set or compositional object of the nature of the original data
 #' used for creating the genDiag object.
+#' @include gmAnisotropy.R
 #' @export
+#' @method predict genDiag
 #' @family generalised Diagonalisations
 #' @examples
 #' data("jura", package="gstat")
@@ -459,7 +461,6 @@ predict.genDiag = function (object, newdata=NULL, ...) {
 }
 
 
-
 #' Colored biplot for gemeralised diagonalisations
 #' Colored biplot method for objects of class genDiag
 #' @param x a generalized diagonalisation object, as obtained from a call to
diff --git a/R/geostats.R b/R/geostats.R
index 32f495f..46bed0e 100644
--- a/R/geostats.R
+++ b/R/geostats.R
@@ -65,8 +65,7 @@ gsi.calcCgram <- function(X,Y,vgram,ijEqual=FALSE) {
            A = checkDouble(t(apply(vgram$M,1,gsi.powM,alpha=-1/2)),k*m*m),
            Sill=checkDouble(vgram$sill,k*d*d),
            moreCgramData=checkDouble(vgram$data,k),
-           ijEqual=checkLogical(ijEqual,1),
-           PACKAGE = "gmGeostats"
+           ijEqual=checkLogical(ijEqual,1)#,     PACKAGE = "gmGeostats"
   )
   structure(erg$C,dim=c(d*nX,d*nY))
 }
@@ -127,8 +126,7 @@ gsi.TurningBands <- function(X,vgram,nBands,nsim=NULL) {
             typeCgram = checkInt(vgram$type,k),
             A = checkDouble(t(apply(vgram$M,1,gsi.powM,alpha=-1/2)),k*m*m),
             sqrtSill=checkDouble(t(apply(vgram$sill,1,gsi.powM,alpha=1/2)),k*d*d),
-            moreCgramData=checkDouble(vgram$data,k),
-            PACKAGE = "gmGeostats"
+            moreCgramData=checkDouble(vgram$data,k)#,     PACKAGE = "gmGeostats"
     )
       erg = cbind(X[, 1:m0],Z=t(structure(erg$Z,dim=c(d,nrow(X)))))
       return(erg) 
@@ -146,8 +144,7 @@ gsi.TurningBands <- function(X,vgram,nBands,nsim=NULL) {
             typeCgram = checkInt(vgram$type,k),
             A = checkDouble(t(apply(vgram$M,1,gsi.powM,alpha=-1/2)),k*m*m),
             sqrtSill=checkDouble(t(apply(vgram$sill,1,gsi.powM,alpha=1/2)),k*d*d),
-            moreCgramData=checkDouble(vgram$data,k),
-            PACKAGE = "gmGeostats"
+            moreCgramData=checkDouble(vgram$data,k)#,     PACKAGE = "gmGeostats"
     )
     erg = aperm(structure(erg$Z,dim=c(d,nrow(X),nsim)),c(2,1,3))
     return(erg) 
@@ -226,8 +223,7 @@ gsi.CondTurningBands <- function(Xout, Xin, Zin, vgram,
             Sill=checkDouble(vgram$sill,k*d*d),
             moreCgramData=checkDouble(vgram$data,k),
             cbuf=numeric(d*d*nin),
-            dbuf=numeric(d*nin*nsim),
-            PACKAGE = "gmGeostats"
+            dbuf=numeric(d*nin*nsim)#,     PACKAGE = "gmGeostats"
   )
   if( normal ) {
     cbind(Xout,t(structure(erg$Z,dim=dimZ))[-(1:nrow(Xin)),])
diff --git a/R/gmSimulation.R b/R/gmSimulation.R
index 88db6fd..cdbb135 100644
--- a/R/gmSimulation.R
+++ b/R/gmSimulation.R
@@ -20,7 +20,7 @@ gsi.DS4CoDa <- function(n, f, t, n_realiz, nx_TI, ny_TI, nx_SimGrid, ny_SimGrid,
                         SimGrid_mask = ncol(SimGrid_input), invertMask = TRUE
                         ){
   .Deprecated(new = "gsi.DS", 
-              msg="gsi.DS4CoDa is deprecated; use gsi.DS or go via make.gm*-functions followed by DSpars and predict.gmSpatialModel")
+              msg="gsi.DS4CoDa is deprecated; use gsi.DS or go via make.gm*-functions followed by DSpars and predict for gmSpatialModel")
   ### extract elements
   # mask
   if(length(SimGrid_mask)==1){
@@ -226,7 +226,7 @@ gsi.DS4CoDa <- function(n, f, t, n_realiz, nx_TI, ny_TI, nx_SimGrid, ny_SimGrid,
 #' grid is simulated. The '@data' slot of these objects contains a [DataFrameStack()] with the stacking dimension
 #' running through the realisations. It is safer to use this functionality through the interface
 #' [make.gmCompositionalMPSSpatialModel()], then request a direct simulation with [DSpars()] and
-#' finally run it with [predict.gmSpatialModel()].
+#' finally run it with [predict_gmSpatialModel].
 #' @export
 #' @importFrom utils txtProgressBar setTxtProgressBar
 #' @importFrom stats complete.cases lm 
diff --git a/R/gmSpatialMethodParameters.R b/R/gmSpatialMethodParameters.R
index dec722f..5934247 100644
--- a/R/gmSpatialMethodParameters.R
+++ b/R/gmSpatialMethodParameters.R
@@ -16,7 +16,7 @@
 #' @param ... further arguments, currently ignored
 #'
 #' @return an S3-list of class "gmKrigingNeighbourhood" containing the six elements given as arguments 
-#' to the function. This is just a compact way to provide further functions such as [predict.gmSpatialModel()]
+#' to the function. This is just a compact way to provide further functions such as [predict_gmSpatialModel]
 #' with appropriate triggers for choosing a prediction method or another, in this case for triggering
 #' cokriging (if alone) or eventually sequential simulation (see [SequentialSimulation()]).
 #' @export
@@ -64,7 +64,7 @@ KrigingNeighbourhood <- function(nmax=Inf, nmin=0, omax=0, maxdist=Inf, force=FA
 #' @param ... further parameters, not used
 #'
 #' @return an S3-list of class "gmDirectSamplingParameters" containing the six elements given as arguments 
-#' to the function. This is just a compact way to provide further functions such as [predict.gmSpatialModel()]
+#' to the function. This is just a compact way to provide further functions such as [predict_gmSpatialModel]
 #' with appropriate triggers for choosing a prediction method or another, in this case for triggering 
 #' direct sampling.
 #' @export
@@ -96,7 +96,7 @@ DSpars <- DirectSamplingParameters <- function(nsim=1, scanFraction=0.25, patter
 #' @param ... further parameters, currently ignored
 #'
 #' @return an S3-list of class "gmSequentialSimulation" containing the four elements given as arguments 
-#' to the function. This is just a compact way to provide further functions such as [predict.gmSpatialModel()]
+#' to the function. This is just a compact way to provide further functions such as [predict_gmSpatialModel]
 #' with appropriate triggers for choosing a prediction method or another, in this case for triggering 
 #' sequential Gaussian simulation.
 #' @export
@@ -128,7 +128,7 @@ SequentialSimulation = function(nsim=1, ng=NULL, rank=Inf, debug.level=1, ...){
 #' @param ... further parameters, currently ignored
 #'
 #' @return an S3-list of class "gmTurningBands" containing the few elements given as arguments 
-#' to the function. This is just a compact way to provide further functions such as [predict.gmSpatialModel()]
+#' to the function. This is just a compact way to provide further functions such as [predict_gmSpatialModel]
 #' with appropriate triggers for choosing a prediction method or another, in this case for triggering 
 #' turning bands simulation.
 #' @export
@@ -152,7 +152,7 @@ TurningBands = function(nsim=1, nBands=1000,  ...){
 #' @param ... further parameters, currently ignored
 #'
 #' @return an S3-list of class "gmCholeskyDecomposition" containing the few elements given as arguments 
-#' to the function. This is just a compact way to provide further functions such as [predict.gmSpatialModel()]
+#' to the function. This is just a compact way to provide further functions such as [predict_gmSpatialModel]
 #' with appropriate triggers for choosing a prediction method or another, in this case for triggering 
 #' LU or Cholesky decomposition simulation.
 #' @export
diff --git a/R/gmSpatialModel.R b/R/gmSpatialModel.R
index d4cb2d2..89bd9e1 100644
--- a/R/gmSpatialModel.R
+++ b/R/gmSpatialModel.R
@@ -64,7 +64,7 @@ setClass("gmSpatialModel", contains="SpatialPointsDataFrame",
 #' @export
 #' @family gmSpatialModel 
 #' @seealso [SequentialSimulation()], [TurningBands()] or [CholeskyDecomposition()] for specifying the exact 
-#' simulation method and its parameters, [predict.gmSpatialModel()] for running predictions or simulations
+#' simulation method and its parameters, [predict_gmSpatialModel] for running predictions or simulations
 #'
 #' @examples
 #' data("jura", package="gstat")
@@ -120,7 +120,7 @@ make.gmMultivariateGaussianSpatialModel <- function(
 #' @export
 #' @family gmSpatialModel
 #' @seealso [SequentialSimulation()], [TurningBands()] or [CholeskyDecomposition()] for specifying the exact 
-#' simulation method and its parameters, [predict.gmSpatialModel()] for running predictions or simulations
+#' simulation method and its parameters, [predict_gmSpatialModel] for running predictions or simulations
 #'
 #' @examples
 #' data("jura", package="gstat")
@@ -181,7 +181,7 @@ make.gmCompositionalGaussianSpatialModel <- function(
 #' @export
 #' @family gmSpatialModel 
 #' @seealso [DirectSamplingParameters()] for specifying a direct simulation method parameters,
-#' [predict.gmSpatialModel()] for running the simulation
+#' [predict_gmSpatialModel] for running the simulation
 make.gmCompositionalMPSSpatialModel = function(
   data, coords = attr(data, "coords"), V="ilr", prefix=NULL,  ## data trench
   model=NULL)   ## model trench
@@ -401,15 +401,18 @@ as.gmSpatialModel.gstat = function(object, V=NULL, ...){
 NULL
 
 #' @rdname predict_gmSpatialModel
-#' @export
-setGeneric("predict", function(object, newdata, pars, ...){
-  standardGeneric("predict")
-})
+setMethod("predict", signature(object="gmSpatialModel"),
+          function(object, newdata, pars, ...){
+            Predict(object, newdata, pars, ...)
+          })
+          
 
 
 #' @rdname predict_gmSpatialModel
-#' @export
-setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="missing"),
+#' @include gmAnisotropy.R
+#' @include preparations.R
+#'  @export
+setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="missing"),
           function(object, newdata, pars, ...){
             if(is.null(object$pars)) object$pars = KrigingNeighbourhood()
             predict(object, newdata, pars = object$pars, ...)
@@ -418,12 +421,13 @@ setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="miss
 
 
 #' @rdname predict_gmSpatialModel
+#' @include gmSpatialMethodParameters.R
 #' @export
-setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmNeighbourhoodSpecification"),
+setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmNeighbourhoodSpecification"),
           function(object, newdata, pars, ...){
             cat("starting cokriging \n")
             object@parameters = pars
-            out = gstat:::predict.gstat(as.gstat(object), newdata=newdata)
+            out = predict(as.gstat(object), newdata=newdata, ...)
             return(out)
           }
 )
@@ -432,7 +436,7 @@ setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmNe
 
 #' @rdname predict_gmSpatialModel
 #' @export
-setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmTurningBands"),
+setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmTurningBands"),
           function(object, newdata, pars, ...){
             stop("Turning Bands method not yet interfaced here; use")
             cat("starting turning bands \n")
@@ -442,7 +446,7 @@ setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmTu
 
 #' @rdname predict_gmSpatialModel
 #' @export
-setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmCholeskyDecomposition"),
+setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmCholeskyDecomposition"),
           function(object, newdata, pars, ...){
             stop("Choleski decomposition method not yet implemented")
             cat("starting Choleski decomposition \n")
@@ -454,7 +458,7 @@ setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmCh
 
 #' @rdname predict_gmSpatialModel
 #' @export
-setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmSequentialSimulation"),
+setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmSequentialSimulation"),
           function(object, newdata, pars, ...){
             cat("starting SGs \n")
             object@parameters = pars$ng
@@ -475,7 +479,7 @@ setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmSe
 
 #' @rdname predict_gmSpatialModel
 #' @export
-setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmDirectSamplingParameters"),
+setMethod("Predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmDirectSamplingParameters"),
           function(object, newdata, pars, ...){
             
             cat("starting direct sampling \n")
@@ -533,7 +537,7 @@ setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmDi
             dt.nw = as(newdata,"SpatialGridDataFrame")@data
             dt.ti = as(object@model,"SpatialGridDataFrame")@data
             if(all(gt.nw@cellsize!=gt.ti@cellsize)) 
-              stop("predict.gmSpatialModel with gmDirectSamplingParameters: inferred grid topologies for newdata, conditioning data and model do not coincide")
+              stop("predict for gmSpatialModel with gmDirectSamplingParameters: inferred grid topologies for newdata, conditioning data and model do not coincide")
             if(is.null(mask)) mask = rep(TRUE, nrow(dt.nw))
             #erg = gsi.DS4CoDa(n=pars$patternSize, f=pars$scanFraction, t=pars$gof, n_realiz=pars$nsim, 
             #            nx_TI=gt.ti@cells.dim[1], ny_TI=gt.ti@cells.dim[2], 
@@ -559,6 +563,7 @@ setMethod("predict",signature(object="gmSpatialModel", newdata="ANY", pars="gmDi
 
 #' @describeIn gmSpatialModel Compute a variogram, see [variogram_gmSpatialModel()] for details
 #' @inheritParams variogram_gmSpatialModel
+#' @include variograms.R
 #' @export
 setMethod("variogram", signature=c(object="gmSpatialModel"), 
           def=variogram_gmSpatialModel)
@@ -569,6 +574,7 @@ setMethod("variogram", signature=c(object="gmSpatialModel"),
 #' objects
 #' @inheritParams logratioVariogram
 #' @inheritParams logratioVariogram_gmSpatialModel
+#' @include compositionsCompatibility.R
 #' @export
 setMethod("logratioVariogram", "gmSpatialModel", logratioVariogram_gmSpatialModel)
 
@@ -577,6 +583,7 @@ setMethod("logratioVariogram", "gmSpatialModel", logratioVariogram_gmSpatialMode
 #' for details
 #' @inheritParams as.gstat
 #' @export
+#' @include gstatCompatibility.R
 setMethod("as.gstat", signature="gmSpatialModel", def=as.gstat.gmSpatialModel)
 
 
diff --git a/R/gmValidationStrategy.R b/R/gmValidationStrategy.R
index 11cd8bb..256a2bf 100644
--- a/R/gmValidationStrategy.R
+++ b/R/gmValidationStrategy.R
@@ -45,9 +45,6 @@ LeaveOneOut = function(){
   return(res)
 }
 
-setOldClass("LeaveOneOut")
-setOldClass("NfoldCrossValidation")
-
 
 #### validation methods -------------
 #' Validate a spatial model
@@ -64,7 +61,7 @@ setOldClass("NfoldCrossValidation")
 #' @export
 #' @family validation functions 
 #' @family accuracy functions
-#' @name validate_gmSpatialModel
+#'
 #' @examples
 #' data("Windarling")
 #' X = Windarling[,c("Easting","Northing")]
@@ -77,41 +74,38 @@ setOldClass("NfoldCrossValidation")
 #' vs2 = NfoldCrossValidation(nfolds=sample(1:10, nrow(X), replace=TRUE))
 #' vs2
 #' \dontrun{ v2 = validate(gs, strategy=vs2) # quite slow }
-NULL
-
-#' @rdname validate_gmSpatialModel
-#' @export
-setGeneric("validate", function(object, strategy, ...){
-  standardGeneric("validate")
-})
+validate <- function(object, strategy, ...) UseMethod("validate", strategy)
 
 
-#' @rdname validate_gmSpatialModel
+#' @describeIn validate Validate a spatial model
+#' @method validate LeaveOneOut
 #' @export
-setMethod("validate",signature(object="gstat", strategy="LeaveOneOut"),
-       function(object, strategy, ...){
-  n = nrow(object$data[[1]]$data)
+validate.LeaveOneOut = function(object, strategy, ...){
+  if("gstat" %in% class(object)){
+    n = nrow(object$data[[1]]$data)
+  }else if(is(object, "gmSpatialModel")){
+    n = nrow(object@data)
+  }else{
+    object = try(as.gmSpatialModel(object))
+    if(class(object)=="try-error") stop("validate.LeaveOneOut: provided object not interpretable")
+    n = nrow(object@data)
+  }
   v = validate(object, NfoldCrossValidation(nfolds=n, doAll=TRUE))
   return(v)
-})
-
-#' @rdname validate_gmSpatialModel
-#' @export
-setMethod("validate",signature(object="gmSpatialModel", strategy="LeaveOneOut"),
-  function(object, strategy, ...){
-    n = nrow(object@data)
-    v = validate(object, NfoldCrossValidation(nfolds=n, doAll=TRUE))
-    return(v)
-})
-
-
-
+}
 
-#' @rdname validate_gmSpatialModel
+#' @describeIn validate Validate a spatial model
+#' @method validate NfoldCrossValidation
 #' @export
-setMethod("validate",signature(object="gmSpatialModel", strategy="NfoldCrossValidation"),
-          function(object, strategy, ...){
+validate.NfoldCrossValidation = function(object, strategy, ...){
+  # manage "gstat" case
+  if("gstat" %in% class(object)){
+    warning("validate: object provided is of class 'gstat', attempting 'gstat.cv(..., remove.all=TRUE, all.residuals=TRUE)'")
+    return(gstat::gstat.cv(object, nfold=strategy$nfolds, remove.all = TRUE, all.residuals = TRUE))
+  }
   # manage "gmSpatialModel" case
+  object = try(as.gmSpatialModel(object))
+  if(class(object)=="try-error") stop("validate.NfoldCrossValidation: provided object not interpretable")
   # interpret the information about the n-folds provided
   n = strategy$nfolds
   m = nrow(object@data)
@@ -144,19 +138,7 @@ setMethod("validate",signature(object="gmSpatialModel", strategy="NfoldCrossVali
   # reorder and return
   res = res[nms,] 
   return(res)
-})
-
-
-
-
+}
 
-#' @rdname validate_gmSpatialModel
-#' @export
-setMethod("validate",signature(object="gstat", strategy="NfoldCrossValidation"),
-          function(object, strategy, ...){
-            # manage "gstat" case
-              warning("validate: object provided is of class 'gstat', attempting 'gstat.cv(..., remove.all=TRUE, all.residuals=TRUE)'")
-              return(gstat::gstat.cv(object, nfold=strategy$nfolds, remove.all = TRUE, all.residuals = TRUE))
-          })
 
 
diff --git a/R/grids.R b/R/grids.R
index 8651541..0e78c3c 100644
--- a/R/grids.R
+++ b/R/grids.R
@@ -391,7 +391,7 @@ spatialGridRmult = function(coords, data, dimcomp=2, dimsim=NA){
 #' for "spatialGridAcomp" and "spatialGridRmult", but the default method is able to 
 #' deal with "SpatialPointsDataFrame", "SpatialPixelsDataFrame" and "SpatialGridDataFrame"
 #' objects, and with the "data.frame" output of [gstat::predict.gstat()] and 
-#' [predict.gmSpatialModel()]
+#' [predict_gmSpatialModel]
 #' @param ... generic functionality, currently ignored
 #'
 #' @return Invisibly, a list with elements `breaks` and `col` containing the breaks 
diff --git a/R/gstatCompatibility.R b/R/gstatCompatibility.R
index c7ab6ce..bcf7b7a 100644
--- a/R/gstatCompatibility.R
+++ b/R/gstatCompatibility.R
@@ -1,5 +1,5 @@
 #### gstat easy/easier interface for multivariate data
-setOldClass("gstat")
+# setOldClass("gstat") ### breaks down
 
 #' Fit an LMC to an empirical variogram
 #' 
diff --git a/R/preparations.R b/R/preparations.R
new file mode 100644
index 0000000..b64ccf5
--- /dev/null
+++ b/R/preparations.R
@@ -0,0 +1,47 @@
+### generic S4 functions --------------
+#' @export
+#' @rdname predict_gmSpatialModel
+setGeneric("Predict", function(object, newdata, pars, ...){
+  standardGeneric("Predict")
+})
+
+#' @export
+#' @rdname predict_gmSpatialModel
+setGeneric("predict", function(object,...) standardGeneric("predict"))  
+
+
+
+#' @export
+setGeneric("variogram", function(object,...) standardGeneric("variogram"))  
+
+
+# @include gmAnisotropy.R
+# @include gmValidationStrategy.R 
+# @include variograms.R
+setOldClass("gmKrigingNeighbourhood")
+setOldClass("gmDirectSamplingParameters")
+setOldClass("gmTurningBands")
+setOldClass("gmSequentialSimulation")
+setOldClass("gmCholeskyDecomposition")
+setOldClass("LeaveOneOut")
+setOldClass("NfoldCrossValidation")
+
+
+
+# @include gstatCompatibility.R
+# @include compositionsCompatibility.R
+setOldClass("gmCgram")
+setOldClass("LMCAnisCompo")
+setOldClass("variogramModelList")
+setOldClass("variogramModel")
+setOldClass("genDiag")
+
+
+# @include gstatCompatibility.R
+# @include compositionsCompatibility.R
+setOldClass("gmEVario")
+setOldClass("logratioVariogram")
+setOldClass("logratioVariogramAnisotropy")
+setOldClass("gstatVariogram")
+
+
diff --git a/R/variograms.R b/R/variograms.R
index 0026e92..5de2167 100644
--- a/R/variograms.R
+++ b/R/variograms.R
@@ -272,6 +272,7 @@ as.function.gmCgram = function(x,...){
 #' @describeIn as.function.gmCgram predict a gmCgram object on some lag vector coordinates
 #' @param object gmCgram object
 #' @param newdata matrix, data.frame or Spatial object containing coordinates
+#' @include gmAnisotropy.R
 #' @method predict gmCgram
 #' @export
 predict.gmCgram = function(object, newdata, ...){
@@ -553,9 +554,6 @@ variogram_gmSpatialModel <-  function(object, methodPars=NULL, ...){
 #  return(variogram_gmSpatialModel(object, ...))
 #}
 
-  
-#' @export
-setGeneric("variogram", function(object,...) standardGeneric("variogram"))  
 
 
 
@@ -1066,17 +1064,15 @@ gsi.midValues.azimuthInterval <- function(x){ (x[[1]]+x[[2]])/2 }
 ## theoretical structural functions
 # S3 -> S4 classes  
 # cat("creating variogram model classes\n")
-#' @include gstatCompatibility.R
-#' @include compositionsCompatibility.R
-setOldClass("gmCgram")
-setOldClass("LMCAnisCompo")
-setOldClass("variogramModelList")
-setOldClass("variogramModel")
+
 
 # abstract classes
 #' @title Structural function model specification
 #' @description Abstract class, containing any specification of a variogram (or covariance) model
 #' @export
+#' @include compositionsCompatibility.R
+#' @include gstatCompatibility.R
+#' @include preparations.R
 setClassUnion(name="ModelStructuralFunctionSpecification", 
               members=c("NULL","gmCgram", "LMCAnisCompo", "variogramModelList", "variogramModel"))
 
@@ -1143,15 +1139,14 @@ setClass("gmGaussianModel",
 ## empirical structural functions
 # S3 -> S4 classes
 # cat("creating empirical variogram classes\n")
-setOldClass("gmEVario")
-setOldClass("logratioVariogram")
-setOldClass("logratioVariogramAnisotropy")
-setOldClass("gstatVariogram")
 
 
 # abstract classes
 #' @title Empirical structural function specification
 #' @description Abstract class, containing any specification of an empirical variogram (or covariance function, or variations)
 #' @export
+#' @include compositionsCompatibility.R
+#' @include gstatCompatibility.R
+#' @include preparations.R
 setClassUnion(name="EmpiricalStructuralFunctionSpecification", members=c("NULL","gmEVario", "logratioVariogram", "logratioVariogramAnisotropy", "gstatVariogram"))
 
diff --git a/src/gmGeostats.c b/src/gmGeostats.c
index e003d2c..4ad5886 100644
--- a/src/gmGeostats.c
+++ b/src/gmGeostats.c
@@ -103,7 +103,7 @@ vgramFunctionPtr cgramFunctions[]={
 
 static R_NativePrimitiveArgType CcalcCgram_t[] = {
   /* dimX  LDX   X       dimY    LDY    Y       dimC   C     Nugget  nCgr  typeCgr  A       Sill    moreC   ijEq*/
-  INTSXP,INTSXP,REALSXP,INTSXP,INTSXP,REALSXP,INTSXP,REALSXP,REALSXP,INTSXP,REALSXP,REALSXP,REALSXP,REALSXP,INTSXP
+  INTSXP,INTSXP,REALSXP,INTSXP,INTSXP,REALSXP,INTSXP,REALSXP,REALSXP,INTSXP,INTSXP,REALSXP,REALSXP,REALSXP,LGLSXP
 };
 
 void CcalcCgram(
@@ -878,14 +878,14 @@ static R_CMethodDef cMethods[] = {
 //{"CMVTurningBands2", (DL_FUNC) &CMVTurningBands2, 11, CMVTurningBands2_t},
 
 
-/*
+
 void R_init_gmGeostats(DllInfo *info)
 {
   R_registerRoutines(info, cMethods, NULL, NULL, NULL);
   R_useDynamicSymbols(info, FALSE);
   R_forceSymbols(info, TRUE);
 }
-*/
+
 
 
 
-- 
GitLab