Second order stationary random functions

Selection of the covariance model. It can only be: - GAUSSIAN - EXPONENTIAL - MATERN

idx_type = 3
type_cov <- c("GAUSSIAN", "EXPONENTIAL", "MATERN")[idx_type]
nu <- c(10, 1/2, 3/4)[idx_type]
ranges = c(5, 10, 15)
angles = c(30., 0, 0)

# initialization of the simulator
nb <- 1000 # number of spectral components
ns <- 20   # number of simulations
seed <- 13112 # Seed for the first simulation

# Define the title (used for plots)
title = paste0(type_cov, " model")
if (type_cov == "MATERN") {title = paste0(title, " (nu = ", nu, ")")}

1D simulation of a random function

ndim = 1
# definition of the Space and creation of the grid
err = defineDefaultSpace(ESpaceType_RN(), ndim)
nx   = c(10000)
dx   = c(1.)
grid = DbGrid_create(nx = nx, dx = dx)

# specification of the model using Model
mod  = Model_createFromParam(type = ECov_fromKey(type_cov), 
                             ranges = ranges[1:ndim], param = nu, flagRange=FALSE)
# using the *SimuSpectral* interface
sim = SimuSpectralRN(mod$getCov())
err = sim$simulate(ns = nb, seed = 123, verbose = TRUE)
## Simulation of the spectrum
## - Space dimension   = 1
## - Number of variables  = 1
## - Number of spectral components = 1000
## Simulation of the spectrum
## - Space dimension   = 1
## - Number of variables  = 1
## - Number of spectral components = 1000
# 
gamma = sim$getGamma()$toTL()
omega = sim$getOmega()$toTL()
phi   = sim$getPhi()
coor  = grid$getColumnsAsMatrix(names = c("x1"), useSel = FALSE)$toTL()
val   = as.numeric(t(gamma) %*% cos(omega %*% t(coor) + phi))

# # using the interface
err = sim$compute(dbout = grid, verbose = TRUE, namconv = NamingConvention(prefix))
## Spectral Simulation on a set of Isolated Points
## - Number of samples = 10000
## - Space dimension   = 1
## - Number of variables = 1
## Spectral Simulation on a set of Isolated Points
## - Number of samples = 10000
## - Space dimension   = 1
## - Number of spectral components = 1000
## - Number of variables = 1
# 
stopifnot (all(abs(val - grid[paste(prefix, "V1", "simu", sep = ".")]) < 1e-6))

# using the simuSpectral interface
err  = simuSpectral(NULL, dbout = grid, cova = mod$getCov(), nbsimu = ns, ns = nb, seed=seed, cov0 = NULL, namconv = NamingConvention(prefix))
# statistics
knitr::kable(dbStatisticsMono(db = grid, names = paste(prefix, "*", sep = "."), opers = opers)$toTL(),
             digits = 3, caption = title)
MATERN model (nu = 0.75)
Number Minimum Maximum Mean St. Dev.
Simu.V1.simu 10000 -4.368 3.657 0.047 0.993
Simu.V1.S1 10000 -3.666 3.998 -0.049 1.007
Simu.V1.S2 10000 -4.020 3.618 -0.033 1.018
Simu.V1.S3 10000 -3.287 4.378 0.049 1.016
Simu.V1.S4 10000 -3.657 3.930 0.055 0.984
Simu.V1.S5 10000 -3.737 3.749 -0.014 1.061
Simu.V1.S6 10000 -3.425 3.458 0.063 0.980
Simu.V1.S7 10000 -3.489 3.426 -0.022 0.991
Simu.V1.S8 10000 -4.073 4.776 0.014 0.993
Simu.V1.S9 10000 -3.688 3.582 0.030 1.027
Simu.V1.S10 10000 -3.838 3.235 -0.101 1.030
Simu.V1.S11 10000 -3.977 3.918 -0.074 0.984
Simu.V1.S12 10000 -3.377 4.153 0.020 0.996
Simu.V1.S13 10000 -3.744 3.302 0.001 0.984
Simu.V1.S14 10000 -3.328 3.415 0.002 1.020
Simu.V1.S15 10000 -3.367 3.403 -0.083 0.969
Simu.V1.S16 10000 -3.838 3.474 -0.005 0.987
Simu.V1.S17 10000 -3.251 3.353 0.029 1.018
Simu.V1.S18 10000 -3.188 3.374 -0.008 0.996
Simu.V1.S19 10000 -4.585 3.880 0.007 1.028
Simu.V1.S20 10000 -5.182 3.945 -0.065 1.020 
# histogram
p = plot.init() + 
  plot.hist(db = grid, name = var_name, bins = 25, col = "skyblue") +
  plot.decoration(title = title)
