simtub.cpp File Reference

#include "geoslib_f.h"
#include "geoslib_old_f.h"
#include "geoslib_f_private.h"
#include "Enum/EJustify.hpp"
#include "Enum/ECov.hpp"
#include "Enum/EProcessOper.hpp"
#include "Enum/ERule.hpp"
#include "Calculators/CalcMigrate.hpp"
#include "Gibbs/GibbsUMultiMono.hpp"
#include "Gibbs/GibbsUPropMono.hpp"
#include "Gibbs/GibbsMMulti.hpp"
#include "Gibbs/GibbsFactory.hpp"
#include "Morpho/Morpho.hpp"
#include "Covariances/CovAniso.hpp"
#include "Basic/MathFunc.hpp"
#include "Basic/String.hpp"
#include "Basic/NamingConvention.hpp"
#include "Basic/Utilities.hpp"
#include "Basic/Law.hpp"
#include "Basic/File.hpp"
#include "Basic/OptDbg.hpp"
#include "Basic/OptCustom.hpp"
#include "LithoRule/PropDef.hpp"
#include "LithoRule/Rule.hpp"
#include "LithoRule/RuleShift.hpp"
#include "LithoRule/RuleShadow.hpp"
#include "LithoRule/RuleProp.hpp"
#include "Db/Db.hpp"
#include "Model/Model.hpp"
#include "Neigh/NeighUnique.hpp"
#include "Neigh/NeighMoving.hpp"
#include "Simulation/CalcSimuTurningBands.hpp"
#include "Simulation/SimuBoolean.hpp"
#include "Simulation/SimuBooleanParam.hpp"
#include "Simulation/SimuSpherical.hpp"
#include "Simulation/SimuSphericalParam.hpp"
#include "Simulation/SimuFFTParam.hpp"
#include "Simulation/SimuRefineParam.hpp"
#include "Simulation/SimuRefine.hpp"
#include "Simulation/CalcSimuEden.hpp"
#include "Simulation/CalcSimuFFT.hpp"
#include "Space/SpaceRN.hpp"
#include <math.h>
#include <string.h>

