GEONE - DEESSE - Simulations with probability constraints

Main points addressed

deesse simulation with probability (proportion) constraints:
- categorical or continuous variable
- local or global probability constraints

Note: if global probability constraints are used and if deesse is launched in parallel (with more than one thread), the reproducibility is not guaranteed.

Import what is required

[1]:

import numpy as np
import matplotlib.pyplot as plt
import time
import os

# import package 'geone'
import geone as gn

[2]:

# Show version of python and version of geone
import sys
print(sys.version_info)
print('geone version: ' + gn.__version__)

sys.version_info(major=3, minor=13, micro=7, releaselevel='final', serial=0)
geone version: 1.3.0

1. Categorical simulations - global probability constraints

Training image (TI)

Read the training image. Source of the image: D. Allard, D. D’or, and R. Froidevaux, An efficient maximum entropy approach for categorical variable prediction, EUROPEAN JOURNAL OF SOIL SCIENCE, 62(3):381-393, JUN 2011,doi:10.1111/j.1365-2389.2011.01362.x.

[3]:

# Read file
data_dir = 'data'
filename = os.path.join(data_dir, 'ti2.txt')
ti = gn.img.readImageTxt(filename)

# Values in the TI
ti.get_unique()

[3]:

array([0., 1., 2.])

[4]:

# Setting for categories / colors
categ_val = [0, 1, 2]
categ_col = ['lightblue', 'darkgreen', 'orange']

plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(ti, categ=True, categVal=categ_val, categCol=categ_col, title='TI')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_6_0.png

Simulation grid

Define the simulation grid (number of cells in each direction, cell unit, origin).

[5]:

nx, ny, nz = 300, 300, 1         # number of cells
sx, sy, sz = ti.sx, ti.sy, ti.sz # cell unit
ox, oy, oz = 0.0, 0.0, 0.0       # origin (corner of the "first" grid cell)

Define the classes of values

First, the classes of values for which the proportions will be specified have to be defined.

A class is defined as an interval or a union of intervals:

cl = [a, b], with a < b, define the class cl of any numerical value $v$ such that $a\leqslant v < b$ ( $a$ is included and $b$ excluded);
cl = [[a1, b1], [a2, b2]], with a1 < b1, a2 < b2 , define the class cl of any numerical value $v$ such that $a1\leqslant v < b1$ or $a2\leqslant v < b2$ ; more than two sub-intervals can be given.

In categorical case, a class for a category has to be defined as an interval of lower bound $a$ (included) and upper bound $b$ (excluded) that contains the value $v$ used for the category. Adapting the bounds allows to exclude other categories.

A class can be defined as a union of intervals to gather several categories whose the category values are not “adjacent”.

[6]:

nclass = 3
class1 = [-0.5, 0.5] # interval [-0.5, 0.5[ (for facies code 0)
class2 = [ 0.5, 1.5] # interval [ 0.5, 1.5[ (for facies code 1)
class3 = [ 1.5, 2.5] # interval [ 1.5, 2.5[ (for facies code 2)
# classx = [[-0.5, 0.5],[ 1.5, 2.5]] # for the union [-0.5, 0.5[ U [1.5, 2.5[, containing facies codes 0 and 2
list_of_classes = [class1, class2, class3]

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

To save time and to avoid noisy simulations, probability constraints can be deactivated when the last pattern node (the farest away from the central cell) is at a distance less than a given value (deactivationDistance). Note that the distance is computed according to the units defined for the search neighborhood ellipsoid.

[7]:

global_pdf = [0.2, 0.5, 0.3]
sp = gn.deesseinterface.SoftProbability(
    probabilityConstraintUsage=1,   # global probability constraints
    nclass=nclass,                  # number of classes of values
    classInterval=list_of_classes,  # list of classes
    globalPdf=global_pdf,           # global target PDF (list of length nclass)
    comparingPdfMethod=5,           # method for comparing PDF's (see doc:
                                    #    help(gn.deesseinterface.SoftProbability))
    deactivationDistance=4.0,       # deactivation distance (checking PDF is deactivated for narrow patterns)
    constantThreshold=1.e-3)        # acceptation threshold

Fill the input structure for deesse and launch deesse

[8]:

nreal = 1
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='categ',
    TI=ti,
    softProbability=sp,      # set probability constraints
    distanceType='categorical',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.25,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input, nthreads=8)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 8 thread(s)]
