#include "geoslib_f.h"
#include "geoslib_f_private.h"
#include "geoslib_old_f.h"
#include "geoslib_enum.h"
#include "Enum/ELoadBy.hpp"
#include "Matrix/MatrixFactory.hpp"
#include "Matrix/MatrixSquareGeneral.hpp"
#include "Matrix/MatrixSparse.hpp"
#include "Matrix/LinkMatrixSparse.hpp"
#include "Matrix/NF_Triplet.hpp"
#include "LinearOp/ProjMatrix.hpp"
#include "Model/NoStatArray.hpp"
#include "Mesh/MeshEStandard.hpp"
#include "Covariances/CovAniso.hpp"
#include "Covariances/ACovAnisoList.hpp"
#include "Covariances/CovFactory.hpp"
#include "Basic/AException.hpp"
#include "Basic/Utilities.hpp"
#include "Basic/Law.hpp"
#include "Basic/MathFunc.hpp"
#include "Basic/File.hpp"
#include "Basic/String.hpp"
#include "Basic/VectorHelper.hpp"
#include "Db/Db.hpp"
#include "Model/Model.hpp"
#include "Mesh/AMesh.hpp"
#include "Mesh/MeshETurbo.hpp"
#include "Calculators/CalcMigrate.hpp"
#include "Space/SpaceSN.hpp"
#include "Geometry/GeometryHelper.hpp"
#include <math.h>
#include <string.h>
#include "csparse_f.h"

Macros
#define	NBLIN_TERMS 10

#define	SPDE_MAX_NGRF 2

Functions
void	spde_option_update (SPDE_Option &s_option, const String &triswitch)

SPDE_Option	spde_option_alloc (void)

SPDE_Matelem &	spde_get_current_matelem (int icov)

void	simu_define_func_transf (void(st_simu_transf)(Db , int, int, int))

void	simu_define_func_update (void(st_simu_update)(Db , int, int, int))

void	simu_define_func_scale (void(st_simu_scale)(Db , int, int))

QChol *	qchol_manage (int mode, QChol *QC)

double	spde_compute_correc (int ndim, double param)

int	spde_attach_model (Model *model)

int	spde_build_stdev (double *vcur)

double *	_spde_get_mesh_dimension (AMesh *amesh)

MatrixSparse *	_spde_fill_S (AMesh amesh, Model model, double *units)

VectorDouble	_spde_fill_TildeC (AMesh amesh, double units)

VectorDouble	_spde_fill_Lambda (Model model, AMesh amesh, const VectorDouble &TildeC)

MatrixSparse *	_spde_build_Q (MatrixSparse S, const VectorDouble &Lambda, int nblin, double blin)

int	spde_build_matrices (Model *model, int verbose)

int	spde_chebychev_operate (MatrixSparse S, Cheb_Elem cheb_elem, const VectorDouble &lambda, const VectorDouble &x, VectorDouble &y)

Cheb_Elem *	spde_cheb_manage (int mode, int verbose, double power, const VectorDouble &blin, MatrixSparse S, Cheb_Elem cheb_old)

Cheb_Elem *	_spde_cheb_duplicate (Cheb_Elem *cheb_in)

int	spde_process (Db dbin, Db dbout, SPDE_Option &s_option, int nbsimu, int ngibbs_burn, int ngibbs_iter, int ngibbs_int)

int	spde_external_copy (SPDE_Matelem &matelem, int icov0)

AMesh *	spde_mesh_load (Db dbin, Db dbout, const VectorDouble &gext, SPDE_Option &s_option, bool verbose)

void	spde_external_mesh_define (int icov0, AMesh *mesh)

void	spde_external_mesh_undefine (int icov0)

MatrixSparse *	spde_external_A_define (int icov0, MatrixSparse *A)

MatrixSparse *	spde_external_A_undefine (int icov0)

MatrixSparse *	spde_external_Q_define (int icov0, MatrixSparse *Q)

MatrixSparse *	spde_external_Q_undefine (int icov0)

void	spde_mesh_assign (AMesh *amesh, int ndim, int ncorner, int nvertex, int nmesh, const VectorInt &arg_meshes, const VectorDouble &arg_points, int verbose)

int	spde_prepar (Db dbin, Db dbout, const VectorDouble &gext, SPDE_Option &s_option)

int	spde_posterior ()

void	spde_free_all (void)

