1.2.1
CCC
 
spde.cpp File Reference
#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_Matelemspde_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))
 
QCholqchol_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_Elemspde_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)
 
AMeshspde_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)
 
MatrixSparsespde_external_A_define (int icov0, MatrixSparse *A)
 
MatrixSparsespde_external_A_undefine (int icov0)
 
MatrixSparsespde_external_Q_define (int icov0, MatrixSparse *Q)
 
MatrixSparsespde_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)
 
MatrixSparsedb_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

#define NBLIN_TERMS   10
#define SPDE_MAX_NGRF   2

Function Documentation

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]SShift operator
[in]LambdaLambda vector
[in]nblinNumber of blin coefficients
[in]blinArray of coefficients for Linear combinaison
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_inInput Cheb_Eleme structure
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]modelModel structure
[in]ameshMeshEStandard structure
[in]TildeCVector TildeC
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]ameshMeshEStandard_Mesh structure
[in]modelModel structure
[in]unitsArray containing the mesh dimensions
VectorDouble _spde_fill_TildeC ( AMesh amesh,
double *  units 
)

Fill the vector TildeC (Dimension: nvertex)

Returns
Error return code
Parameters
[in]ameshMeshEStandard structure
[in]unitsArray containing the element units
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]ameshMeshEStandard structure
Remarks
The array returned by this function must be deallocated
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]dbDb structure
[in]ameshAMesh structure
[in]flag_exactType of test for intersection (See remarks)
[in]radiusNeighborhood radius
[out]nactive_argNumber of active samples from the Db
[out]ranks_argRanks 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.
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]dbinDb input structure
[in]modelModel structure
[in]s_optionSPDE_Option structure
[in]verboseVerbose option
[out]nmesh_argNumber of triangles generated
[out]nvertex_argNumber of vertices
[out]meshes_argArray of triangulate vertex indices
[out]points_argArray 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
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]dbinDb input structure
[in]dboutDb output structure
[in]modelModel structure
[in]flag_ed1 if External Drit is used
[in]nlayerNumber of layers
[in]niterNumber of iterations
[in]seedSeed for random number generator
[in]nbsimuNumber of simulaations
[in]icol_pinchAddress of the variable containing the pinchout
[in]flag_drift1 to return the drift only 0 the simulations
[in]flag_ce1 if the conditional expectation should be returned instead of simulations
[in]flag_cstd1 if the conditional standard deviation should be returned instead of simulations
[in]verboseVerbose 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* qchol_manage ( int  mode,
QChol QC 
)

Manage the Sparse matrix and its cholesky decomposition

Returns
Pointer to the QChol structure
Parameters
[in]modeManagement mode 1 : Allocation -1 : Total deallocation (NULL is returned)
[in]QCPointer to the QChol structure (when mode < 0)
void simu_define_func_scale ( void(*)(Db *, int, int)  st_simu_scale)

Define the function to scale the Modification arrays

Parameters
[in]st_simu_scalePointer to the scaling function
void simu_define_func_transf ( void(*)(Db *, int, int, int)  st_simu_transf)

Define the function to transform a simulation

Parameters
[in]st_simu_transfPointer to the transformation function
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_updatePointer to the update function
int spde_attach_model ( Model model)

Attach the model

Returns
Error returned code
Parameters
[in]modelModel structure
int spde_build_matrices ( Model model,
int  verbose 
)

Build all matrices needed for establishing the Q sparse matrix

Returns
Error return code
Parameters
[in]modelModel structure
[in]verboseVerbose option
Remarks
Contents of SP_MAT (sparse matrices or vectors) is allocated here
It must be freed by the calling functions
int spde_build_stdev ( double *  vcur)

Perform the calculation of the Standard Deviation of Estimation Error

Returns
Error return code
Parameters
[out]vcurOutput array
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]mode1 for allocation; -1 for deallocation
[in]verboseVerbose flag
[in]powerParameter passed to Chebychev function
[in]blinArray of coefficients for Linear combinaison
[in]SShift operator
[in]cheb_oldCheb_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
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]SShift operator
[in]cheb_elemCheb_Elem structure
[in]lambdaScaling vector
[in]xInput array (Dimension: nvertex)
[out]yOutput array (Dimension: nvertex)
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]dbinPointer to the input Db
[in]dboutPointer to the output Db
[in]model1Model structure (first)
[in]model2Model structure (second)
[in]verboseVerbose flag
[in]gextArray of domain dilation
[in]mesh_dbinTrue if Input points must participate to meshing
[in]mesh_dboutTrue if Output points must participate to meshing
[in]flag_advancedTrue for advanced calculus (estimation or simulation) False if only matrices are required
[in]flag_estTrue for estimation
[in]flag_stdTrue for standard deviation
[in]flag_gibbsTrue for Gibbs sampler
[in]flag_modifTrue for post-processing simulations
Remarks
This function initiates the Matelem structures
double spde_compute_correc ( int  ndim,
double  param 
)
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]blinArray of coefficients for Linear combinaison
[in]SShift operator
[in]LambdaVector Lambda
[in]TildeCVector TildeC
[in]powerParameter used in the Chebychev approximation
[in]xInput array
[out]yOutput array
MatrixSparse* spde_external_A_define ( int  icov0,
MatrixSparse A 
)
MatrixSparse* spde_external_A_undefine ( int  icov0)
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]matelemOutput SPDE_Matelem structure
[in]icov0Rank of the current Covariance
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]icov0Rank of the current covariance (from 0 to 2)
[in]meshAMesh pointer
void spde_external_mesh_undefine ( int  icov0)
MatrixSparse* spde_external_Q_define ( int  icov0,
MatrixSparse Q 
)
MatrixSparse* spde_external_Q_undefine ( int  icov0)
void spde_free_all ( void  )

Free all memory used in SPDE

SPDE_Matelem& spde_get_current_matelem ( int  icov)

Get the pointer to the current SPDE_Matelem structure

Parameters
[in]icovRank of the target Covariance (or -1)typedef struct
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]ndimSpace dimension
[in]ncornerNumber of vertices per element
[in]nvertexNumber of points
[in]nmeshNumber of meshes
[in]arg_meshesArray containing the meshes
[in]arg_pointsArray containing the vertex coordinates
[in]verboseVerbose flag
[in,out]ameshPointer to AMesh to be assigned
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]dbinDb structure for the conditioning data
[in]dboutDb structure of the grid
[in]gextArray of domain dilation
[in]s_optionSPDE_Option structure
[in]verboseVerbose option
SPDE_Option spde_option_alloc ( void  )

Manage the SPDE_Option structure

void spde_option_update ( SPDE_Option s_option,
const String triswitch 
)

Add a new option item for an additional covariance structure

Parameters
[in]s_optionPointer to SPDE_Option to be freed (if mode<0)
[in]triswitchString defining the meshing characteristics
int spde_posterior ( )

Cleaning operation after SPDE

Returns
Error return code
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]dbinDb structure for the conditioning data
[in]dboutDb structure of the grid
[in]gextArray of domain dilation
[in]s_optionSPDE_Option structure
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]dbinInput Db grid structure (optional)
[in]dboutOutput Db grid structure
[in]s_optionSPDE_Option structure
[in]nbsimuNumber of simulations
[in]ngibbs_burnNumber of iterations (Burning step)
[in]ngibbs_iterNumber of iterations
[in]ngibbs_intNumber 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