deesseRun: DeeSse run complete
deesseRun: warnings encountered (1 times in all):
#  1: WARNING 00111: reproducibility guaranteed when using multiple threads and global probability constraint (performance possibly limited), but results differ from the serial version
Elapsed time: 2.2 sec

Retrieve the results (and display)

[9]:

# Retrieve the realization
sim = deesse_output['sim']

# Display
plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(sim[0], categ=True, categVal=categ_val, categCol=categ_col,
                       title='Sim. using global target proportions')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_16_0.png

Compare facies proportions (TI, simulation, target)

[10]:

ti.get_prop()

[10]:

(array([0., 1., 2.]), array([[0.51492767, 0.23114805, 0.25392428]]))

[11]:

sim[0].get_prop()

[11]:

(array([0., 1., 2.]), array([[0.19292222, 0.50857778, 0.2985    ]]))

[12]:

global_pdf # target pdf

[12]:

[0.2, 0.5, 0.3]

2. Categorical simulations - local probability constraints

Target proportions can be specified locally. For each cell, target proportions (for each class) in a region around the cell is considered. Hence, proportion maps are required as well as a support radius: the region is defined as the ensemble of the cells in the search neighborhood ellipsoid and at a distance to the central (simulated) cell inferior or equal to the given support radius. Note that the distance are computed according to the units defined for the search neighborhood ellipsoid.

Build local target proportions maps

[13]:

# Set an image with simulation grid geometry defined above, and no variable
im = gn.img.Img(nx, ny, nz, sx, sy, sz, ox, oy, oz, nv=0)

# Get the x and y coordinates of the centers of grid cell (meshgrid)
xx = im.xx()[0]
yy = im.yy()[0]

# Equivalent:
## xg, yg: coordinates of the centers of grid cell
#xg = ox + 0.5*sx + sx*np.arange(nx)
#yg = oy + 0.5*sy + sy*np.arange(ny)
#xx, yy = np.meshgrid(xg, yg) # create meshgrid from the center of grid cells

# Define proportion maps for each class
c = 0.6
p1 = xx + yy
p1 = c * (p1 - np.min(p1))/ (np.max(p1) - np.min(p1))
p2 = c - p1
p0 = 1.0 - p1 - p2 # constant map (1-c = 0.4)

local_pdf = np.zeros((nclass, nz, ny, nx))
local_pdf[0,0,:,:] = p0
local_pdf[1,0,:,:] = p1
local_pdf[2,0,:,:] = p2

# Set the nclass=3 variables of local_pdf in image im
im.append_var(local_pdf, varname=['p0', 'p1', 'p2'])

# Display
plt.subplots(1, 3, figsize=(17,5), sharey=True) # 1 x 3 sub-plots

plt.subplot(1,3,1)
gn.imgplot.drawImage2D(im, iv=0, vmin=0, vmax=c, title='Porportion for facies 0')

plt.subplot(1,3,2)
gn.imgplot.drawImage2D(im, iv=1, vmin=0, vmax=c, title='Porportion for facies 1')

plt.subplot(1,3,3)
gn.imgplot.drawImage2D(im, iv=2, vmin=0, vmax=c, title='Porportion for facies 2')

plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_22_0.png

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

[14]:

sp = gn.deesseinterface.SoftProbability(
    probabilityConstraintUsage=2,   # local probability constraints
    nclass=nclass,                  # number of classes of values
    classInterval=list_of_classes,  # list of classes
    localPdf=local_pdf,             # local target PDF
    localPdfSupportRadius=12.,      # support radius
    comparingPdfMethod=5,           # method for comparing PDF's (see doc:
                                    #    help(gn.deesseinterface.SoftProbability))
    deactivationDistance=4.0,       # deactivation distance (checking PDF is deactivated for narrow patterns)
    constantThreshold=1.e-3)        # acceptation threshold

Fill the input structure for deesse and launch deesse

[15]:

nreal = 1
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='categ',
    TI=ti,
    softProbability=sp,      # set probability constraints
    distanceType='categorical',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.25,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input, nthreads=8)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 8 thread(s)]
deesseRun: DeeSse run complete
Elapsed time: 2.2 sec

[16]:

# Retrieve the realization
sim = deesse_output['sim']

# Display
plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(sim[0], categ=True, categVal=categ_val, categCol=categ_col,
                       title='Sim. using local target proportions')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_27_0.png

2b. Categorical simulations - local probability constraints based on rejection

An alternative to deal with local probability constraints consists in applying rejection. The given probability maps are used to compute rejection probabilities as follows. Note that these maps will not correspond to target proportions. The rejection is applied during the scan of the training image; two rejection modes are available:

rejectionMode = 0: rejection first (before checking pattern (and other constraint)) according to acceptation probabilities proportional to p[i]/q[i] (for class i), where
- q is the marginal pdf of the scanned training image
- p is the given local pdf at the simulated node
rejectionMode = 1: rejection last (after checking pattern (and other constraint)) according to acceptation probabilities proportional to p[i] (for class i), where
- p is the given local pdf at the simulated node

Remarks:

using rejectionMode = 0, the given probability maps would correspond to proportions over an ensemble of realizations if the pattern (and other constraint) was not accounted for, i.e. if random pixels were drawn from the training image and only its marginal proportions mattered
using rejectionMode = 1, the pdf proportional to $q_{MPS}[i] \cdot p[i]$ is sampled, where $q_{MPS}$ denotes the pdf resulting from multiple-point statistics (checking pattern (and other constraint))

Rejection first ( `rejectionMode = 0`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

[17]:

sp = gn.deesseinterface.SoftProbability(
    probabilityConstraintUsage=3,   # local probability constraints based on rejection
    nclass=nclass,                  # number of classes of values
    classInterval=list_of_classes,  # list of classes
    localPdf=local_pdf,             # local target PDF
    rejectionMode=0,                # rejection first
    deactivationDistance=4.0)       # deactivation distance (checking PDF is deactivated for narrow patterns)

Fill the input structure for deesse and launch deesse

[18]:

nreal = 1
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='categ',
    TI=ti,
    softProbability=sp,      # set probability constraints
    distanceType='categorical',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.25,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input, nthreads=8)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 8 thread(s)]
deesseRun: DeeSse run complete
Elapsed time: 2.8 sec

[19]:

# Retrieve the realization
sim = deesse_output['sim']

# Display
plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(sim[0], categ=True, categVal=categ_val, categCol=categ_col,
                       title='Sim. using local target proportions')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_33_0.png

Rejection last ( `rejectionMode = 1`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

[20]:

sp = gn.deesseinterface.SoftProbability(
    probabilityConstraintUsage=3,   # local probability constraints based on rejection
    nclass=nclass,                  # number of classes of values
    classInterval=list_of_classes,  # list of classes
    localPdf=local_pdf,             # local target PDF
    rejectionMode=1,                # rejection last
    deactivationDistance=4.0)       # deactivation distance (checking PDF is deactivated for narrow patterns)

Fill the input structure for deesse and launch deesse

[21]:

nreal = 1
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='categ',
    TI=ti,
    softProbability=sp,      # set probability constraints
    distanceType='categorical',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.25,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input, nthreads=8)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 8 thread(s)]
deesseRun: DeeSse run complete
Elapsed time: 1.5 sec

[22]:

# Retrieve the realization
sim = deesse_output['sim']

# Display
plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(sim[0], categ=True, categVal=categ_val, categCol=categ_col,
                       title='Sim. using local target proportions')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_38_0.png

3. Continuous simulation - global probability constraints

Training image (TI)

Source of the image: T. Zhang, P. Switzer, and A. Journel, Filter-based classification of training image patterns for spatial simulation, MATHEMATICAL GEOLOGY, 38(1):63-80, JAN 2006,doi:10.1007/s11004-005-9004-x.

[23]:

# Read file
data_dir = 'data'
filename = os.path.join(data_dir, 'tiContinuous.txt')
ti = gn.img.readImageTxt(filename)

# Color settings
cmap='gray'
vmin, vmax = ti.vmin(), ti.vmax()

# Plot
plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(ti, cmap=cmap, vmin=vmin, vmax=vmax)
plt.title(f'TI variable: {ti.varname[0]}')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_40_0.png

Simulation grid

Define the simulation grid (number of cells in each direction, cell unit, origin).

[24]:

nx, ny, nz = 200, 200, 1         # number of cells
sx, sy, sz = ti.sx, ti.sy, ti.sz # cell unit
ox, oy, oz = 0.0, 0.0, 0.0       # origin (corner of the "first" grid cell)

Define the classes of values

Set the number of classes, and for each class define the ensemble of values as a (union of) interval(s).

[25]:

vmin, vmax = 0., 256.
nclass = 10
breaks = np.linspace(vmin, vmax, nclass+1)
list_of_classes = [np.array([[breaks[i], breaks[i+1]]]) for i in range(nclass)]

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

[26]:

global_pdf = np.repeat(1./nclass, nclass) # global pdf (proportion for each class), uniform