int	spde_check (const Db dbin, const Db dbout, Model model1, Model model2, bool verbose, const VectorDouble &gext, bool mesh_dbin, bool mesh_dbout, bool flag_advanced, bool flag_est, bool flag_std, bool flag_gibbs, bool flag_modif)

int	kriging2D_spde (Db dbin, Model model, SPDE_Option &s_option, int verbose, int nmesh_arg, int nvertex_arg, VectorInt &meshes_arg, VectorDouble &points_arg)

int	spde_eval (const VectorDouble &blin, MatrixSparse *S, const VectorDouble &Lambda, const VectorDouble &TildeC, double power, VectorDouble &x, VectorDouble &y)

MatrixSparse *	db_mesh_neigh (const Db db, AMesh amesh, double radius, int flag_exact, bool, int nactive_arg, int *ranks_arg)

int	m2d_gibbs_spde (Db dbin, Db dbout, Model *model, int flag_ed, int nlayer, int niter, int seed, int nbsimu, int icol_pinch, int flag_drift, int flag_ce, int flag_cstd, int verbose)

Macro Definition Documentation

◆ NBLIN_TERMS

#define NBLIN_TERMS 10

◆ SPDE_MAX_NGRF

#define SPDE_MAX_NGRF 2

Function Documentation

◆ _spde_build_Q()

MatrixSparse* _spde_build_Q	(	MatrixSparse *	S,
		const VectorDouble &	Lambda,
		int	nblin,
		double *	blin
	)

Construct the final sparse matrix Q from the Model

Returns: Error return code

Parameters

[in]	S	Shift operator
[in]	Lambda	Lambda vector
[in]	nblin	Number of blin coefficients
[in]	blin	Array of coefficients for Linear combinaison

◆ _spde_cheb_duplicate()

Cheb_Elem* _spde_cheb_duplicate ( Cheb_Elem * cheb_in )

Duplicate a Cheb_Elem structure

Returns: Pointer to the newly created Cheb_Eleme structure

Parameters

[in] cheb_in Input Cheb_Eleme structure

◆ _spde_fill_Lambda()

VectorDouble _spde_fill_Lambda	(	Model *	model,
		AMesh *	amesh,
		const VectorDouble &	TildeC
	)

Fill the vector for sill correction factors Works for both stationary and non-stationary cases

Parameters

[in]	model	Model structure
[in]	amesh	MeshEStandard structure
[in]	TildeC	Vector TildeC

◆ _spde_fill_S()

MatrixSparse* _spde_fill_S	(	AMesh *	amesh,
		Model *	model,
		double *	units
	)

Fill the sparse matrix S linked to mesh vertices

Returns: G sparse matrix

Parameters

[in]	amesh	MeshEStandard_Mesh structure
[in]	model	Model structure
[in]	units	Array containing the mesh dimensions

◆ _spde_fill_TildeC()

VectorDouble _spde_fill_TildeC	(	AMesh *	amesh,
		double *	units
	)

Fill the vector TildeC (Dimension: nvertex)

Returns: Error return code

Parameters

[in]	amesh	MeshEStandard structure
[in]	units	Array containing the element units

◆ _spde_get_mesh_dimension()

double* _spde_get_mesh_dimension ( AMesh * amesh )

Calculate the array of dimensions of the meshes

Returns: Pointer to the newly allocated array containing mesh dimensions

Parameters

[in] amesh MeshEStandard structure

Remarks: The array returned by this function must be deallocated

◆ db_mesh_neigh()

MatrixSparse* db_mesh_neigh	(	const Db *	db,
		AMesh *	amesh,
		double	radius,
		int	flag_exact,
		bool	,
		int *	nactive_arg,
		int **	ranks_arg
	)

Returns the projection matrix of a set of points (contained in a Db) onto a meshing

Returns: Pointer to the newly created sparse matrix (or NULL)

Parameters

[in]	db	Db structure
[in]	amesh	AMesh structure
[in]	flag_exact	Type of test for intersection (See remarks)
[in]	radius	Neighborhood radius
[out]	nactive_arg	Number of active samples from the Db
[out]	ranks_arg	Ranks of the active samples

Remarks: The calling function must free the argument 'ranks'; When flag_exact is TRUE, for each active sample of Db, a vertex; of the mesh is active as soon as it lies within the vicinity; of the sample.; If flag_exact is FALSE, all vertices of a mesh are considered as; active as soon as the mesh intersects the ball around a sample.; The vicinity is defined as any point located at a distance; from the sample smaller than 'radius'. The distance is calculated; as the Euclidean distance over the space whose dimension is; if the smallest value between the Db et Mesh space dimensions.

