\frametitle{Preamble} 
rm(list=ls())
library(gstlearn)
library(ggplot2)
library(ggpubr)

Defining some essential parameters:

ndim = 2
defineDefaultSpace(ESpaceType_RN(), ndim);
## NULL
nbsimu = 20
nbcc   = 4
cmap   = c('red', 'blue', 'yellow')

Downloading the data base (from the distribution Exdemo_PGS.db) and creating the output Grid, the Model (Cubic) and the Neighborhood (Unique):

dlfile = "https://soft.minesparis.psl.eu/gstlearn/data/PluriGaussian/Data.NF"
fileNF = "Scotland_Temperatures.NF"
download.file(dlfile, fileNF)
data = Db_createFromNF(fileNF)

grid = DbGrid_create(nx=c(110,110))

model = Model_createFromParam(type=ECov_CUBIC(), ranges=c(50,30))

neigh = NeighUnique()
\frametitle{Defining Internal Display Function}

Defining the internal function plot.dat used to visualize the data set with convenient color representation:

\frametitle{Display Data} 
# TODO: Utiliser CMAP
plot.setDefaultGeographic(dims=c(6,6))
p = ggDefaultGeographic()
p = p + plot(data, name_color="facies", name_size="connect")
#         edgecolors='black', sizmin=40, cmap=cmap)
p = p + plot.decoration(title="Conditioning Information")
ggPrint(p)

\frametitle{Proportions and Lithotype Rule} 
props = dbStatisticsFacies(data)
rule = Rule_createFromNames(c("S","S","F1","F2","F3"))

# TODO: utiliser cmap
p = ggplot()
p = p + plot(rule, proportions = props, maxG=3)
p = p + plot.decoration(title="Lithotype Rule")
ggPrint(p)

\frametitle{Model of Underlying GRF}

Calculate the Experimental Variogram of the Underlying Gaussian Random Function and fit the Model (used in PGS).

varioparam = VarioParam_createOmniDirection(dpas=5, npas=20)
ruleprop = RuleProp_createFromRule(rule, props)
vario = variogram_pgs(data, varioparam, ruleprop)

model_gaus = Model()
err = model_gaus$fit(vario, types=ECov_fromKeys(c("CUBIC")), 
                     constraints=Constraints(1.))

ggplot() + plot.varmod(vario, model_gaus, asCov=TRUE)

\frametitle{PluriGaussian Simulation} 
err = simpgs(dbin=data, dbout=grid, ruleprop=ruleprop, model1=model_gaus, 
             neigh=neigh, nbsimu=nbsimu,
             namconv = NamingConvention("SimuPGS"))
p1 = ggDefaultGeographic()
p1 = p1 + plot(grid, "SimuPGS.1")
p1 = p1 + plot(data, name_color="facies", name_size="connect")
p2 = ggDefaultGeographic()
p2 = p2 + plot(grid, "SimuPGS.2")
p2 = p2 + plot(data, name_color="facies", name_size="connect")
p3 = ggDefaultGeographic()
p3 = p3 + plot(grid, "SimuPGS.3")
p3 = p3 + plot(data, name_color="facies", name_size="connect")
p4 = ggDefaultGeographic()
p4 = p4 + plot(grid, "SimuPGS.4")
p4 = p4 + plot(data, name_color="facies", name_size="connect")
ggarrange(p1,p2,p3,p4, nrow=2, ncol=2)

\frametitle{Acceptation Function}

Acceptation internal function: Select a Target Facies and build its Connected Components. For each simulation outcome, check the ranks of the connected component(s) at constraining wells and accept the simulation if all ranks are similar.

accept <- function(data, grid, name, verbose=FALSE, transBinary=TRUE, faccc=2)
{
  # Get the indices of samples which should be connected (starting from 0)
  rankData = which(data["connect"] == 1) - 1
  rankGrid = grid$locateDataInGrid(data, rankData)
    if (verbose)
    {
      cat("Number of conditioning data (connected) =",length(rankData),"\n")
    cat("Their ranks in the input Data Base =",rankData,"\n")
    cat("Their ranks in the output Data Base =",rankGrid,"\n")
    }
  
  # Perform the labelling into connected components
  err = grid$setLocator(name, ELoc_Z(), cleanSameLocator=TRUE)
  err = dbMorpho(grid, EMorpho_CC(), vmin=faccc-0.5, vmax=faccc+0.5)
  cc_list = grid[rankGrid,"Morpho*"]
  if (verbose)
    cat("List of their connected components indices =",cc_list,"\n")
  
  # Check that the data points belong to the same connected component
  number = length(unique(cc_list))
  retval = (number == 1)
  if (verbose)
    cat("Acceptation score =",retval,"\n")
        
  # Convert the valid Simulation outcome into a binary image
  if (retval && transBinary)
    grid[name] = (grid["Morpho*"] == cc_list[1])
    
  grid$deleteColumn("Morpho*")
  retval
}
\frametitle{Experiment the Acceptation Function on one Simulation outcome} 
isValid = accept(data, grid, "SimuPGS.1", TRUE)
## Number of conditioning data (connected) = 4 
## Their ranks in the input Data Base = 4 12 15 16 
## Their ranks in the output Data Base = 6000 1270 6648 9712 
## List of their connected components indices = 1 NA NA 1 
## Acceptation score = FALSE
cat("Connectivity for Simulation #1 :",isValid,"\n")
## Connectivity for Simulation #1 : FALSE
\frametitle{Probability Map}

This operation provides the probability that each pixel belongs to the Target Facies, calculated over all simulations that fulfill the Connectivity Constraint.

nb.valid = 0
for (i in 1:nbsimu)
{
  name = paste("SimuPGS.",i, sep="")
  cat("Simulation",name)
  isValid = accept(data, grid, name)
  if (isValid)
  {
    cat(" is valid\n")
    nb.valid = nb.valid + 1
  }
  else
  {
    cat(" is invalid\n")
    grid$deleteColumn(name)
  }
}
## Simulation SimuPGS.1 is invalid
## Simulation SimuPGS.2 is valid
## Simulation SimuPGS.3 is invalid
## Simulation SimuPGS.4 is valid
## Simulation SimuPGS.5 is invalid
## Simulation SimuPGS.6 is invalid
## Simulation SimuPGS.7 is invalid
## Simulation SimuPGS.8 is invalid
## Simulation SimuPGS.9 is invalid
## Simulation SimuPGS.10 is valid
## Simulation SimuPGS.11 is invalid
## Simulation SimuPGS.12 is invalid
## Simulation SimuPGS.13 is invalid
## Simulation SimuPGS.14 is valid
## Simulation SimuPGS.15 is invalid
## Simulation SimuPGS.16 is invalid
## Simulation SimuPGS.17 is invalid
## Simulation SimuPGS.18 is invalid
## Simulation SimuPGS.19 is invalid
## Simulation SimuPGS.20 is invalid
cat("Number of valid simulations =", nb.valid, "out of", nbsimu, "\n")
## Number of valid simulations = 4 out of 20
\frametitle{Probability Map}

Derive the Probability Map

err = grid$statistics(c("SimuPGS*"),
                      EStatOption_fromKeys(c("MEAN")), flagStoreInDb=TRUE)

Display

p = ggDefaultGeographic()
p = p + plot(grid,"Stats.MEAN")
p = p + plot(data,name_color="facies", name_size="connect")
p = p + plot.decoration(title="Probability of Connecting Wells")
ggPrint(p)