Functions
void	simu_func_categorical_transf (Db *db, int verbose, int isimu, int nbsimu)
void	simu_func_continuous_update (Db *db, int verbose, int isimu, int nbsimu)
void	simu_func_categorical_update (Db *db, int verbose, int isimu, int nbsimu)
void	simu_func_continuous_scale (Db *db, int verbose, int nbsimu)
void	simu_func_categorical_scale (Db *db, int verbose, int nbsimu)
void	check_mandatory_attribute (const char *method, Db *db, const ELoc &locatorType)
int	get_rank_from_propdef (PropDef *propdef, int ipgs, int igrf)
int	simpgs (Db *dbin, Db *dbout, RuleProp *ruleprop, Model *model1, Model *model2, ANeigh *neigh, int nbsimu, int seed, int flag_gaus, int flag_prop, int flag_check, int flag_show, int nbtuba, int gibbs_nburn, int gibbs_niter, double percent, const NamingConvention &namconv)
int	simbipgs (Db *dbin, Db *dbout, RuleProp *ruleprop, Model *model11, Model *model12, Model *model21, Model *model22, ANeigh *neigh, int nbsimu, int seed, int flag_gaus, int flag_prop, int flag_check, int flag_show, int nbtuba, int gibbs_nburn, int gibbs_niter, double percent, const NamingConvention &namconv)
int	db_simulations_to_ce (Db *db, const ELoc &locatorType, int nbsimu, int nvar, int *iptr_ce_arg, int *iptr_cstd_arg)
int	gibbs_sampler (Db *dbin, Model *model, int nbsimu, int seed, int gibbs_nburn, int gibbs_niter, bool flag_moving, bool flag_norm, bool flag_multi_mono, bool flag_propagation, bool flag_sym_neigh, int gibbs_optstats, double percent, bool flag_ce, bool flag_cstd, bool verbose, const NamingConvention &namconv)
int	simtub_constraints (Db *dbin, Db *dbout, Model *model, ANeigh *neigh, int seed, int nbtuba, int nbsimu_min, int nbsimu_quant, int niter_max, VectorInt &cols, int(*func_valid)(int flag_grid, int nDim, int nech, int *nx, double *dx, double *x0, double nonval, double percent, VectorDouble &tab))
int	simmaxstable (Db *dbout, Model *model, double ratio, int seed, int nbtuba, int flag_simu, int flag_rank, int verbose)
int	simRI (Db *dbout, Model *model, int ncut, double *zcut, double *wcut, int seed, int nbtuba, int verbose)
int	simpgs_spde (Db *dbin, Db *dbout, RuleProp *ruleprop, Model *model1, Model *model2, const String &triswitch, const VectorDouble &gext, int flag_gaus, int flag_prop, int flag_check, int flag_show, int nfacies, int seed, int nbsimu, int gibbs_nburn, int gibbs_niter, int ngibbs_int, int verbose, double percent)
int	simcond (Db *dbin, Db *dbout, Model *model, int seed, int nbsimu, int nbtuba, int gibbs_nburn, int gibbs_niter, int flag_check, int flag_ce, int flag_cstd, int verbose)
int	simsph (DbGrid *db, Model *model, const SimuSphericalParam &sphepar, int seed, bool verbose, const NamingConvention &namconv)
VectorDouble	simsph_mesh (MeshSpherical *mesh, Model *model, const SimuSphericalParam &sphepar, int seed, int verbose)
DbGrid *	simfine (DbGrid *dbin, Model *model, const SimuRefineParam &param, int seed)
MatrixRectangular	fluid_extract (DbGrid *dbgrid, const String &name_facies, const String &name_fluid, const String &name_poro, const String &name_date, int nfacies, int nfluids, int facies0, int fluid0, int ntime, double time0, double dtime, bool verbose)

Function Documentation

check_mandatory_attribute()

void check_mandatory_attribute	(	const char *	method,
		Db *	db,
		const ELoc &	locatorType
	)

Check for the presence of mandatory attributes

Parameters

[in]	method	Name of the method
[in]	db	Db structure
[in]	locatorType	Mandatory attribute type

db_simulations_to_ce()

int db_simulations_to_ce	(	Db *	db,
		const ELoc &	locatorType,
		int	nbsimu,
		int	nvar,
		int *	iptr_ce_arg,
		int *	iptr_cstd_arg
	)

Convert series of simulations to conditional expectation and variance

Returns

Error return code

Parameters

[in]	db	Db structure
[in]	locatorType	Type of pointer containing the simulations
[in]	nbsimu	Number of simulations
[in]	nvar	Number of variables
[out]	iptr_ce_arg	Pointer to the Conditional Expectation attributes
[out]	iptr_cstd_arg	Pointer to the Conditional St. Dev. attributes

fluid_extract()

MatrixRectangular fluid_extract	(	DbGrid *	dbgrid,
		const String &	name_facies,
		const String &	name_fluid,
		const String &	name_poro,
		const String &	name_date,
		int	nfacies,
		int	nfluids,
		int	facies0,
		int	fluid0,
		int	ntime,
		double	time0,
		double	dtime,
		bool	verbose
	)

Extract time charts from the fluid propagation block

Returns

The returned matrix

Parameters

[in]	dbgrid	Db grid structure
[in]	name_facies	Name of variable containing Facies
[in]	name_fluid	Name of variable containing Fluid
[in]	name_poro	Name of variable containing Porosity (optional)
[in]	name_date	Name of variable containing Date
[in]	nfacies	number of facies (facies 0 excluded)
[in]	nfluids	number of fluids
[in]	facies0	Value of the target facies
[in]	fluid0	Value of the target fluid
[in]	ntime	Number of Time intervals
[in]	time0	Starting time
[in]	dtime	Time interval
[in]	verbose	1 for a verbose option