◆ kriging2D_spde()

int kriging2D_spde	(	Db *	dbin,
		Model *	model,
		SPDE_Option &	s_option,
		int	verbose,
		int *	nmesh_arg,
		int *	nvertex_arg,
		VectorInt &	meshes_arg,
		VectorDouble &	points_arg
	)

Perform Kriging using SPDE on the set of constructed 2-D triangles

Returns: Error return code

Parameters

[in]	dbin	Db input structure
[in]	model	Model structure
[in]	s_option	SPDE_Option structure
[in]	verbose	Verbose option
[out]	nmesh_arg	Number of triangles generated
[out]	nvertex_arg	Number of vertices
[out]	meshes_arg	Array of triangulate vertex indices
[out]	points_arg	Array containing the 2-D vertices

Remarks: It is the responsability of the calling function to free; the returned arrays 'points_arg' and 'meshes_arg'.; This function assumes that there is only one Model defined; Hence the test on the number of models

◆ m2d_gibbs_spde()

int m2d_gibbs_spde	(	Db *	dbin,
		Db *	dbout,
		Model *	model,
		int	flag_ed,
		int	nlayer,
		int	niter,
		int	seed,
		int	nbsimu,
		int	icol_pinch,
		int	flag_drift,
		int	flag_ce,
		int	flag_cstd,
		int	verbose
	)

Perform Gibbs on a multilayer setup

Returns: Error return code

Parameters