plot.end(p)

# experimental variogram for the first simulation and simulated model
nlag = 50
varioPar = VarioParam_createMultipleFromGrid(grid, nlag = nlag)
err = grid$setLocator(name = var_name, locatorType = ELoc_Z(), 
                      cleanSameLocator = TRUE)
vario = Vario_computeFromDb(varioPar, grid)
p = plot.init() + 
  plot.varmod(vario = vario, model = mod, flagLegend=TRUE) +
  plot.decoration(title = title)
plot.end(p)

# variograms on all simulations
vl = list()
for (s in 1:ns) {
  err = grid$setLocator(name = paste(prefix, "V1", paste0("S", s), sep = "."), 
                        locatorType = ELoc_Z(), cleanSameLocator = TRUE)
  vl[[1 + length(vl)]] <- Vario_computeFromDb(varioPar, grid)
}
err = compare_variograms(var_list = vl, title = title)
hh = as.matrix((0:nlag)*dx)
gg = mod$sample(hh, mode=CovCalcMode_create(asVario=TRUE))
lines(hh, gg, col = "skyblue" )

2D simulation of a random function

ndim = 2
err = defineDefaultSpace(ESpaceType_RN(), ndim)

# specification of the model (theoretical range is provided in 'a')
mod  = Model_createFromParam(type = ECov_fromKey(type_cov), 
                             ranges = ranges[1:ndim], param = nu, angles=angles[1:ndim],
                             flagRange=FALSE)

# Two intersecting grids
grid1 = DbGrid_create(nx = c(256, 256))
grid2 = DbGrid_create(nx = c(100, 100), x0 = c(200, 200), dx = c(1,1))


# computing the simulation on two intersecting grids
seed <- 3132
err = simuSpectral(NULL, grid1, mod$getCov(), 1, nb, seed=seed , namconv = NamingConvention(prefix))
err = simuSpectral(NULL, grid2, mod$getCov(), 1, nb, seed=seed , namconv = NamingConvention(prefix))

p = plot.init(asp=1) +
  plot.raster(grid1, name = var_name, palette = "Spectral", limits = 3.5 * c(-1,1)) +
  plot.raster(grid2, name = var_name, palette = "Spectral", limits = 3.5 * c(-1,1)) + 
  plot.decoration(xlab = "Easting", ylab = "Northing", title = title)
plot.end(p)

# Control that the values on the intersection of the two grids are identical
err = migrate(dbin = grid1, dbout = grid2, name = var_name, 
              namconv = NamingConvention("grid1"))
p = plot.init() +
  plot.correlation(grid2, namex = var_name, namey = paste("grid1", var_name, sep = "."), asPoint = TRUE,
                   flagDiag = TRUE, flagSameAxes = TRUE) +
  plot.decoration(xlab = "simulation on grid1", ylab = "simulation on grid2",
                  title = title)
plot.end(p)

# computing the simulation on scattered points (extracted from grid1)
grid1["sel"] <- VectorHelper_simulateBernoulli(grid1$getNSample(), proba=0.005)
err = grid1$setLocator("sel", ELoc_SEL(), cleanSameLocator = TRUE)
db = Db()
db["x1"] <- grid1$getColumn(name = "x1", useSel = TRUE)
db["x2"] <- grid1$getColumn(name = "x2", useSel = TRUE)
db[paste("grid1", var_name, sep = ".")] <- grid1$getColumn(name = var_name, useSel = TRUE)
err = db$setLocators(c("x1", "x2"), ELoc_X(), cleanSameLocator = TRUE)
err = grid1$clearSelection()