get_rank_from_propdef()

int get_rank_from_propdef	(	PropDef *	propdef,
		int	ipgs,
		int	igrf
	)

Give the rank of a "variable" for a given GRF and PGS

Returns

Returned rank

Parameters

[in]	propdef	PropDef structure
[in]	ipgs	Rank of the GS
[in]	igrf	Rank of the Gaussian

gibbs_sampler()

int gibbs_sampler	(	Db *	dbin,
		Model *	model,
		int	nbsimu,
		int	seed,
		int	gibbs_nburn,
		int	gibbs_niter,
		bool	flag_moving,
		bool	flag_norm,
		bool	flag_multi_mono,
		bool	flag_propagation,
		bool	flag_sym_neigh,
		int	gibbs_optstats,
		double	percent,
		bool	flag_ce,
		bool	flag_cstd,
		bool	verbose,
		const NamingConvention &	namconv
	)

Perform the Gibbs sampler

Returns

Error return code

Parameters

[in]	dbin	Db structure
[in]	model	Model structure
[in]	nbsimu	Number of simulations
[in]	seed	Seed for random number generator
[in]	gibbs_nburn	Initial number of iterations for bootstrapping
[in]	gibbs_niter	Maximum number of iterations
[in]	flag_moving	1 for Moving
[in]	flag_norm	1 if the Model must be normalized
[in]	flag_multi_mono	1 for the Multi_mono algorithm
[in]	flag_propagation	1 for the propagation algorithm
[in]	flag_sym_neigh	Deprecated argument
[in]	gibbs_optstats	0: No stats - 1: Print - 2: Save Neutral file
[in]	percent	Amount of nugget effect added to too continuous model (expressed in percentage of total variance)
[in]	flag_ce	1 if the conditional expectation should be returned instead of simulations
[in]	flag_cstd	1 if the conditional standard deviation should be returned instead of simulations
[in]	verbose	Verbose flag
[in]	namconv	Naming convention

simbipgs()

int simbipgs	(	Db *	dbin,
		Db *	dbout,
		RuleProp *	ruleprop,
		Model *	model11,
		Model *	model12,
		Model *	model21,
		Model *	model22,
		ANeigh *	neigh,
		int	nbsimu,
		int	seed,
		int	flag_gaus,
		int	flag_prop,
		int	flag_check,
		int	flag_show,
		int	nbtuba,
		int	gibbs_nburn,
		int	gibbs_niter,
		double	percent,
		const NamingConvention &	namconv
	)

Perform the conditional or non-conditional Bi Pluri-gaussian simulations

Returns

Error return code

Parameters

[in]	dbin	Input Db structure (optional)
[in]	dbout	Output Db structure
[in]	ruleprop	Ruleprop definition
[in]	model11	First Model structure for First Lithotype Rule
[in]	model12	Second Model structure for First Lithotype Rule
[in]	model21	First Model structure for Second Lithotype Rule
[in]	model22	Second Model structure for Second Lithotype Rule
[in]	neigh	ANeigh structure
[in]	nbsimu	Number of simulations
[in]	seed	Seed for random number generator
[in]	flag_gaus	1 gaussian results; otherwise facies
[in]	flag_prop	1 for facies proportion
[in]	flag_check	1 if the facies at data must be checked against the closest simulated grid node
[in]	flag_show	1 if the grid node which coincides with the data should be represented with the data facies (only if flag_cond && !flag_gaus)
[in]	nbtuba	Number of turning bands
[in]	gibbs_nburn	Number of bootstrap iterations
[in]	gibbs_niter	Maximum number of iterations
[in]	percent	Amount of nugget effect added to too continuous model (expressed in percentage of the total variance)
[in]	namconv	Naming convention

Remarks

When conditional, the two first variables in the input Db

should correspond to the two facies indices (starting from 1)

The argument 'dbin' is optional: it must be defined only for

conditional simulations

The proportions (nfac1 * nfac2) must be ordered as follows:

f1af2a, f1bf2a, f1cf2a, ..., f1bf2a, f1bf2b, ..., f1nf2m

simcond()

int simcond	(	Db *	dbin,
		Db *	dbout,
		Model *	model,
		int	seed,
		int	nbsimu,
		int	nbtuba,
		int	gibbs_nburn,
		int	gibbs_niter,
		int	flag_check,
		int	flag_ce,
		int	flag_cstd,
		int	verbose
	)

Perform the conditional simulations under inequality constraints

Returns

Error return code

Parameters

[in]	dbin	Input Db structure (optional)
[in]	dbout	Output Db structure
[in]	model	Model structure
[in]	seed	Seed for random number generator
[in]	nbsimu	Number of simulations
[in]	nbtuba	Number of turning bands
[in]	gibbs_nburn	Initial number of iterations for bootstrapping
[in]	gibbs_niter	Maximum number of iterations
[in]	flag_check	1 to check the proximity in Gaussian scale
[in]	flag_ce	1 if the conditional expectation should be returned instead of simulations
[in]	flag_cstd	1 if the conditional standard deviation should be returned instead of simulations
[in]	verbose	Verbose flag

Remarks

The Neighborhood does not have to be defined as this method

only functions using a Unique Neighborhood. For consistency

it is generated internally.

simfine()

DbGrid* simfine	(	DbGrid *	dbin,
		Model *	model,
		const SimuRefineParam &	param,
		int	seed
	)

Refine the simulation

Returns

Newly refined Grid.

Parameters

[in]	dbin	Input grid Db structure
[in]	model	Model structure
[in]	param	SimuRefineParam structure
[in]	seed	Seed for the random number generator

Remarks

For each dimension of the space, if N stands for the number of

nodes in the input grid, the number of nodes of the output grid

will be (N-1) * 2^p + 1 where p is the param.getNmult()

simmaxstable()

int simmaxstable	(	Db *	dbout,
		Model *	model,
		double	ratio,
		int	seed,
		int	nbtuba,
		int	flag_simu,
		int	flag_rank,
		int	verbose
	)

Perform the non-conditional simulation of the Max-Stable Model

Returns

Error return code

Parameters

[in]	dbout	Output Db structure
[in]	model	Model structure
[in]	ratio	Ratio modifying the range at each iteration
[in]	seed	Seed for random number generator
[in]	nbtuba	Number of turning bands
[in]	flag_simu	1 if the simulation must be stored
[in]	flag_rank	1 if the iteration rank must be stored
[in]	verbose	Verbose flag

Remarks

This function uses a threshold that can be defined using

keypair mechanism with keyword "MaxStableThresh".

simpgs()

int simpgs	(	Db *	dbin,
		Db *	dbout,
		RuleProp *	ruleprop,
		Model *	model1,
		Model *	model2,
		ANeigh *	neigh,
		int	nbsimu,
		int	seed,
		int	flag_gaus,
		int	flag_prop,
		int	flag_check,
		int	flag_show,
		int	nbtuba,
		int	gibbs_nburn,
		int	gibbs_niter,
		double	percent,
		const NamingConvention &	namconv
	)

Perform the conditional or non-conditional Pluri-gaussian simulations

Returns

Error return code

Parameters

[in]	dbin	Input Db structure (optional)
[in]	dbout	Output Db structure
[in]	ruleprop	RuleProp structure
[in]	model1	First Model structure
[in]	model2	Second Model structure (optional)
[in]	neigh	ANeigh structure
[in]	nbsimu	Number of simulations
[in]	seed	Seed for random number generator
[in]	flag_gaus	1 if results must be gaussian; otherwise facies
[in]	flag_prop	1 for facies proportion
[in]	flag_check	1 if the facies at data must be checked against the closest simulated grid node
[in]	flag_show	1 if the grid node which coincides with the data should be represented with the data facies (only if flag_cond && !flag_gaus)
[in]	nbtuba	Number of turning bands
[in]	gibbs_nburn	Number of bootstrap iterations
[in]	gibbs_niter	Maximum number of iterations
[in]	percent	Amount of nugget effect added to too much continous model (expressed in percentage of the total variance)
[in]	namconv	Naming convention