[in]	dbin	Db input structure
[in]	dbout	Db output structure
[in]	model	Model structure
[in]	flag_ed	1 if External Drit is used
[in]	nlayer	Number of layers
[in]	niter	Number of iterations
[in]	seed	Seed for random number generator
[in]	nbsimu	Number of simulaations
[in]	icol_pinch	Address of the variable containing the pinchout
[in]	flag_drift	1 to return the drift only 0 the simulations
[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 option

Remarks: In 'dbin':; - the lower and upper bounds must be defined for each datum; (set to the locator ELoc::L and ELoc::U; In 'dbout':; - the trend (if flag_ed is 1) must be defined and set to; the locator ELoc::F; When defined, the pinchout should be defined as a grid variable; with values ranging between 0 and 1 (FFFF are admitted).; It will serve as a multiplier to the Mean thickness maps.

◆ qchol_manage()

QChol* qchol_manage	(	int	mode,
		QChol *	QC
	)

Manage the Sparse matrix and its cholesky decomposition

Returns: Pointer to the QChol structure

Parameters

[in]	mode	Management mode 1 : Allocation -1 : Total deallocation (NULL is returned)
[in]	QC	Pointer to the QChol structure (when mode < 0)

◆ simu_define_func_scale()

void simu_define_func_scale ( void(*)(Db *, int, int) st_simu_scale )

Define the function to scale the Modification arrays

Parameters

[in] st_simu_scale Pointer to the scaling function

◆ simu_define_func_transf()

void simu_define_func_transf ( void(*)(Db *, int, int, int) st_simu_transf )

Define the function to transform a simulation

Parameters

[in] st_simu_transf Pointer to the transformation function

◆ simu_define_func_update()

void simu_define_func_update ( void(*)(Db *, int, int, int) st_simu_update )

Define the function to account for the current simulation outcome in the calculation of the Modification arrays

Parameters

[in] st_simu_update Pointer to the update function

◆ spde_attach_model()

int spde_attach_model ( Model * model )

Attach the model

Returns: Error returned code

Parameters

[in] model Model structure

◆ spde_build_matrices()

int spde_build_matrices	(	Model *	model,
		int	verbose
	)

Build all matrices needed for establishing the Q sparse matrix

Returns: Error return code

Parameters

[in]	model	Model structure
[in]	verbose	Verbose option

Remarks: Contents of SP_MAT (sparse matrices or vectors) is allocated here; It must be freed by the calling functions

◆ spde_build_stdev()

int spde_build_stdev ( double * vcur )

Perform the calculation of the Standard Deviation of Estimation Error

Returns: Error return code

Parameters

[out] vcur Output array

◆ spde_cheb_manage()

Cheb_Elem* spde_cheb_manage	(	int	mode,
		int	verbose,
		double	power,
		const VectorDouble &	blin,
		MatrixSparse *	S,
		Cheb_Elem *	cheb_old
	)

Manage Cheb_Elem structure

Returns: Error return code

Parameters

[in]	mode	1 for allocation; -1 for deallocation
[in]	verbose	Verbose flag
[in]	power	Parameter passed to Chebychev function
[in]	blin	Array of coefficients for Linear combinaison
[in]	S	Shift operator
[in]	cheb_old	Cheb_Elem to be freed (only for mode=-1)

Remarks: Arguments 'power', 'nblin', 'blin' and 'B' are used if mode=1; Argument 'cheb_old' is used if mode=-1

◆ spde_chebychev_operate()

int spde_chebychev_operate	(	MatrixSparse *	S,
		Cheb_Elem *	cheb_elem,
		const VectorDouble &	lambda,
		const VectorDouble &	x,
		VectorDouble &	y
	)

Perform the Chebychev polynomial procedure on an input vector

Returns: Error return code

Parameters

[in]	S	Shift operator
[in]	cheb_elem	Cheb_Elem structure
[in]	lambda	Scaling vector
[in]	x	Input array (Dimension: nvertex)
[out]	y	Output array (Dimension: nvertex)

◆ spde_check()

int spde_check	(	const Db *	dbin,
		const Db *	dbout,
		Model *	model1,
		Model *	model2,
		bool	verbose,
		const VectorDouble &	gext,
		bool	mesh_dbin,
		bool	mesh_dbout,
		bool	flag_advanced,
		bool	flag_est,
		bool	flag_std,
		bool	flag_gibbs,
		bool	flag_modif
	)

Define the main options

Returns: Error return code

Parameters

[in]	dbin	Pointer to the input Db
[in]	dbout	Pointer to the output Db
[in]	model1	Model structure (first)
[in]	model2	Model structure (second)
[in]	verbose	Verbose flag
[in]	gext	Array of domain dilation
[in]	mesh_dbin	True if Input points must participate to meshing
[in]	mesh_dbout	True if Output points must participate to meshing
[in]	flag_advanced	True for advanced calculus (estimation or simulation) False if only matrices are required
[in]	flag_est	True for estimation
[in]	flag_std	True for standard deviation
[in]	flag_gibbs	True for Gibbs sampler
[in]	flag_modif	True for post-processing simulations

Remarks: This function initiates the Matelem structures

◆ spde_compute_correc()

double spde_compute_correc	(	int	ndim,
		double	param
	)

◆ spde_eval()

int spde_eval	(	const VectorDouble &	blin,
		MatrixSparse *	S,
		const VectorDouble &	Lambda,
		const VectorDouble &	TildeC,
		double	power,
		VectorDouble &	x,
		VectorDouble &	y
	)

Perform the product of a vector by the inverse of the power of a sparse matrix using the Chebychev Polynomial procedure

Returns: Error return code

Parameters

[in]	blin	Array of coefficients for Linear combinaison
[in]	S	Shift operator
[in]	Lambda	Vector Lambda
[in]	TildeC	Vector TildeC
[in]	power	Parameter used in the Chebychev approximation
[in]	x	Input array
[out]	y	Output array

◆ spde_external_A_define()

MatrixSparse* spde_external_A_define	(	int	icov0,
		MatrixSparse *	A
	)

◆ spde_external_A_undefine()

MatrixSparse* spde_external_A_undefine ( int icov0 )

◆ spde_external_copy()

int spde_external_copy	(	SPDE_Matelem &	matelem,
		int	icov0
	)

Copy the contents of the internal S_EXTERNAL_AQ into an output Matelem

Error return code

Parameters

[in]	matelem	Output SPDE_Matelem structure
[in]	icov0	Rank of the current Covariance

◆ spde_external_mesh_define()

void spde_external_mesh_define	(	int	icov0,
		AMesh *	mesh
	)

Manage the contents of the External AMesh structure used for storing external meshing information

Parameters

[in]	icov0	Rank of the current covariance (from 0 to 2)
[in]	mesh	AMesh pointer

◆ spde_external_mesh_undefine()

void spde_external_mesh_undefine ( int icov0 )

◆ spde_external_Q_define()

MatrixSparse* spde_external_Q_define	(	int	icov0,
		MatrixSparse *	Q
	)

◆ spde_external_Q_undefine()

MatrixSparse* spde_external_Q_undefine ( int icov0 )

◆ spde_free_all()

void spde_free_all ( void )

Free all memory used in SPDE

◆ spde_get_current_matelem()

SPDE_Matelem& spde_get_current_matelem ( int icov )

Get the pointer to the current SPDE_Matelem structure

Parameters

[in] icov Rank of the target Covariance (or -1)typedef struct

◆ spde_mesh_assign()

void spde_mesh_assign	(	AMesh *	amesh,
		int	ndim,
		int	ncorner,
		int	nvertex,
		int	nmesh,
		const VectorInt &	arg_meshes,
		const VectorDouble &	arg_points,
		int	verbose
	)

Assign fields of the AMesh structure which would have been calculated elsewhere

Parameters

[in]	ndim	Space dimension
[in]	ncorner	Number of vertices per element
[in]	nvertex	Number of points
[in]	nmesh	Number of meshes
[in]	arg_meshes	Array containing the meshes
[in]	arg_points	Array containing the vertex coordinates
[in]	verbose	Verbose flag
[in,out]	amesh	Pointer to AMesh to be assigned

◆ spde_mesh_load()

AMesh* spde_mesh_load	(	Db *	dbin,
		Db *	dbout,
		const VectorDouble &	gext,
		SPDE_Option &	s_option,
		bool	verbose
	)

Load the AMesh structure

Returns: Pointer on the newly allocated AMesh

Parameters

[in]	dbin	Db structure for the conditioning data
[in]	dbout	Db structure of the grid
[in]	gext	Array of domain dilation
[in]	s_option	SPDE_Option structure
[in]	verbose	Verbose option

◆ spde_option_alloc()

SPDE_Option spde_option_alloc ( void )

Manage the SPDE_Option structure

◆ spde_option_update()

void spde_option_update	(	SPDE_Option &	s_option,
		const String &	triswitch
	)

Add a new option item for an additional covariance structure

Parameters

[in]	s_option	Pointer to SPDE_Option to be freed (if mode<0)
[in]	triswitch	String defining the meshing characteristics

◆ spde_posterior()

int spde_posterior ( )

Cleaning operation after SPDE

Returns: Error return code

◆ spde_prepar()

int spde_prepar	(	Db *	dbin,
		Db *	dbout,
		const VectorDouble &	gext,
		SPDE_Option &	s_option
	)

Preparation using SPDE (for all GRF and COV)

Returns: Error return code

Parameters

[in]	dbin	Db structure for the conditioning data
[in]	dbout	Db structure of the grid
[in]	gext	Array of domain dilation
[in]	s_option	SPDE_Option structure

◆ spde_process()

int spde_process	(	Db *	dbin,
		Db *	dbout,
		SPDE_Option &	s_option,
		int	nbsimu,
		int	ngibbs_burn,
		int	ngibbs_iter,
		int	ngibbs_int
	)

Perform the main Simulations steps after Q has been constructed

Returns: Error return code

Parameters

[in]	dbin	Input Db grid structure (optional)
[in]	dbout	Output Db grid structure
[in]	s_option	SPDE_Option structure
[in]	nbsimu	Number of simulations
[in]	ngibbs_burn	Number of iterations (Burning step)
[in]	ngibbs_iter	Number of iterations
[in]	ngibbs_int	Number of iterations internal to Gibbs (SPDE)

Remarks: If the number of simulations 'nbsimu' is set to 0,; the simulation algorithm is turned into a kriging one

Macros

Functions

Macro Definition Documentation

◆ NBLIN_TERMS

◆ SPDE_MAX_NGRF

Function Documentation

◆ _spde_build_Q()

◆ _spde_cheb_duplicate()

◆ _spde_fill_Lambda()

◆ _spde_fill_S()

◆ _spde_fill_TildeC()

◆ _spde_get_mesh_dimension()

◆ db_mesh_neigh()

◆ kriging2D_spde()

◆ m2d_gibbs_spde()

◆ qchol_manage()

◆ simu_define_func_scale()

◆ simu_define_func_transf()

◆ simu_define_func_update()

◆ spde_attach_model()

◆ spde_build_matrices()

◆ spde_build_stdev()

◆ spde_cheb_manage()

◆ spde_chebychev_operate()

◆ spde_check()

◆ spde_compute_correc()

◆ spde_eval()

◆ spde_external_A_define()

◆ spde_external_A_undefine()

◆ spde_external_copy()

◆ spde_external_mesh_define()

◆ spde_external_mesh_undefine()

◆ spde_external_Q_define()

◆ spde_external_Q_undefine()

◆ spde_free_all()

◆ spde_get_current_matelem()

◆ spde_mesh_assign()

◆ spde_mesh_load()

◆ spde_option_alloc()

◆ spde_option_update()

◆ spde_posterior()

◆ spde_prepar()

◆ spde_process()