https://github.com/cran/RandomFields
Tip revision: c751cbeafa6661ec449a5a53a8b3e659f558be70 authored by Martin Schlather on 27 May 2005, 00:00:00 UTC
version 1.2.16
version 1.2.16
Tip revision: c751cbe
addownfctns.cc
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/timeb.h>
#include <assert.h>
#include <string.h>
#include "RFsimu.h"
#include "RFCovFcts.h"
#include "MPPFcts.h"
#include <unistd.h>
#include "addownfctns.h"
/* In the following an example is given for an extension of
a current implementation of a covariance function.
If stablelocal and stableS are coded correctly, then
a stable random field can be generated by Stein's
(2002) local circulant embedding method using
GaussRF(...., model="stable", method="local")
Probably, newcheckstable should also be changed.
*/
double XXstablelocal(double *x, double *p, int effectivdim) {
// cf. Stein (2002), "K(r)"
// and twodimfractalBrownianlocal in RFCovFcts.cc
assert(false); // to be replaced by the true code
}
double XXstableS(double *p) {
// cf. Stein (2002), "c_2"
// and twodimfractalBrownianS in RFCovFcts.cc
assert(false); // to be replaced by the true code
}
double CIRCEMBED_MINDIAMETER=0.0;
void SetMinDiam(double *x) {
// to be used for experimental purposes in R by
// .C("SetMinDiam", as.double( *** ))
// if the value is less than or equal to zero, CIRCEMBED_MINDIAMETER
// will be ignored, see below. Note, that here we have interactions
// between fields size, scale parameter and CIRCEMBED_MINDIAMETER
CIRCEMBED_MINDIAMETER = *x;
}
double stableGetTheta(double *p){
double theta;
theta = p[DIAMETER]; //p[DIAMETER] ist der Durchmesser des zu simulierenden
// Feldes in Einheiten des Scalenparameters von
// GaussRF(..., param(*,*,*,scale))
if (CIRCEMBED_MINDIAMETER > 0.0) {
if (theta<CIRCEMBED_MINDIAMETER) theta = CIRCEMBED_MINDIAMETER;
} else {
if (p[KAPPA]<1.5) {
if (theta<1.0) theta=1.0; // derzeitige Mindestanforderung, wird man
// vermutlich abmildern koennen
} else if (p[KAPPA] <= 1.99) {
if (theta<4.0) theta=4.0;
} else assert(false); // noch keine Daten vorhanden,
// diese Weiche ist in newcheckstable abzufangen!
}
return theta;
}
double stablelocal(double *x, double *p, int effectivdim) {
double tha, theta, y;
theta = stableGetTheta(p);
tha = pow(theta, p[KAPPA]);
y = fabs(*x);
return (y <= theta) ?
(tha * p[KAPPA] * (y * y / theta / theta - 1.0) - 2.0) / (2.0*exp(tha))
+ exp(-pow(y, p[KAPPA]))
: 0.0; // irgendwie so aehnlich, bitte ueberpruefen
// da ich in der Theorie nicht so weit drin bin:
// stellt die obige Definition die Varianz 1 des Feldes sicher?
// wenn nein, muss entsprechende korrigiert werden
}
double stableS(double *p) {
double tha;
// cf. Stein (2002), "c_2"
// and twodimfractalBrownianS in RFCovFcts.cc
tha = pow(stableGetTheta(p), p[KAPPA]);
return p[KAPPA] * tha / (2.0 * exp(tha));
}
int newcheckstable(double *param, int timespacedim, SimulationType method){
if ((param[KAPPA]<=0) || (param[KAPPA]>2.0)) {
strcpy(ERRORSTRING_OK,"0<kappa<=2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
if (method==CircEmbedIntrinsic && param[KAPPA]>1.1) {
strcpy(ERRORSTRING_OK,"0<=kappa<=1.1 if method=local CE");
sprintf(ERRORSTRING_WRONG,"%f", param[KAPPA]);
return ERRORCOVFAILED;
}
return 0;
}
void addownfunctions() {
/* compare with InitModelList() in RFinitNerror.cc!
actually the definition of "stable" in InitModelList is overwritten,
cf. GetModelNr(char **name, int *n, int *nr) in RFgetNset.cc
Note that models with overwritten names may be called only by *full* name.
If unhappy, choose any new name instead of overwritting an old one!
See RFgetNset.cc or RFsimu.h for the definition of IncludeModel, addCov,
addTBM, etc.
*/
int nr;
// NEW !!! the following line make the call of PrintModelList()
// unnecessary, so
// library(RandomFields); .C("addownfunctions")
// is sufficient.
if (currentNrCov==-1) InitModelList();
// new stable
nr=IncludeModel("stable", // name that can also be chosen freely
1,
newcheckstable, // NEW
FULLISOTROPIC, false,
infostable, rangestable);
addCov(nr,stable, Dstable, Scalestable);
// addLocal(...);
addTBM(nr, NULL,TBM3stable, CircEmbed, NULL);
/*
int nr, i=1;
char stable[]="stable";
// new stable
GetModelNr(&stable, &i, &nr); // implicit call of InitModelList
ModifyModel(nr, newcheckstable, // NEW
infostable, rangestable);
addCov(nr,stable,Scalestable,
stablelocal, // NEW
stableS // NEW
);
*/
}