Remarks

When conditional, the unique variable in the input Db structure

should correspond to the facies index (starting from 1)

The argument 'dbin' is optional: it must be defined only for

conditional simulations

simpgs_spde()

int simpgs_spde	(	Db *	dbin,
		Db *	dbout,
		RuleProp *	ruleprop,
		Model *	model1,
		Model *	model2,
		const String &	triswitch,
		const VectorDouble &	gext,
		int	flag_gaus,
		int	flag_prop,
		int	flag_check,
		int	flag_show,
		int	nfacies,
		int	seed,
		int	nbsimu,
		int	gibbs_nburn,
		int	gibbs_niter,
		int	ngibbs_int,
		int	verbose,
		double	percent
	)

Perform the conditional Pluri-gaussian simulations using spde

Returns

Error return code

Parameters

[in]	dbin	Input Db structure (optional)
[in]	dbout	Output Db structure
[in]	ruleprop	RuleProp definition
[in]	model1	First Model structure
[in]	model2	Second Model structure (optional)
[in]	triswitch	Meshing option
[in]	gext	Array of domain dilation
[in]	flag_gaus	1 if results must be gaussian; otherwise facies
[in]	flag_prop	1 for facies proportion
[in]	flag_check	1 if the facies at data must be checked against the closest simulated grid node
[in]	flag_show	1 if the grid node which coincides with the data should be represented with the data facies
[in]	nfacies	Number of facies
[in]	seed	Seed for random number generator
[in]	nbsimu	Number of simulations
[in]	gibbs_nburn	Number of iterations (Burning step)
[in]	gibbs_niter	Maximum number of iterations
[in]	ngibbs_int	Number of iterations internal to Gibbs (SPDE)
[in]	verbose	Verbose flag
[in]	percent	Amount of nugget effect added to too continous model (expressed in percentage of the total variance)

Remarks

When conditional, the unique variable in the input Db structure

should correspond to the facies index (starting from 1)

simRI()

int simRI	(	Db *	dbout,
		Model *	model,
		int	ncut,
		double *	zcut,
		double *	wcut,
		int	seed,
		int	nbtuba,
		int	verbose
	)

Perform the non-conditional simulation of the Orthogonal Residual Model

Returns

Error return code

Parameters

[in]	dbout	Output Db structure
[in]	model	Model structure
[in]	ncut	Number of cutoffs
[in]	zcut	Array of cutoffs
[in]	wcut	Array of weights
[in]	seed	Seed for random number generator
[in]	nbtuba	Number of turning bands
[in]	verbose	Verbose flag

simsph()

int simsph	(	DbGrid *	db,
		Model *	model,
		const SimuSphericalParam &	sphepar,
		int	seed,
		bool	verbose,
		const NamingConvention &	namconv
	)

Simulates the random function on the sphere

Parameters

[in]	db	Data base containing the coordinates of target points These coordinates must be expressed in long/lat
[in]	model	Model (defined in Euclidean space) to be used
[in]	sphepar	SimuSphericalParam structure
[in]	seed	Seed for random number generation
[in]	verbose	Verbose flag
[in]	namconv	Naming convention

simsph_mesh()

VectorDouble simsph_mesh	(	MeshSpherical *	mesh,
		Model *	model,
		const SimuSphericalParam &	sphepar,
		int	seed,
		int	verbose
	)

Simulates the random function on the sphere

Returns

The Vector simulated values

Parameters

[in]	mesh	MeshSpherical object
[in]	model	Model (defined in Euclidean space) to be used
[in]	sphepar	SimuSphericalParam structure
[in]	seed	Seed for random number generation
[in]	verbose	Verbose flag

simtub_constraints()

