.Rbuildignore

@@ -14,3 +14,4 @@ svn
 jurapred
 ^README\.Rmd$
 ^.*\.vscode$
+^revdep$

.gitignore

@@ -4,4 +4,5 @@ inst/doc
 1_oldcode
 R-future/
 tests0/
-experiments/
+experiments/
\ No newline at end of file
+.Rbuildignore

DESCRIPTION

 Package: gmGeostats
-Version: 0.11.0-9002
-Date: 2021-12-14
+Version: 0.11.1
+Date: 2022-05-03
 Title: Geostatistics for Compositional Analysis
 Authors@R: c(person(given = "Raimon", 
         family = "Tolosana-Delgado", 

NEWS.md

+# gmGeostats 0.11.1
+
+* (2022-05-03) minor changes to adapt to more stringent standards for class determination, for variable definition in C code, and for calls to FORTRAN routines from C code involving string arguments.
+
 # gmGeostats 0.11.0-9002
 
 * (2021-12-14) bugs in turning bands corrected: getUnitVec was producing a wrong sequence of directions, and bands for exponential and Gaussian variograms did not correct for the difference between parametric and effective range in the right way; 

R/compositionsCompatibility.R

@@ -61,7 +61,7 @@ logratioVariogram_acomp <- function(data=comp, loc, ..., azimuth=0, azimuth.tol=
 logratioVariogram_gmSpatialModel <- function(data, ..., azimuth=0, azimuth.tol=180/length(azimuth)){
   coords = sp::coordinates(data)
   comp = tryCatch(gsi.orig(data@data))
-  if(class(comp)=="try-catch") stop("logratioVariogram.gmSpatialModel: provided data is not compositional!")
+  if(inherits(comp,"try-catch")) stop("logratioVariogram.gmSpatialModel: provided data is not compositional!")
   if(length(azimuth)>1) return(gsi.logratioVariogramAnisotropy(comp, coords, ..., azimuths=azimuth, azimuth.tol=azimuth.tol))
   logratioVariogram.default(comp, loc=coords, ...)    
 }

R/gmSimulation.R

@@ -419,7 +419,7 @@ gsi.DS <- function(n, f, t, n_realiz,
   }else if(is.character(ivars_TI)){
      varnames_out = ivars_TI
   }
-  if(length(varnames_out)!=D | class(varnames_out)=="try-error") varnames_out = paste("v", 1:D, sep="")
+  if(length(varnames_out)!=D | inherits(varnames_out,"try-error")) varnames_out = paste("v", 1:D, sep="")
   dm = list(loc=1:length(mask), var=varnames_out, sim=paste("sim", 1:n_realiz, sep="") )
   
   SimGrid_list = lapply(SimGrid_list, function(x){

R/gmSpatialMethodParameters.R

@@ -36,7 +36,7 @@
 KrigingNeighbourhood <- function(nmax=Inf, nmin=0, omax=0, maxdist=Inf, force=FALSE, anisotropy=NULL, ...){
   if(!is.null(anisotropy)){
     anisotropy = try(as.AnisotropyScaling(anisotropy))
-    if(class(anisotropy)=="try-error") stop("krigingNeighbourhood: anisotropy description provided cannot be parsed")
+    if(inherits(anisotropy,"try-error")) stop("krigingNeighbourhood: anisotropy description provided cannot be parsed")
   }
   # here space for checking that parameters are sensible
   
@@ -61,6 +61,8 @@ KrigingNeighbourhood <- function(nmax=Inf, nmin=0, omax=0, maxdist=Inf, force=FA
 #' @param scanFraction maximum fraction of the training image to be scanned on each iteration
 #' @param patternSize number of observations used for conditioning the simulation
 #' @param gof maximum acceptance discrepance between a data event in the training image and the conditioning data event
+#' @param seed an object specifying if and how the random number generator should be 
+#' initialized, see `?simulate` in base "stats" package
 #' @param ... further parameters, not used
 #'
 #' @return an S3-list of class "gmDirectSamplingParameters" containing the six elements given as arguments 
@@ -73,7 +75,7 @@ KrigingNeighbourhood <- function(nmax=Inf, nmin=0, omax=0, maxdist=Inf, force=FA
 #' @examples
 #' (dsp = DSpars(nsim=100, scanFraction=75, patternSize=6, gof=0.05))
 #' ## then run predict(..., pars=dsp)
-DSpars <- DirectSamplingParameters <- function(nsim=1, scanFraction=0.25, patternSize=10, gof=0.05, ...){
+DSpars <- DirectSamplingParameters <- function(nsim=1, scanFraction=0.25, patternSize=10, gof=0.05, seed=NULL, ...){
   ll = list(nsim=nsim, scanFraction=scanFraction, patternSize=patternSize, gof=gof, ...)
   class(ll) = "gmDirectSamplingParameters"
   return(ll)
@@ -93,6 +95,8 @@ DSpars <- DirectSamplingParameters <- function(nsim=1, scanFraction=0.25, patter
 #' @param rank currently ignored (future functionality: obtain a reduced-rank simulation)
 #' @param debug.level degree of verbosity of results; negative values produce a progress bar; values can be
 #' extracted from [gstat::predict.gstat()] 
+#' @param seed an object specifying if and how the random number generator should be 
+#' initialized, see `?simulate` in base "stats" package
 #' @param ... further parameters, currently ignored
 #'
 #' @return an S3-list of class "gmSequentialSimulation" containing the four elements given as arguments 
@@ -109,7 +113,7 @@ DSpars <- DirectSamplingParameters <- function(nsim=1, scanFraction=0.25, patter
 #' ng_local = KrigingNeighbourhood(maxdist=1, nmin=4, omax=5, force=TRUE)
 #' (sgs_local = SequentialSimulation(nsim=100, ng=ng_local, debug.level=-1))
 #' ## then run predict(..., pars=sgs_local)
-SequentialSimulation = function(nsim=1, ng=NULL, rank=Inf, debug.level=1, ...){
+SequentialSimulation = function(nsim=1, ng=NULL, rank=Inf, debug.level=1, seed=NULL, ...){
   if(is.null(ng)) warning("SequentialSimulation: local neighbourhood is required; calculations will be stopped if the spatial model object does not include it")
   res = list(nsim=nsim, ng=ng, rank=rank, debug.level=debug.level, ...)
   class(res) = "gmSequentialSimulation"
@@ -125,6 +129,8 @@ SequentialSimulation = function(nsim=1, ng=NULL, rank=Inf, debug.level=1, ...){
 #' 
 #' @param nsim number of realisations desired
 #' @param nBands number of bands desired for the decomposition of the 2D or 3D space in individual signals 
+#' @param seed an object specifying if and how the random number generator should be 
+#' initialized, see `?simulate` in base "stats" package
 #' @param ... further parameters, currently ignored
 #'
 #' @return an S3-list of class "gmTurningBands" containing the few elements given as arguments 
@@ -136,7 +142,7 @@ SequentialSimulation = function(nsim=1, ng=NULL, rank=Inf, debug.level=1, ...){
 #' @examples
 #' (tbs_local = TurningBands(nsim=100, nBands=300))
 #' ## then run predict(..., pars=tbs_local)
-TurningBands = function(nsim=1, nBands=1000,  ...){
+TurningBands = function(nsim=1, nBands=1000, seed=NULL, ...){
   res = list(nsim=nsim, nBands=nBands, ...)
   class(res) = "gmTurningBands"
   return(res)
@@ -149,6 +155,8 @@ TurningBands = function(nsim=1, nBands=1000,  ...){
 #' mostly for covariance or variogram-based gaussian random fields.
 #' 
 #' @param nsim number of realisations desired
+#' @param seed an object specifying if and how the random number generator should be 
+#' initialized, see `?simulate` in base "stats" package
 #' @param ... further parameters, currently ignored
 #'
 #' @return an S3-list of class "gmCholeskyDecomposition" containing the few elements given as arguments 
@@ -160,7 +168,7 @@ TurningBands = function(nsim=1, nBands=1000,  ...){
 #' @examples
 #' (chols_local = CholeskyDecomposition(nsim=100, nBands=300))
 #' ## then run predict(..., pars=chols_local)
-CholeskyDecomposition = function(nsim=1, ...){
+CholeskyDecomposition = function(nsim=1, seed=NULL, ...){
   res = list(nsim=nsim,  ...)
   class(res) = "gmCholeskyDecomposition"
   return(res)

R/gmSpatialModel.R

@@ -369,6 +369,7 @@ as.gmSpatialModel.gstat = function(object, V=NULL, ...){
   stop("as.gmSpatialModel.gstat: not yet implemented")
 }
 
+## predict and simulate methods -------------
 
 #' Predict method for objects of class 'gmSpatialModel'
 #' 
@@ -682,4 +683,3 @@ setMethod("as.gstat", signature="gmSpatialModel", def=as.gstat.gmSpatialModel)
 #      }
 #)
 
-

R/gmValidationStrategy.R

@@ -87,7 +87,7 @@ validate.LeaveOneOut = function(object, strategy, ...){
     n = nrow(object@data)
   }else{
     object = try(as.gmSpatialModel(object))
-    if(class(object)=="try-error") stop("validate.LeaveOneOut: provided object not interpretable")
+    if(inherits(object,"try-error")) stop("validate.LeaveOneOut: provided object not interpretable")
     n = nrow(object@data)
   }
   v = validate(object, NfoldCrossValidation(nfolds=n, doAll=TRUE))
@@ -105,7 +105,7 @@ validate.NfoldCrossValidation = function(object, strategy, ...){
   }
   # manage "gmSpatialModel" case
   object = try(as.gmSpatialModel(object))
-  if(class(object)=="try-error") stop("validate.NfoldCrossValidation: provided object not interpretable")
+  if(inherits(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)

R/grids.R

@@ -527,7 +527,7 @@ image_cokriged.spatialGridRmult<- function(x, ivar=1, isim=NULL, breaks=10, mask
     }else{
       X = x[,isim, ivar]
     }
-    if(!is.null(mask) & class(mask)=="mask"){
+    if(!is.null(mask) & inherits(mask,"mask")){
       if(!is.null(attr(mask, "fullgrid"))){
         X = unmask(data.frame(X), mask=mask)[,1]
         coords = sp::coordinates(attr(mask, "fullgrid"))

R/gstatCompatibility.R

@@ -256,7 +256,7 @@ variogramModelPlot.gstatVariogram =
         vrnames = levels(model$id)[1]
       }else{
         ggtr = tryCatch(as.variogramModel(model))
-        if(class(ggtr)=="try-error"){
+        if(inherits(ggtr,"try-error")){
           stop("argument 'gg' must either be a variogramModel, a gstat object with a variogramModel, a variogramModelList object, or an object convertible to one") 
         }else{
           return(variogramModelPlot(vg=vg,  gg = ggtr, col = col, commonAxis = commonAxis, newfig = newfig, ...))

R/mask.R

@@ -80,12 +80,12 @@ constructMask = function(grid, method="maxdist", maxval=NULL, x=NULL){
       x = data.frame(sp::coordinates(x), x@data)
     }
     x = try(as.data.frame(x))
-    if(class(x)=="try-error") stop("constructMask: provided object x should be a data.frame or convertible to it for method 'maxdist'")
+    if(inherits(x,"try-error")) stop("constructMask: provided object x should be a data.frame or convertible to it for method 'maxdist'")
     out = gsi.masking.nearest(grid0, x, maxdist=maxval)
   }else if(m=="sillprop"){
     if(is.null(x))
       stop("constructMask: sillprop method requires a variogram model")
-    if(class(x)=="gstat") x = x$model
+    if(inherits(x,"gstat")) x = x$model
     if(is(x,"gmSpatialModel")) x = x@model@structure
     if(is(x, "ModelStructuralFunctionSpecification")) as.variogramModel(x)
     maxval = ifelse(is.null(maxval), 0.99, maxval)
@@ -166,7 +166,7 @@ gsi.masking.nearest = function(grid, x, maxdist){
 gsi.masking.polygon = function(grid, poly){
   requireNamespace("sp", quietly = TRUE)
   poly = try(as(poly, "SpatialPolygons"))
-  if(class(poly)=="try-error")
+  if(inherits(poly,"try-error"))
     stop("object 'poly' cannot be coerced to SpatialPolygons")
   FUN = function(i){
     poly = poly@polygons[[i]]@Polygons[[1]]@coords
@@ -260,7 +260,7 @@ setMask <- function(x,...) UseMethod("setMask", x)
 #' to specify the matrix of spatial coordinates (all `setMask` methods including it)
 setMask.default <- function(x, mask, coordinates = 1:2, ...){
   x = as.data.frame(x)
-  if(class(mask)=="mask") attributes(mask) = NULL
+  if(inherits(mask,"mask")) attributes(mask) = NULL
   if(is.null(dim(coordinates) )){
     fullgrid = x[,coordinates]
   }else{
@@ -281,7 +281,7 @@ setMask.data.frame <- setMask.default
 #' @method setMask DataFrameStack
 #' @export
 setMask.DataFrameStack <- function(x, mask, coordinates=attr(x, "coordinates"), ...){
-  if(class(mask)=="mask") attributes(mask) = NULL
+  if(inherits(mask,"mask")) attributes(mask) = NULL
   cc = coordinates
   x = x[mask,,drop=FALSE]
   attr(mask, "fullgrid") = cc
@@ -298,7 +298,7 @@ setMask.SpatialGrid <- function(x, mask, ...){
   r = order(+cc[,2],+cc[,1])
   o = 1:nrow(cc)
   r = o[r]
-  if(class(mask)=="mask") attributes(mask) = NULL
+  if(inherits(mask,"mask")) attributes(mask) = NULL
   maskaux = mask[o]
   cc = cc[maskaux,, drop=FALSE]
   cc = sp::SpatialPoints(coords = cc, proj4string = sp::CRS(sp::proj4string(x)), 
@@ -331,7 +331,7 @@ setMask.GridTopology <- function(x, mask, ...){
 #' @importClassesFrom sp SpatialPoints
 setMask.SpatialPoints <- function(x, mask, ...){
   cc = sp::coordinates(x)
-  if(class(mask)=="mask") attributes(mask) = NULL
+  if(inherits(mask,"mask")) attributes(mask) = NULL
   cc = cc[mask,,drop=FALSE]
   if("data" %in% slotNames(x)){
     dt = x@data
@@ -433,7 +433,7 @@ unmask.SpatialPixels <- function(x, mask=NULL, fullgrid =attr(mask, "fullgrid"),
   if(is(fullgrid, "SpatialGrid")) fullgrid = sp::getGridTopology(fullgrid)
   # compute number of points
   npoints <- try( prod(fullgrid@cells.dim))
-  if(class(npoints)=="try-error") stop("unmask.SpatialPixels: provided fullgrid could not be intepreted as a grid")
+  if(inherits(npoints,"try-error")) stop("unmask.SpatialPixels: provided fullgrid could not be intepreted as a grid")
   # construct mask
   if(is.null(mask)){
     mask = rep(FALSE, npoints)

cran-comments.md

 ## Test environments
   
 * Linux: Ubuntu 18.04.4LTS  (local)
-  - R 3.4.4 (production)
-  - R devel 2020-07-28 r78930 (test)
-  - R patched 4.0.2 r78930 (test)
-* Windows 7 Enterprise (test)
-  - R 3.5.1 (test)
-  - R 4.0.2 (test)
-  
+  - R 3.6.3 (production)
+  - R devel 2022-04-23 r82240 (test)
+
   
 ## R CMD check results
   
@@ -18,23 +14,19 @@
   
 ```
 Dear maintainer,
-  
+
 Please see the problems shown on
-  <https://cran.r-project.org/web/checks/check_results_gmGeostats.html>.
-  
-  Please correct before 2020-10-03 to safely retain your package on CRAN.
-  
+<https://www.stats.ox.ac.uk/pub/bdr/BLAS/gmGeostats.out>
+
+For reproduction details see
+<https://www.stats.ox.ac.uk/pub/bdr/BLAS/README.txt> .
+
+Please correct before 2022-05-05 to safely retain your package on CRAN.
+
 The CRAN Team
 ```
   
-... the sympthom being an error while #include <omp.h> in r-release-macos-x86_64, we identified the problem as a lack of appropriate linkage of OpenMP. 
-  
-Solution: 
-  1.- the header inclusion and all usages of #pragma directives were enclosed in conditional structures
-#ifdef _OPENMP 
-   ...
-#endif 
-  
-   2. `$(SHLIB_OPENMP_CFLAGS)` was added to PKG_CFLAGS and PKG_LIBS fields in Makevars files
+... related to the way classes were assessed (now using `inherits(...)`) and the calls of FORTRAN routines from C code that involve string arguments (now following the instructions given in <https://www.stats.ox.ac.uk/pub/bdr/BLAS/README.txt> ).
   
-  NOTE: we do not have access to any MacOS machine to test that our solution is effective. Our apologies if this is not working!
+
\ No newline at end of file
+best regards
\ No newline at end of file

src/gmGeostats.c

@@ -9,6 +9,9 @@
 #ifdef _OPENMP
 #include <omp.h>
 #endif
+
+#define USE_FC_LEN_T
+
 #include <Rinternals.h>
 
 #define inR   // attention: this must be uncommented if not compiling
@@ -20,6 +23,10 @@
 #include <R_ext/Lapack.h> 
 #endif
 
+#ifndef FCONE
+# define FCONE
+#endif
+
 #define maxIntervals 1000
 short int binBuf[maxIntervals];
 double doubleBuf[maxIntervals];
@@ -140,7 +147,7 @@ const int *ijEqual
   if( m<1 || m>3)
     error("Can not handel spatial dimensions outside 1-3");
   int outBufSize=d*d*nX*nY;
-  int i,j,k,ev,lx,ly,s;
+  int i,j,k,lx,ly,s; /* 20220424: ev removed */
   double delta[3];
   double v[3];
   double h2,h,val;
@@ -328,7 +335,7 @@ void getUnitvec(
 ) {
   /* weak discrepancy sequence of pseudorandom directions in 2D or 3D:
    * Lantuejoul (2002), page 194, after Freulon (1992) */			  
-  int i;
+  /* 20220424: removed int i; */ 
   double d1,d2,d3;
   if(dimX>3)
     error("no expression for unit vectors in dimension larger than 3");
@@ -387,9 +394,9 @@ void CMVTurningBands(
   const int maxCgramType=2;
   const int nsim=dimZ[2];
   int i,j,k,s,ss,ev;
-  double d1,d2,d3;
+  double d1,d2; /* 20220424: removed ,d3; */
   const double sqrtNBands=sqrt((double) *nBands);
-  double phase,amp;
+  double phase; /* 20220424: removed ,amp; */
   const int n=dimX[0];
   const int m=dimX[1];
   const int d=dimZ[0];
@@ -412,50 +419,50 @@ void CMVTurningBands(
       for(j=0;j<d;j++)
         Z[d*i+j]=0.0;
     for(s=0;s<*nBands;s++){/* band */
-  getUnitvec(3, s+1, &(omega[0])); /* obtain a direction; always in 3D, in order for the spherical variogram to be correct */
-  //getUnitvec(m, s+1, omega); /* obtain a direction */
-  for(ss=0;ss< *nCgrams;ss++) { /* variogram structure */
-  if( typeCgram[ss]<0 || typeCgram[ss] > maxCgramType )
-    error("CMVTurningBands: Unknown variogram type");
-  /* project all data onto the direction */
-  for(i=0;i<n;i++){ /* location */
-  for(j=0;j<m;j++) {
-    v[j]=0;
-    for(k=0;k<m;k++)
-      v[j]+=A[ss+ *nCgrams *(j+m*k)]*X[i+n*k];
-  }
-  projs[i]=0;
-    for(j=0;j<m;j++){ /* spatial dimension */
-  projs[i]+=omega[j]*v[j];
-    }
-  }
-  /* for each eigenvalue, ... */
-  for(ev=0;ev<d;ev++) { /* eigenvector */
-  /* ... obtain a curve at all proj points following the covariance model */  
-  (*bandSim[typeCgram[ss]])(n,projs,band,1.0,moreCgramData+ss);  
-    /* this function takes the projs and returns on band the  the curve */
-    /* ... multiply the eigenvector by the curve, and accumulate */
+      getUnitvec(3, s+1, &(omega[0])); /* obtain a direction; always in 3D, in order for the spherical variogram to be correct */
+      //getUnitvec(m, s+1, omega); /* obtain a direction */
+      for(ss=0;ss< *nCgrams;ss++) { /* variogram structure */
+        if( typeCgram[ss]<0 || typeCgram[ss] > maxCgramType )
+           error("CMVTurningBands: Unknown variogram type");
+        /* project all data onto the direction */
+        for(i=0;i<n;i++){ /* location */
+         for(j=0;j<m;j++) {
+          v[j]=0;
+          for(k=0;k<m;k++)
+           v[j]+=A[ss+ *nCgrams *(j+m*k)]*X[i+n*k];
+         }
+         projs[i]=0;
+         for(j=0;j<m;j++){ /* spatial dimension */
+          projs[i]+=omega[j]*v[j];
+         }
+        }
+        /* for each eigenvalue, ... */
+        for(ev=0;ev<d;ev++) { /* eigenvector */
+         /* ... obtain a curve at all proj points following the covariance model */  
+         (*bandSim[typeCgram[ss]])(n,projs,band,1.0,moreCgramData+ss);  
+         /* this function takes the projs and returns on band the  the curve */
+         /* ... multiply the eigenvector by the curve, and accumulate */
 #ifdef _OPENMP
 #pragma omp parallel for		             \
     if(!omp_in_parallel()&&0)		        \
       num_threads(omp_get_num_procs())	\
       default(shared) private(i,j,d2) 
-    for(i=0;i<n;i++){ /* location */
-    for(j=0;j<d;j++){ /* variable */
-    d2 = sqrtSill[ss + *nCgrams *(j+d*ev)]*band[i]; 
-      Z[d*i+j]+=d2;  
-    }	  
-    }
+         for(i=0;i<n;i++){ /* location */
+          for(j=0;j<d;j++){ /* variable */
+           d2 = sqrtSill[ss + *nCgrams *(j+d*ev)]*band[i]; 
+           Z[d*i+j]+=d2;  
+          }	  
+         }
 #else
-    for(i=0;i<n;i++){ /* location */
-    for(j=0;j<d;j++){ /* variable */
-    d2 = sqrtSill[ss + *nCgrams *(j+d*ev)]*band[i]; 
-      Z[d*i+j]+=d2;  
-    }	  
-    }
+         for(i=0;i<n;i++){ /* location */
+          for(j=0;j<d;j++){ /* variable */
+           d2 = sqrtSill[ss + *nCgrams *(j+d*ev)]*band[i]; 
+           Z[d*i+j]+=d2;  
+          }	  
+         }
 #endif
-  }
-  }
+        }
+      }
     } 
     /* Rescale*/
     for(i=0;i<n;i++)
@@ -505,7 +512,7 @@ void CCondSim(
     double *cbuf,  /* BUF: Buffer of length d*d*nin */
     double *dbuf  /* BUF: Buffer of length d*nin*nsim */
 ) {
-  const int maxCgramType=2;
+  /* 20220424 removed: , const int maxCgramType=2; */
   const int n=dimX[0];
   const int m=dimX[1];
   const int d=dimZin[0];
@@ -520,9 +527,9 @@ void CCondSim(
   const double zero=0.0;
   const double one=1.0;
   const double minus1=-1.0;
-  int i,j,k,s,ss,ev,l;
+  int i,j,k; /* 20220424 removed: ,s,ss,ev,l; */
   int sim,shift;
-  double d1,d2,d3,cv;
+  /* 20220424 removed: double d1,d2,d3,cv; */
   const int dimXin[2] = {nin,m}; 
   const int dimXout[2] = {1,m};
   const int dimCbuf[4] = {d,nin,d,1};
@@ -593,7 +600,7 @@ void CCondSim(
                &oneI,
                &zero,
                dbuf+dmnin*sim,
-               &oneI);
+               &oneI FCONE);
     }
 #else
     for(sim=0;sim<nsim;sim++) {
@@ -607,7 +614,7 @@ void CCondSim(
                &oneI,
                &zero,
                dbuf+dmnin*sim,
-               &oneI);
+               &oneI FCONE);
     }
 #endif
     // /* points in*/
@@ -672,7 +679,7 @@ void CCondSim(
                    &oneI,
                    &one,
                    Z+i*d+nd*sim,
-                   &oneI);
+                   &oneI FCONE);
         }
 #else
         for(sim=0;sim<nsim;sim++) {
@@ -686,7 +693,7 @@ void CCondSim(
                    &oneI,
                    &one,
                    Z+i*d+nd*sim,
-                   &oneI);
+                   &oneI FCONE);
         }
 #endif
     }