From 67a78c277c4d827199b63260297f45ed13162f01 Mon Sep 17 00:00:00 2001
From: Raimon Tolosana-Delgado <tolosa53@fwg206.ad.fz-rossendorf.de>
Date: Tue, 3 May 2022 14:54:14 +0200
Subject: [PATCH] class determination now with is(...) or inherits(...); calls
to FORTRAN within C with string arguments platform independent
---
.Rbuildignore | 1 +
R/compositionsCompatibility.R | 2 +-
R/gmSimulation.R | 2 +-
R/gmSpatialMethodParameters.R | 18 ++++--
R/gmSpatialModel.R | 2 +-
R/gmValidationStrategy.R | 4 +-
R/grids.R | 2 +-
R/gstatCompatibility.R | 2 +-
R/mask.R | 14 ++---
src/gmGeostats.c | 103 ++++++++++++++++++----------------
10 files changed, 83 insertions(+), 67 deletions(-)
diff --git a/.Rbuildignore b/.Rbuildignore
index f0a1a18..f083d49 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -14,3 +14,4 @@ svn
jurapred
^README\.Rmd$
^.*\.vscode$
+^revdep$
diff --git a/R/compositionsCompatibility.R b/R/compositionsCompatibility.R
index b8fd531..0846f3b 100644
--- a/R/compositionsCompatibility.R
+++ b/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, ...)
}
diff --git a/R/gmSimulation.R b/R/gmSimulation.R
index a62837f..6762513 100644
--- a/R/gmSimulation.R
+++ b/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){
diff --git a/R/gmSpatialMethodParameters.R b/R/gmSpatialMethodParameters.R
index 5934247..6f5ec7f 100644
--- a/R/gmSpatialMethodParameters.R
+++ b/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)
diff --git a/R/gmSpatialModel.R b/R/gmSpatialModel.R
index 66255e0..b11abdc 100644
--- a/R/gmSpatialModel.R
+++ b/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)
# }
#)
-
diff --git a/R/gmValidationStrategy.R b/R/gmValidationStrategy.R
index 256a2bf..2db7e55 100644
--- a/R/gmValidationStrategy.R
+++ b/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)
diff --git a/R/grids.R b/R/grids.R
index 0e78c3c..c5aa3c3 100644
--- a/R/grids.R
+++ b/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"))
diff --git a/R/gstatCompatibility.R b/R/gstatCompatibility.R
index 725e19a..cb9b9fc 100644
--- a/R/gstatCompatibility.R
+++ b/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, ...))
diff --git a/R/mask.R b/R/mask.R
index 5168516..7bfe084 100644
--- a/R/mask.R
+++ b/R/mask.R
@@ -80,7 +80,7 @@ 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))
@@ -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)
diff --git a/src/gmGeostats.c b/src/gmGeostats.c
index 7a90c72..5d6b720 100644
--- a/src/gmGeostats.c
+++ b/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
}
--
GitLab