int simtub_constraints	(	Db *	dbin,
		Db *	dbout,
		Model *	model,
		ANeigh *	neigh,
		int	seed,
		int	nbtuba,
		int	nbsimu_min,
		int	nbsimu_quant,
		int	niter_max,
		VectorInt &	cols,
		int(*)(int flag_grid, int nDim, int nech, int *nx, double *dx, double *x0, double nonval, double percent, VectorDouble &tab)	func_valid
	)

Perform a set of valid conditional or non-conditional simulations

Returns

Error return code

Parameters

[in]	dbin	Input Db structure (optional)
[in]	dbout	Output Db structure
[in]	model	Model structure
[in]	neigh	ANeigh structure (optional)
[in]	seed	Seed for random number generator
[in]	nbtuba	Number of turning bands
[in]	nbsimu_min	Minimum number of simulations
[in]	nbsimu_quant	Additional quantum of simulations
[in]	niter_max	Maximum number of iterations
[in]	cols	Vector of column indices
[in]	func_valid	Testing function

Remarks

This function calls a simulation outcome.

It provides a return code:

-1: the simulation outcome is not valid and the process must be

interrupted gently with no error (case of non-convergence)

0: the simulation outcome is not valid

1: the simulation outcome is valid and must be kept

2: the simulation outcome is valid and its returned version

must be kept (this is usefull if func_valid() has modified it).

The func_valid prototype has the following arguments:

• ndim Space dimension
• nx Array of number of grid meshes along all space direction
• dx Array of grid mesh along all space direction
• x0 Array of grid origin along all space direction
• nonval Value corresponding to a missing grid value
• percent Percentage of the simulations already validated
• tab Array containing the simulated outcome

The following lines give an example of func_valid() which considers

a simulation outcome as valid if more than 50% of the valid samples

have a positive value.

 int func_valid(int flag_grid,int ndim,int nech,
                int *nx,double *dx,double *x0,
                double nonval, double percent, double *tab)
{
  double ratio;
  int i,npositive,nvalid;
 
  for (i=0; i<nech; i++)
    {
       if (tab[i] == nonval) continue;
       nvalid++;
       if (tab[i] > 10.) npositive++;
    }
    ratio = (nvalid > 0) ? npositive / nvalid : 0.;
    return(ratio > 0.5);
 }

simu_func_categorical_scale()

void simu_func_categorical_scale	(	Db *	db,
		int	verbose,
		int	nbsimu
	)

Scaling function for the Modification categorical case

Parameters

[in]	db	Db structure
[in]	verbose	1 for the verbose flag
[in]	nbsimu	Number of simulations

simu_func_categorical_transf()

void simu_func_categorical_transf	(	Db *	db,
		int	verbose,
		int	isimu,
		int	nbsimu
	)

Transformation function for the Modification categorical case

Parameters

[in]	db	Db structure
[in]	verbose	1 for the verbose flag
[in]	isimu	Rank of the current simulation
[in]	nbsimu	Number of simulations

simu_func_categorical_update()

void simu_func_categorical_update	(	Db *	db,
		int	verbose,
		int	isimu,
		int	nbsimu
	)

Updating function for the Modification categorical case

Parameters

[in]	db	Db structure
[in]	verbose	1 for the verbose flag
[in]	isimu	Rank of the current simulation
[in]	nbsimu	Number of simulations (stored)

simu_func_continuous_scale()

void simu_func_continuous_scale	(	Db *	db,
		int	verbose,
		int	nbsimu
	)

Scaling function for the Modification continuous case

Parameters

[in]	db	Db structure
[in]	verbose	1 for the verbose flag
[in]	nbsimu	Number of simulations

simu_func_continuous_update()

void simu_func_continuous_update	(	Db *	db,
		int	verbose,
		int	isimu,
		int	nbsimu
	)

Updating function for the Modification continuous case

Parameters

[in]	db	Db structure
[in]	verbose	1 for the verbose flag
[in]	isimu	Rank of the current simulation
[in]	nbsimu	Number of simulations (stored)