#ifndef RFsimu_H #define RFsimu_H 1 // #define RFBOOL 1 #include "RFsimu.public.h" #include "GSLvsR.h" #include // global definitions used: // RF_FLOAT, RF_GSL #ifdef RF_FLOAT // we use Real instead of double or float, to get // flexibility (double->precision, float:reduced memory allocation) // R needs compilation with Real==double // tolerances are defined correpondingly #define EPSILON 0.0000001 #define EPSILON1000 0.0001 #else #define EPSILON 0.00000000001 #define EPSILON1000 0.00000001 #endif #define METHODMAXCHAR 31 /* max number of character to describe a Method (including \0)*/ #define COVMAXCHAR 14 /* max number of characters for covariance name*/ #define MAXDIM 3 /* if explicite simulation method is not given by user, the programme will first try a methods according to "first"; 1 then it tries the other methods according to a default list: (CircEmbed only if the simulation is on a grid) 1 dim: CircEmbed, Direct 2 dim: CircEmbed, TBM2, spectralTBM, TBM3, AdditiveMpp, Direct 3 dim: CircEmbed, TBM3, Direct the method actually used is stored in actually_used. */ #define MAXERRORSTRING 100 #define MAXKEYS 10 #define MAXNRCOVFCTS 30 #define GetNotPrint 1 /* if parameters are questionned, then also printed if =0 */ #define TWOPI 6.2831853071795864769252867666 #define INVSQRTTWO 0.70710678118654752440084436210 #define INVSQRTTWOPI 0.39894228040143270286 #define SQRTTWOPI 2.5066282746310002416 #define SQRTPI 1.7724538509055158819 #define SQRTTWO 1.414213562373095048801688724209698 #define INVLOG2 1.442695040888963 //#define TWOSQRTPI 3.5449077018110317638 // Error codes & messages // positive values are error codes // negative values are messages #define USEOLDINIT -1 /* not an error, just a message !*/ #define NOERROR 0 #define ERRORNOTDEFINED 1 /* the specification for the covariance and method is not given/known, e.g. TBM2 for many covariance functions */ #define ERRORMETHODNOTALLOWED 2 /* wrong method for the specified covariance function or grid */ #define ERRORNOTPROGRAMMED 3 #define ERRORCOVNOTALLOWED 4 /* wrong dimension for the specified covariance function */ #define ERRORFAILED 5 /* method didn't work for the specified parameters */ #define ERRORMEMORYALLOCATION 6 /* malloc returned NULL pointer */ #define ERRORNOTINITIALIZED 7 /* key.active==false in DoSimulateRF; is only checked there !!! */ #define ERRORKEYNROUTOFRANGE 8 /* wrong keynr */ #define ERRORDECOMPOSITION 9 /* direct simulation; matrix inversion failed */ #define ERRORPRECISION 10 /* direct simulation; matrix inversion didn't reach the desired presicion */ #define ERRORRESCALING 11 /* attempt to rescale where not possible */ #define ERRORVARIOGRAMONLY 12 /* attempt to get the covariance whereas only the variogram is defined */ #define ERRORFOURIER 13 /* specific error in Fourier transformation */ #define ERRORCOVFAILED 14 /* covariance model not defined for specified parameters */ #define ERRORREGISTER 15 /* wrong register number */ #define ERRORCOORDINATES 16 /* coordinates are not given or invalid grid specification */ #define ERRORNEGATIVEVAR 17 #define ERRORPARAMNUMBER 18 #define ERRORDIM 19 /* dim<1 or dim>MAXDIM */ #define ERRORNEGATIVESCALE 20 #define ERRORWAVING 21 /* numerical determination of practical range fails */ #define ERRORWRONGINIT 22 /* dosimulate and initsimulate do not match */ #define ERRORNN 23 /* nn in TBM too large */ #define ERRORDUMMY 99 /* no error, only for tracing purposes */ #define ERRORSILLNULL 101 /* extremes : sill must be positive */ #define ERRORVARMEAN 201 /* Poisson distr.: mean<>variance */ /* parameters: */ #define MEAN 0 #define VARIANCE 1 /* without NUGGET */ #define NUGGET 2 #define SCALE 3 #define KAPPA 4 /* Kappas should be always the last ones of transferred parameters */ #define KAPPAX 5 #define KAPPAXX 6 #define LASTKAPPA 7 #define INVSCALE 8 /* 1.0 / SCALE, calculated by InitSimulateRF */ #define SILL 9 /* VARIANCE + NUGGET */ #define TOTAL_PARAM 10 /* check CheckCovariance, InitSimulateRF (both C & R), and, hyperbolic, TBM3hyperbolic, twodimfractalbrownian R if changed! */ /* Natural Scaling -- in contrast to the RFparameters(PracticalRange) GENERAL_NATURALSCALING may take the following values (furthermore, 0: no rescaling, +10 numerical determination allowed */ #define NATSCALE_EXACT 1 #define NATSCALE_APPROX 2 #define NATSCALE_MLE 3 /* check mleRF when changing !! */ /* check with InitSimulateRF in case of any changes !! both InitSimulateRF and DoSimulateRF */ #define DISTR_GAUSS 0 #define DISTR_POISSON 1 #define DISTR_MAXSTABLE 2 #define XSTART 0 #define XEND 1 #define XSTEP 2 extern int currentNrCov; extern char ERRORSTRING_OK[MAXERRORSTRING],ERRORSTRING_WRONG[MAXERRORSTRING]; typedef enum InversionMethod { Cholesky, SVD, NoFurtherInversionMethod } InversionMethod; extern char METHODNAMES[][METHODMAXCHAR]; typedef enum SimulationType { CircEmbed, /* Circulant embedding */ CircEmbedLocal,/* Circulant embedding on an enlarged grid; (fract. Brown.), not implemented yet */ TBM2, TBM3, /* Turning Bands, performed in 2 or 3 dimensions */ SpectralTBM, /* Spectral turning bands, only 2 dimensional */ Direct, /* directly, by matrix inversion */ Nugget, /* just simulate independent variables */ AdditiveMpp, /* "additive MPP (random coins)", will be only in the plane*/ Hyperplane, /* by superposition of Poisson hyperplane tessellations; only in the plane! Not implemented yet*/ Special, /* Not implemented yet, any other (exoctic) method not included in the above ones */ Nothing, /* must always be the last of the list of methods for simulation Gaussian random fields */ MaxMpp, Forbidden /* must always be the last one */ } SimulationType; typedef enum action_type{Store,Storing,UseOldOne,NoAction} action_type; typedef struct cov_fct; // in the following structure, user specifications and the intermediate results // by the specific init_procedures is stored typedef void (*destructor)(void **); typedef struct key_type { bool grid, /* simulation on a grid ? */ active, /* has the init_procedure been called successfully? */ nuggetincluded; /* is the method (see below) able to simulate an additive nugget effect? */ SimulationType method, /* the current method (out of SimulationsType) which is tried at the moment or which has been successfully initialized */ tbm_method; /* if method==tbm2/3, then the method for the line is specified, this method can only be chosen by the programmer of the covariance functions, not by the user (except by global settings)! maybe, a search algorithm among all 1-dim simulation methods should be implemented, in case the programmers specification fails */ Real *x[MAXDIM], /* the coordinates */ param[TOTAL_PARAM]; /* parameters */ void *S; /* pointer to intermediate results; delicate : either ==NULL or should be a genuine pointer ! no hanging pointers allowed! */ destructor destruct; /* function called to delete intermediate results */ void *SX; /* pointer to intermediate results, second type, like extremes -- not used yet */ destructor destructX; cov_fct *cov; /* pointer to CovList[covnr], where details about the covariance functions can be found */ int covnr, /* number of the covariance function in CovList; obtained by function GetCovFunction */ lx, /* number of elements in *x[0][] */ length[MAXDIM], /* only if grid==true: number of pixels of each direction */ dim; /* dimension of RF */ long totallength; /* number of pixels in total; == lx if grid==false == lengthx * lengthy * lengthz if grid==true */ int naturalscaling; int distribution; } key_type; extern key_type KEY[MAXKEYS]; typedef struct mpp_storage{ // avoid the same names as in addBool_storage, RFbool.h Real integral, integralsq, integralpos, effectiveRadius, effectivearea, maxheight, addradius; Real min[MAXDIM], length[MAXDIM]; Real c[10]; /* local memory for diverse constants, depending on the specific model */ int dim; } mpp_storage; // how a covariance function, a spectral measure, a function specifying the // additive MPP model should look like typedef Real (*covfct)(Real,Real *); /* h,parameters */ typedef Real (*scalefct)(Real *, int); /* h,parameters */ typedef int (*checkfct)(key_type*); /* h,parameters */ typedef Real (*parameterfct)(Real *,Real,Real *); /* ?? h,parameters */ typedef Real (*randommeasure)(Real * END_WITH_GSL); /* parameters, RANDOM */ typedef Real (*mppmodel)(Real *); /* point of random field */ typedef void (*MPPScales)(key_type *); /* integral, integralOFsq, effectiveRadius */ typedef void (*MPPRandom)(mpp_storage *bs, Real *min, Real *max, mppmodel* /* the random function/model, returned */ END_WITH_RANDOM); // how any simulation method should look like, they return the error value typedef int (*generalSimuInit)(key_type*); typedef void (*generalSimuMethod)(key_type*,Real* END_WITH_GSL); /* key, RF, RANDOM */ // specification for simulation method on line for TBM2/3: typedef void (*simu1dim)(key_type*,Real* END_WITH_GSL); /* h,parameters */ // here all the different method for simulating a RF of a specified covariance // function are stored typedef struct cov_fct { char name[COVMAXCHAR]; covfct cov; /* CircEmbed, direct */ scalefct naturalscale; /* CircEmbedLocal, Brownian motion!: cov_loc : local covariance function Real extension_factor(Real *p, Real x, Real *newp): returns the range of the local covariance function (if true one has range one) furthermore, it returns a new parameter set, that is adapted to the local simulation x is usually the largest distance to be simulated square_factor : see Stein, p.3 (1), parameter $c_2$. */ covfct cov_loc; parameterfct extension_factor;scalefct square_factor; covfct cov_tbm2,cov_tbm3; SimulationType tbm_method; randommeasure spectral; /* spectral tbm */ MPPScales add_mpp_scl; MPPRandom add_mpp_rnd; void *hyperplane; /* not programmed yet */ generalSimuInit initother; /* any other exotic way */ generalSimuMethod other; /* which may be given by a programmer */ SimulationType first[MAXDIM]; /* the preferred method for a specific dimension no distinction between grid/no grid -- lack! */ int kappas; /* number of parameters additional to the standard ones, MEAN..SCALE */ checkfct check; // } cov_fct; void ErrorMessage(SimulationType m, int error); extern void InitModelList(); /* absolutely necessary to call this function at the very beginning ! but done automatically if SimulateRF, etc. is called for the first time (direct call necessary if further covariance functions ate to be added by user!) initiating CovList with a dozen standard covariance models */ // the following function enables the programmer to get a new covariance // function stored extern void IncludeModel(char *name, covfct cov, scalefct naturalscale, covfct cov_loc, parameterfct extension_factor, scalefct square_factor, covfct cov_tbm2, covfct cov_tbm3, SimulationType tbm_method, randommeasure spectral, MPPScales add_mpp_scl, MPPRandom add_mpp_rnd, generalSimuInit initother, generalSimuMethod other, int kappas, SimulationType r1,SimulationType r2,SimulationType r3, checkfct check ); extern void IncludeModel(char *name, covfct cov,scalefct naturalscale, covfct cov_loc, Real extension_factor, covfct cov_tbm2, covfct cov_tbm3, SimulationType tbm_method, randommeasure spectral, int kappas, SimulationType r1,SimulationType r2,SimulationType r3, checkfct check ); extern void IncludeModel(char *name, covfct cov, scalefct naturalscale, covfct cov_tbm2, covfct cov_tbm3, SimulationType tbm_method, randommeasure spectral, int kappas, SimulationType r1,SimulationType r2,SimulationType r3, checkfct check ); // all simulation methods have the following structure // int init_foo(key_type *key): // returns errorcode // This function does the most initialization settings, // often time consuming and errors may occur // intermediate results are stored by means the pointer key->S // the function expects that key->S NULL when called // the function is expected to return key->S on error // // void do_foo(key_type *key, Real *res END_WITH_GSLRNG) // returns the simulation result in res ; // only basics checks like assert(RANDOM!=NULL); // any other error may not occur, as this should be put into init_foo // key_active is set to GENERAL_STORING; key->S is freed and set to NULL // if !key_active // see RFspectral.cc int init_simulatespectral(key_type *key); void do_simulatespectral(key_type *key, Real *res END_WITH_GSLRNG); // see RFdirect.cc int init_directGauss(key_type *key); void do_directGauss(key_type *key, Real *res END_WITH_GSLRNG); // see RFcircembed.cc typedef struct FFT_storage { double* work; int *iwork, nseg; } FFT_storage; int fastfourier(double *data, int *m, int dim, bool first, FFT_storage *S); void FFT_destruct(FFT_storage *S); void FFT_NULL(FFT_storage *S); int internal_init_circulantembedding(key_type * key); int init_circulantembedding(key_type * key); void do_circulantembedding(key_type *key, Real *res END_WITH_GSLRNG); int init_circ_embed_local(key_type * key); void do_circ_embed_local(key_type *key, Real *res END_WITH_GSLRNG); // see RFtbm.cc int init_circulantembedding1dim(key_type *key, Real *p); void do_circulantembedding1dim(key_type *key,Real *res END_WITH_GSLRNG); int init_turningbands(key_type *key); void do_turningbands(key_type *key,Real *res END_WITH_GSLRNG); extern SimulationType TBM_METHOD; // see RFbool.cc int init_mpp(key_type * key); void do_addmpp(key_type *key, Real *res END_WITH_GSLRNG); void printKEY(int keyNr); extern int GENERAL_PRINTLEVEL; extern int GENERAL_NATURALSCALING; extern int GENERAL_STORING; extern char ERRORSTRING_OK[MAXERRORSTRING]; extern char ERRORSTRING_WRONG[MAXERRORSTRING]; extern Real MPP_APPROXZERO; #endif