sp = gn.deesseinterface.SoftProbability(
    probabilityConstraintUsage=1,   # global probability constraints
    nclass=nclass,                  # number of classes of values
    classInterval=list_of_classes,  # list of classes
    globalPdf=global_pdf,           # global target PDF (list of length nclass)
    comparingPdfMethod=5,           # method for comparing PDF's (see doc:
                                    #    help(gn.deesseinterface.SoftProbability))
    deactivationDistance=2.0,       # deactivation distance (checking PDF is deactivated for narrow patterns)
    constantThreshold=0)            # acceptation threshold

Fill the input structure for deesse and launch deesse

[27]:

nreal = 1
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='code',
    TI=ti,
    softProbability=sp,       # set probability constraints
    distanceType='continuous',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.5,
    npostProcessingPathMax=0, # disable post-processing (to avoid loosing target proportions)
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input, nthreads=8)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 8 thread(s)]
deesseRun: DeeSse run complete
deesseRun: warnings encountered (1 times in all):
#  1: WARNING 00111: reproducibility guaranteed when using multiple threads and global probability constraint (performance possibly limited), but results differ from the serial version
Elapsed time: 6.4 sec

[28]:

# Retrieve the realization
sim = deesse_output['sim']

# Display
plt.subplots(1,2, figsize=(15,5)) # 1 x 2 sub-plots

plt.subplot(1,2,1)
gn.imgplot.drawImage2D(sim[0], cmap=cmap, vmin=vmin, vmax=vmax,
                       title='Sim. using global target proportions')

plt.subplot(1,2,2)
plt.hist(ti.val.reshape(-1), bins=breaks, density=True, color='lightblue', alpha=0.5, label='TI')
plt.hist(sim[0].val.reshape(-1), bins=breaks, density=True, color='orange', alpha=0.5, label='SIM')
plt.legend()
plt.title('Proportions of values on the classes')

plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_49_0.png

4. Continuous simulations - local probability constraints

Define new classes of values

[29]:

nclass = 2
class1 = [0., 50.]   # interval [0., 50.[ (low values)
class2 = [50., 256.] # interval [50., 256.[ (high values)
list_of_classes = [class1, class2]

Build local target proportions maps

[30]:

# Set an image with simulation grid geometry defined above, and no variable
im = gn.img.Img(nx, ny, nz, sx, sy, sz, ox, oy, oz, nv=0)

# Get the x and y coordinates of the centers of grid cell (meshgrid)
xx = im.xx()[0]
yy = im.yy()[0]

# Equivalent:
## xg, yg: coordinates of the centers of grid cell
#xg = ox + 0.5*sx + sx*np.arange(nx)
#yg = oy + 0.5*sy + sy*np.arange(ny)
#xx, yy = np.meshgrid(xg, yg) # create meshgrid from the center of grid cells

# Define proportion maps for each class
p0 = (yy - np.min(yy))/ (np.max(yy) - np.min(yy))
p1 = 1.0 - p0

local_pdf = np.zeros((nclass, nz, ny, nx))
local_pdf[0,0,:,:] = p0
local_pdf[1,0,:,:] = p1

# Set the nclass=2 variables of local_pdf in image im
im.append_var(local_pdf, varname=['p0', 'p1'])

# Display
plt.subplots(1, 2, figsize=(14,7), sharey=True) # 1 x 2 sub-plots

plt.subplot(1,2,1)
gn.imgplot.drawImage2D(im, iv=0, vmin=0, vmax=1, title='Porportion for low values')

plt.subplot(1,2,2)
gn.imgplot.drawImage2D(im, iv=1, vmin=0, vmax=1, title='Porportion for high values')

plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_53_0.png

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

[31]:

sp = gn.deesseinterface.SoftProbability(
    probabilityConstraintUsage=2,   # local probability constraints
    nclass=nclass,                  # number of classes of values
    classInterval=list_of_classes,  # list of classes
    localPdf=local_pdf,             # local target PDF
    localPdfSupportRadius=12.,      # support radius
    comparingPdfMethod=5,           # method for comparing PDF's (see doc:
                                    #    help(gn.deesseinterface.SoftProbability))
    deactivationDistance=4.0,       # deactivation distance (checking PDF is deactivated for narrow patterns)
    constantThreshold=0.001)        # acceptation threshold

Fill the input structure for deesse and launch deesse

[32]:

nreal = 1
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='code',
    TI=ti,
    softProbability=sp,      # set probability constraints
    distanceType='continuous',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.5,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input, nthreads=8)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 8 thread(s)]