# computing the simulation on scattered points
err = simuSpectral(NULL, db, mod$getCov(), 1, nb, seed=seed , namconv = NamingConvention(prefix))

p1 = plot.init(asp=1) +
  plot.raster(grid1, name = var_name, palette = "Spectral") +
  plot.symbol(db, color = "red") + 
  plot.decoration(xlab = "Easting", ylab = "Norhting", title = title)

# control that the values on the intersection of the two grids are identical
p2 = plot.init() +
  plot.correlation(db, namex = var_name, namey = paste("grid1", var_name, sep = "."), 
                   asPoint = TRUE, 
                   flagDiag = TRUE, flagSameAxes = TRUE) +
  plot.decoration(xlab = "simulation on db", ylab = "simulation on grid1", 
                  title = title)

ggarrange(p1, p2, nrow = 1, ncol = 2)

# statistics
tab = rbind(
  dbStatisticsMono(db = grid1, names = var_name, opers = opers)$toTL(),
  dbStatisticsMono(db = grid2, names = c(var_name), opers = opers, 
                   flagIso = TRUE)$toTL(),
  dbStatisticsMono(db = db, names = c(var_name), opers = opers, 
                   flagIso = TRUE)$toTL())
rownames(tab) = c("grid1", "grid2", "scattered points")
knitr::kable(tab,  digits = 3, caption = title)
MATERN model (nu = 0.75)
Number Minimum Maximum Mean St. Dev.
grid1 65536 -3.953 3.799 -0.066 0.998
grid2 10000 -2.844 3.455 -0.085 0.973
scattered points 357 -3.108 2.879 -0.031 0.971 
# variogram on the simulations
nlag = 25
coldir = c("green", "blue", "red")
err = grid1$setLocator(var_name, ELoc_Z(), cleanSameLocator = TRUE)
varioPar = VarioParam_createMultipleFromGrid(grid1, nlag=nlag)

# computing the experimental variogram  
vario = Vario_computeFromDb(varioPar, grid1)

p = plot.init() + 
  plot.varmod(vario = vario, model=mod, flagLegend=TRUE) + 
  plot.decoration(title = title)
plot.end(p)

3D simulation of a random function

ndim = 3
err = defineDefaultSpace(ESpaceType_RN(), ndim)

# specification of the model (theoretical range is provided in 'a')
mod  = Model_createFromParam(type = ECov_fromKey(type_cov), 
                             ranges = ranges[1:ndim], param = nu, angles=angles[1:ndim],
                             flagRange=FALSE)

# Two intersecting grids
nx <- c(100, 100, 20)
dx <- c(1., 1., 1.)
x0 <- c(10, 1000, 2)
grid = DbGrid_create(nx = nx, x0 = x0, dx = dx)

# computing the simulation on the grid
err = simuSpectral(NULL, grid, mod$getCov(), 1, nb, seed=seed , namconv = NamingConvention(prefix))
# statistics
knitr::kable(dbStatisticsMono(db = grid, names = var_name, opers = opers)$toTL(),
             digits = 3, caption = title)
MATERN model (nu = 0.75)
Number Minimum Maximum Mean St. Dev.
Simu.V1.S1 2e+05 -3.299 3.247 -0.074 0.952 
# histogram
p = plot.init() + 
  plot.hist(db = grid, name = var_name, bins = 25, col = "skyblue")
plot.end(p)

# map
p = plot.init(asp=1) + 
  plot.raster(grid, name = var_name, palette = "Spectral", limits = 3.5 * c(-1,1)) +
  plot.decoration(xlab = "Easting", ylab = "Northing", title = title)
plot.end(p)

# variogram on the simulations
nlag = 20
coldir = c("green", "red", "blue")
err = grid$setLocator(var_name, ELoc_Z(), cleanSameLocator = TRUE)
varioPar = VarioParam_createMultipleFromGrid(grid, nlag=nlag)

# computing the experimental variogram  
vario = Vario_computeFromDb(varioPar, grid)

p = plot.init() + 
  plot.varmod(vario = vario, model=mod, flagLegend=TRUE) + 
  plot.decoration(title = title)
plot.end(p)