Simulation Post-Processing for calculating Layer Proportions
Xavier Freulon
2023-07-24
Introduction
This R Markdown document tests the utilities developed to compute the layer proportions for each selected cells of a three dimensional grid. The \(N\) layers are defined by the elevation of their top limit.
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
library(gstlearn)## Loading required package: ggplot2## Loading required package: ggpubr## Loading required package: ggnewscale##
## Attaching package: 'gstlearn'## The following objects are masked from 'package:base':
##
## encodeString, messageset.seed(43243)
opers = EStatOption_fromKeys(c("NUM","MINI","MAXI","MEAN","STDV"))Parameters
# The 3D grid is defined on [0,1]^3
# 2D limits of the 3D grid
limits <- list(XP = 1.0*c(0, 1, 1, 0, 0),
YP = 1.0*c(0, 0, 1, 1, 0))
p_lim = PolyElem(limits$XP, limits$YP)
pol_lim = Polygons()
err = pol_lim$addPolyElem(p_lim)
flag.rotation = TRUE
# ----------------------------------------------
# ----------------------------------------------
if(flag.rotation) {
# 3D grid
nx = c(10, 13, 7)
dx = 1/nx
x0 = dx/2
# 2D grid
nx_2d = c(256, 256)
dx_2d = sqrt(2.0) / nx_2d
x0_2d = c(1/2, -1/2)
angles_2d = c(45, 0) # in degrees
} else { # no rotation
# 3D grid
nx = c(10, 10, 10)
dx = 1/nx
x0 = dx/2
# 2D grid
nx_2d = nx[1:2]
dx_2d = dx[1:2]
x0_2d = dx_2d/2
angles_2d = c(0, 0) # in degrees
}
grid = DbGrid_create(nx = nx, dx = dx, x0 = x0)
surfaces = DbGrid_create(nx = nx_2d, dx = dx_2d, x0 = x0_2d, angles = angles_2d)
# plots
p = plot.init(asp=1) +
plot.raster(grid, name = "rank", useSel = TRUE, flagLegend = FALSE) +
plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
plot.decoration(xlab = "Easting", ylab = "Northing", title = paste0("3D grid and limits"))
plot.end(p)
p = plot.init(asp=1) +
plot.raster(surfaces, name = "rank", useSel = TRUE, flagLegend = FALSE) +
plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
plot.decoration(xlab = "Easting", ylab = "Northing", title = paste0("2D grid and limits"))
plot.end(p)
# ----------------------
# selection of some grid cells
# ----------------------
sel = (grid["x1"] >= 0.20)&(grid["x1"] <= 0.70) &
(grid["x3"] >= 0.10)&(grid["x3"] <= 0.70)
sel = as.numeric(sel)
err = grid$setColumn(tab = sel, name = "sel")
err = grid$setLocator(name = "sel", ELoc_SEL(), cleanSameLocator = TRUE)
# ----------------------
# plot of a slice
# ----------------------
var_nm = "rank"
posx = 1; xlab = "Ox"
posy = 3; ylab = "Oz"
k = 3
corner = rep(k - 1, 3);
corner[posx] = 0;
corner[posy] = 0
slice = DbGrid_create(nx = c(nx[posx], nx[posy]), dx = c(dx[posx], dx[posy]), x0 = c(x0[posx], x0[posy]))
val = grid$getOneSlice(name = var_nm, posx = posx - 1, posy = posy - 1 , corner = c(0,0,k-1),
useSel = TRUE)
err = slice$setColumn(tab = val, name = paste0("K = ", k))
p = plot.init(asp=1) +
plot.raster(slice, name = paste0("K = ", k), useSel = TRUE, legendName = var_nm,
flagLegend = TRUE) +
plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
plot.decoration(xlab = xlab, ylab = ylab,
title = paste0("Slice K = ", k))
plot.end(p)
# -----------------------------------
# simulation of the surfaces
# -----------------------------------
N = 3
h_mean = c(0.5, 0.25, 0.25)
nbsim = c(3, 4, 3)
# N = 2
# h_mean = c(0.75, 0.25)
# nbsim = c(3, 4)
model = Model_createFromParam(ECov_GAUSSIAN(), range=0.5, sill=1.0)
for (i in 1:N) {
err = simtub(dbin = NULL, dbout = surfaces, model = model, nbsimu = nbsim[i],
namconv = NamingConvention(paste0("Y_", i)))
}
# conversion into uniform thickness
for (i in 1:N) {
val = matrix(surfaces$getColumns(names = paste0("Y_", i , ".S", 1:nbsim[i]),
useSel = FALSE),
nrow = surfaces$getNSample(useSel = FALSE), ncol = nbsim[i])
val = h_mean[i] * pnorm(val)
for (k in 1:nbsim[i]) {
err = surfaces$setColumn(tab = val[, k], name = paste0("Z_", i , ".", k),
useSel = FALSE)
}
}
# plot of a variable
var_nm = "Z_1.3"
p = plot.init(asp=1) +
plot.raster(surfaces, name = var_nm, useSel = TRUE, legendName = var_nm,
flagLegend = TRUE) +
plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
plot.decoration(xlab = "Easting", ylab = "Northing",
title = paste0("2D Surface: variable = ", var_nm))
plot.end(p)
Computing the proportions
The simulated variables are the thickness \(Z_i^{(k)}(s)\).
Sample selection
For a block selected by \((i, j, k)\) the samples \((x,y,z)\) of the surfaces are selected by:
$|x - (x_0[1]+ dx[1] (i - 1)| dx[1] /2 $
$|y - (x_0[2]+ dx[2] (j - 1)| dx[2] /2 $
# compute the rank from ijk of the values of an array defined by its dimensions nn
# This is a vector version of ijk2rank and rank2ijk for any number of dimensions
# rank = i + nn[1] (j - 1) + nn[1]*nn[2]*(k - 1) + ...
rank2ijkVec <- function(rank, nn) {
stopifnot((range(rank)[1] > 0)&(range(rank)[2] <= prod(nn)))
N = length(nn)
IJK = matrix(NaN, nrow = length(rank), ncol = N)
rr = rank - 1
for (i in 1:N) {
IJK[,i] = 1 + rr %% nn[i]
if(i < N) {rr = (rr - (IJK[,i] - 1)) %/% nn[i]}
}
if(length(rank) == 1) {IJK = as.numeric(IJK)}
IJK
}
# convert the ijk corresponding to the rank of an array defined by its dimensions nn
ijk2rankVec <- function(ijk, nn) {
if(!is.matrix(ijk)){ijk = matrix(ijk, nrow = 1, ncol = length(ijk))}
stopifnot(dim(ijk)[2] == length(nn))
as.numeric(1 + (ijk-1) %*% cumprod(c(1, nn[-length(nn)])))
}
# test on a 4d array
nn = c(5, 13, 7, 31)
rank = 1:prod(nn)
ijk = rank2ijkVec(rank = rank, nn = nn)
stopifnot(max(abs(rank - ijk2rankVec(rank2ijkVec(rank,nn),nn))) == 0)
stopifnot(max(abs(ijk - rank2ijkVec(ijk2rankVec(ijk,nn),nn))) == 0)
# test on a cell
cell = c(3, 4, 5)
rank = ijk2rankVec(ijk = cell, nn = grid$getNXs())
stopifnot(rank == ijk2rankVec(ijk = rank2ijkVec(rank = rank, nn = grid$getNXs()),
nn = grid$getNXs()))
stopifnot(max(abs(cell - rank2ijkVec(rank = ijk2rankVec(ijk = cell, nn = grid$getNXs()), nn = grid$getNXs()))) == 0)
# select the samples in dbin belonging to the grid cell and
# return the number of selected samples
# the selection works only with "not rotated" grid
select_samples <- function(dbin, rank, grid, nameSel = "inBlock") {
stopifnot(max(abs(grid$getAngles())) == 0) # no rotation
cell = rank2ijkVec(rank = rank, nn = grid$getNXs())
x0 = grid$getX0s()
dx = grid$getDXs()
xCell = x0 + dx * (cell - 1)
u1 = dbin$getColumn(name = "x1", useSel = FALSE)
u2 = dbin$getColumn(name = "x2", useSel = FALSE)
sel = as.numeric( (abs(u1 - xCell[1]) <= dx[1]/2)&(abs(u2 - xCell[2]) <= dx[2]/2))
err = dbin$deleteColumn(name = nameSel)
err = dbin$setColumn(tab = sel, name = nameSel, useSel=FALSE)
err = dbin$setLocator(name = nameSel, locatorType = ELoc_SEL(), cleanSameLocator = TRUE)
dbin$getNSample(useSel = TRUE)
}
err = surfaces$deleteColumn("inBlock")
err = surfaces$deleteColumn("sel_*")
print(paste0(">>> number of selected samples = ", select_samples(surfaces, rank = ijk2rankVec(ijk = cell, nn = grid$getNXs()), grid)))## [1] ">>> number of selected samples = 238"# plot of selected samples of the surface grid
err = surfaces$clearSelection()
err = surfaces$setLocator(name = "inBlock", locatorType = ELoc_SEL(),
cleanSameLocator = TRUE)
var_nm = "Z_1.3"
p = plot.init(asp=1) +
plot.raster(surfaces, name = var_nm, useSel = TRUE, legendName = var_nm,
flagLegend = TRUE) +
plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
plot.decoration(xlab = "Easting", ylab = "Northing",
title = paste0("2D surface: variable = ", var_nm))
plot.end(p)
u1 = surfaces$getColumn(name = "x1", useSel = TRUE)
u2 = surfaces$getColumn(name = "x2", useSel = TRUE)
#
X0 = grid$getX0s()
DX = grid$getDXs()
ll = list( x = c(
X0[1] + DX[1] * (cell[1] - 1 - 1/2),
X0[1] + DX[1] * (cell[1] - 1 + 1/2 ),
X0[1] + DX[1] * (cell[1] - 1 + 1/2),
X0[1] + DX[1] * (cell[1] - 1 - 1/2),
X0[1] + DX[1] * (cell[1] - 1 - 1/2)),
y = c(
X0[2] + DX[2] * (cell[2] - 1 - 1/2),
X0[2] + DX[2] * (cell[2] - 1 - 1/2),
X0[2] + DX[2] * (cell[2] - 1 + 1/2),
X0[2] + DX[2] * (cell[2] - 1 + 1/2),
X0[2] + DX[2] * (cell[2] - 1 - 1/2)))
plot(ll$x, ll$y, col = "orange", type = "l")
abline(v = X0[1] + DX[1] * (cell[1] - c(-1,1)/2), col = "gray", lty = 2)
abline(h = X0[2] + DX[2] * (cell[2] - 1 - c(-1,1)/2), col = "gray", lty = 2)
points(u1, u2, col = "red", pch = 19)
Computation of the multivariate simulation
flag.match = FALSE
ranks = grid$getColumn(name = "rank", useSel = TRUE)
# processing some cells
cell_selected = sort(sample(ranks, size = 5, replace = FALSE))
# --------------------------------------------------------------------------
# combining the univariate simulations into multivariate simulations
# --------------------------------------------------------------------------
if (flag.match) {
K = min(nbsim)
idx = matrix(NaN, nrow = K, ncol = N)
for (k in 1:K) {idx[k,] = k}
} else {
K = prod(nbsim)
idx = rank2ijkVec(rank = 1:K, nn = nbsim)
} # flag.match
P = N + 1
print(paste0(">>> computing ", P, " proportions over ", K, " simulations"))## [1] ">>> computing 4 proportions over 36 simulations"for (r in cell_selected) {
cell = rank2ijkVec(rank = r, nn = grid$getNXs())
ns = select_samples(surfaces, r, grid)
# mono-variate simulations of the N variables for the selected samples
H = list()
for (n in 1:N) {
H[[1+length(H)]] <- matrix(surfaces$getColumns(names = paste0("Z_", n, ".*"),
useSel = TRUE),
nrow = surfaces$getNSample(useSel = TRUE),
ncol = nbsim[n]
)
}
z_base = (cell[3]-1)/grid$getNXs()[3]
z_top = (cell[3])/grid$getNXs()[3]
h_max = z_top - z_base
# code = 0 -Inf < z <= H1
# code = 1 H1 < z <= H1 + H2
# code = 2 H1 + H2 < z <= H1 + H2 + H3
# code = 3 H1 + H2 + H3 < z <= + Inf
res = array(0.0, c(ns, K, P))
# bottom to top
for (i in 1:(P-1)) {
for (j in 1:i) {res[,,i] = res[,,i] + H[[j]][,idx[,j]]}
res[,,i] = pmin(pmax(res[,,i] - z_base,0), h_max)
}
res[,,P] = h_max
# bottom to top
for (i in P:2){
res[,,i] = res[,,i] - res[,,i-1]
}
res = res / h_max
# -----------------------------------------------------------------
# up-scaling (averaging the ns samples) to defined a P x K matrix
# -----------------------------------------------------------------
up_res = apply(X = res, MARGIN = c(2,3), FUN = mean)
# -----------------------------------------------------------------
# statistics on the K simulation up-scaling (averaging the ns samples)
# to define a P x K matrix
# -----------------------------------------------------------------
count <- function(u){length(u) - sum(is.na(u))}
opers = c("count", "min", "max", "median", "mean", "sd")
tab = matrix(NaN, nrow = P, ncol = length(opers))
for (j in seq_along(opers)) {
tab[, j] = apply(X = up_res, MARGIN = 2, FUN = opers[j])
}
colnames(tab) <- opers
rownames(tab) <- paste0("Code = ", 0:(P-1))
print(paste0("== Cell #", r, "/", grid$getNSample(useSel=FALSE),
" (regrouping ", ns, " samples)"))
titre = paste0("Statistics on ", ns, " samples",
" for cell ", cell[1], "/", cell[2], "/", cell[3],
" over the ", K, " simulations")
print(knitr::kable(tab, digits = 3, caption = titre))
} # loop over the cells of the grid## [1] "== Cell #213/910 (regrouping 250 samples)"
##
##
## Table: Statistics on 250 samples for cell 3/9/2 over the 36 simulations
##
## | | count| min| max| median| mean| sd|
## |:--------|-----:|---:|----:|------:|-----:|-----:|
## |Code = 0 | 36| 0| 1.00| 0| 0.333| 0.478|
## |Code = 1 | 36| 0| 1.00| 0| 0.247| 0.383|
## |Code = 2 | 36| 0| 0.98| 0| 0.140| 0.276|
## |Code = 3 | 36| 0| 1.00| 0| 0.280| 0.405|
## [1] "== Cell #294/910 (regrouping 247 samples)"
##
##
## Table: Statistics on 247 samples for cell 4/4/3 over the 36 simulations
##
## | | count| min| max| median| mean| sd|
## |:--------|-----:|---:|-----:|------:|-----:|-----:|
## |Code = 0 | 36| 0| 0.984| 0.000| 0.328| 0.470|
## |Code = 1 | 36| 0| 0.543| 0.016| 0.131| 0.199|
## |Code = 2 | 36| 0| 0.919| 0.025| 0.245| 0.299|
## |Code = 3 | 36| 0| 1.000| 0.081| 0.297| 0.394|
## [1] "== Cell #346/910 (regrouping 258 samples)"
##
##
## Table: Statistics on 258 samples for cell 6/9/3 over the 36 simulations
##
## | | count| min| max| median| mean| sd|
## |:--------|-----:|-----:|-----:|------:|-----:|-----:|
## |Code = 0 | 36| 0.000| 0.005| 0| 0.002| 0.002|
## |Code = 1 | 36| 0.000| 0.544| 0| 0.076| 0.173|
## |Code = 2 | 36| 0.000| 0.439| 0| 0.063| 0.126|
## |Code = 3 | 36| 0.112| 1.000| 1| 0.859| 0.253|
## [1] "== Cell #423/910 (regrouping 238 samples)"
##
##
## Table: Statistics on 238 samples for cell 3/4/4 over the 36 simulations
##
## | | count| min| max| median| mean| sd|
## |:--------|-----:|---:|----:|------:|-----:|-----:|
## |Code = 0 | 36| 0| 0.08| 0.000| 0.027| 0.038|
## |Code = 1 | 36| 0| 0.92| 0.000| 0.238| 0.365|
## |Code = 2 | 36| 0| 1.00| 0.032| 0.216| 0.303|
## |Code = 3 | 36| 0| 1.00| 0.648| 0.519| 0.460|
## [1] "== Cell #507/910 (regrouping 260 samples)"
##
##
## Table: Statistics on 260 samples for cell 7/12/4 over the 36 simulations
##
## | | count| min| max| median| mean| sd|
## |:--------|-----:|-----:|-----:|------:|-----:|-----:|
## |Code = 0 | 36| 0.000| 0.000| 0.000| 0.000| 0.000|
## |Code = 1 | 36| 0.000| 0.380| 0.020| 0.082| 0.119|
## |Code = 2 | 36| 0.263| 0.660| 0.422| 0.434| 0.112|
## |Code = 3 | 36| 0.101| 0.729| 0.541| 0.484| 0.212|Using convenient post-simulation processing function
err = surfaces$clearSelection()
err = simuPostPropByLayer(dbin = surfaces, dbout = grid, names = paste0("Z_", 1:length(nbsim), "*"),
flag_match=FALSE,
upscale = EPostUpscale_fromKey("MEAN"),
stats = EPostStat_fromKeys(c("MINI","MAXI","MED","MEAN","STD")),
verbose=TRUE, check_targets = cell_selected, check_level=0
)##
## Simulation Post-Processing
## --------------------------
## Multiplicity order for all variables
## - Variable 1 (Z_1*) = 3
## - Variable 2 (Z_2*) = 4
## - Variable 3 (Z_3*) = 3
## Number of Iterations: 36 (using the 'product' criterion)
## Number of Statistics: 5
## Number of Transform Variables: 4
## Number of Output Variables: 20
##
## == Cell #213/910 (regrouping 250 samples)
## Minimum Maximum Median Mean Std.
## Var 1 0.000 1.000 0.000 0.333 0.478
## Var 2 0.000 1.000 0.000 0.247 0.383
## Var 3 0.000 0.980 0.000 0.140 0.276
## Var 4 0.000 1.000 0.000 0.280 0.405
##
## == Cell #294/910 (regrouping 247 samples)
## Minimum Maximum Median Mean Std.
## Var 1 0.000 0.984 0.000 0.328 0.470
## Var 2 0.000 0.543 0.016 0.131 0.199
## Var 3 0.000 0.919 0.025 0.245 0.299
## Var 4 0.000 1.000 0.081 0.297 0.394
##
## == Cell #346/910 (regrouping 258 samples)
## Minimum Maximum Median Mean Std.
## Var 1 0.000 0.005 0.000 0.002 0.002
## Var 2 0.000 0.544 0.000 0.076 0.173
## Var 3 0.000 0.439 0.000 0.063 0.126
## Var 4 0.112 1.000 1.000 0.859 0.253
##
## == Cell #423/910 (regrouping 238 samples)
## Minimum Maximum Median Mean Std.
## Var 1 0.000 0.080 0.000 0.027 0.038
## Var 2 0.000 0.920 0.000 0.238 0.365
## Var 3 0.000 1.000 0.032 0.216 0.303
## Var 4 0.000 1.000 0.648 0.519 0.460
##
## == Cell #507/910 (regrouping 260 samples)
## Minimum Maximum Median Mean Std.
## Var 1 0.000 0.000 0.000 0.000 0.000
## Var 2 0.000 0.380 0.020 0.082 0.119
## Var 3 0.263 0.660 0.422 0.434 0.112
## Var 4 0.101 0.729 0.541 0.484 0.212Testing the standard deviation
U = runif(10)
# standard deviation on the samples
u = round(sd(U), digits = 6)
# formulae
local_mean <- function(u) {sum(u)/length(u)}
local_varP <- function(u) {sum(u^2)/length(u) - local_mean(u)^2}
local_var <- function(u) {m = local_mean(u); sum((u-m)^2)/(length(u) - 1)}
local_var2 <- function(u) {(sum(u^2) - length(u) * mean(u)^2)/(length(u) - 1)}
stopifnot(abs(sd(U) - sqrt(local_var(U))) < 1.e-15)
stopifnot(abs(sd(U) - sqrt(local_var2(U))) < 1.e-15)