deesseRun: DeeSse run complete
Elapsed time: 6.3 sec

[33]:

# Retrieve the realization
sim = deesse_output['sim']

# Display
plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(sim[0], cmap=cmap, vmin=vmin, vmax=vmax,
                       title='Sim. using local target proportions')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_58_0.png

4b. Continuous simulations - local probability constraints based on rejection

Rejection first ( `rejectionMode = 0`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

[34]:

sp = gn.deesseinterface.SoftProbability(
    probabilityConstraintUsage=3,   # local probability constraints based on rejection
    nclass=nclass,                  # number of classes of values
    classInterval=list_of_classes,  # list of classes
    localPdf=local_pdf,             # local target PDF
    rejectionMode=0,                # rejection first
    deactivationDistance=4.0)       # deactivation distance (checking PDF is deactivated for narrow patterns)

Fill the input structure for deesse and launch deesse

[35]:

nreal = 1
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='code',
    TI=ti,
    softProbability=sp,      # set probability constraints
    distanceType='continuous',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.5,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input, nthreads=8)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 8 thread(s)]
deesseRun: DeeSse run complete
Elapsed time: 8.4 sec

[36]:

# Retrieve the realization
sim = deesse_output['sim']

# Display
plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(sim[0], cmap=cmap, vmin=vmin, vmax=vmax,
                       title='Sim. using local target proportions')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_64_0.png

Rejection last ( `rejectionMode = 1`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

[37]:

sp = gn.deesseinterface.SoftProbability(
    probabilityConstraintUsage=3,   # local probability constraints based on rejection
    nclass=nclass,                  # number of classes of values
    classInterval=list_of_classes,  # list of classes
    localPdf=local_pdf,             # local target PDF
    rejectionMode=1,                # rejection last
    deactivationDistance=4.0)       # deactivation distance (checking PDF is deactivated for narrow patterns)

Fill the input structure for deesse and launch deesse

[38]:

nreal = 1
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='code',
    TI=ti,
    softProbability=sp,      # set probability constraints
    distanceType='continuous',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.5,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input, nthreads=8)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 8 thread(s)]
deesseRun: DeeSse run complete
Elapsed time: 5.9 sec

[39]:

# Retrieve the realization
sim = deesse_output['sim']

# Display
plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(sim[0], cmap=cmap, vmin=vmin, vmax=vmax,
                       title='Sim. using local target proportions')
plt.show()

../_images/notebooks_ex_deesse_06_proba_constraint_69_0.png

GEONE - DEESSE - Simulations with probability constraints

Main points addressed

Import what is required

1. Categorical simulations - global probability constraints

Training image (TI)

Simulation grid

Define the classes of values

Define probability constraints (class geone.deesseinterface.SoftProbability)

Fill the input structure for deesse and launch deesse

Retrieve the results (and display)

Compare facies proportions (TI, simulation, target)

2. Categorical simulations - local probability constraints

Build local target proportions maps

Define probability constraints (class geone.deesseinterface.SoftProbability)

Fill the input structure for deesse and launch deesse

2b. Categorical simulations - local probability constraints based on rejection

Rejection first ( rejectionMode = 0)

Define probability constraints (class geone.deesseinterface.SoftProbability)

Fill the input structure for deesse and launch deesse

Rejection last ( rejectionMode = 1)

Define probability constraints (class geone.deesseinterface.SoftProbability)

Fill the input structure for deesse and launch deesse

3. Continuous simulation - global probability constraints

Training image (TI)

Simulation grid

Define the classes of values

Define probability constraints (class geone.deesseinterface.SoftProbability)

Fill the input structure for deesse and launch deesse

4. Continuous simulations - local probability constraints

Define new classes of values

Build local target proportions maps

Define probability constraints (class geone.deesseinterface.SoftProbability)

Fill the input structure for deesse and launch deesse

4b. Continuous simulations - local probability constraints based on rejection

Rejection first ( rejectionMode = 0)

Define probability constraints (class geone.deesseinterface.SoftProbability)

Fill the input structure for deesse and launch deesse

Rejection last ( rejectionMode = 1)

Define probability constraints (class geone.deesseinterface.SoftProbability)

Fill the input structure for deesse and launch deesse

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

Rejection first ( `rejectionMode = 0`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

Rejection last ( `rejectionMode = 1`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

Rejection first ( `rejectionMode = 0`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)

Rejection last ( `rejectionMode = 1`)

Define probability constraints (class `geone.deesseinterface.SoftProbability`)