Revision 6aa7805ff1eef8ecb881b5099f3cd30f7c2bd637 authored by Martin Schlather on 08 August 1977, 00:00:00 UTC, committed by Gabor Csardi on 08 August 1977, 00:00:00 UTC
1 parent a3c58bc
RFCovFcts.cc
/*
Authors
Martin Schlather, martin.schlather@cu.lu
Definition of correlation functions and derivatives (spectral measures,
tbm operators)
* Never use the below functions directly, but only by the functions indicated
in RFsimu.h, since there is no error check (e.g. initialization of RANDOM)
* when defining your own function, make sure that the covariance function
itself allows for an additional nugget effect (spectral measures and tbm
operators don't)
* VARIANCE, SCALE should not be used
Copyright (C) 2001 -- 2004 Martin Schlather,
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/*
MAKE SURE IN CHECK OF *SPATIALLY* ISOTROPIC (but not fully
isotropic) FUNCTIONS THAT THEY ARE ALWAYS CALLED BY A TIME AND A SPACE
COMPONENT, i.e. x is always a vector of two components !!
dim must always be the true dimension!
*/
#define MINUSINVLOG005 0.3338082006953340674649
#define SQRTINVLOG005 0.5777613700268771079749
#define RANGE_EPSILON 1E-20
#include <math.h>
#include <assert.h>
#include "RFsimu.h"
#include "RFCovFcts.h"
// {min, max} x {theor, pract} x length(param)
static Real range_stable[4] = {0, 2, 0.06, 2};
static Real range_whittle[4]= {0, RF_INF, 1e-2, 10.0};
static Real range_cauchy[4] = {0, RF_INF, 0.09, 10.0};
static Real range_genCauchy[8] = {0, 2, 0.05, 2, 0, RF_INF, 0.05, 10.0};
Real interpolate(Real y, Real *stuetz, int nstuetz, int origin,
Real lambda, int delta)
{
int index,minindex,maxindex,i;
double weights,sum,diff,a;
index = origin + (int) y;
minindex = index - delta; if (minindex<0) minindex=0;
maxindex = index + 1 + delta; if (maxindex>nstuetz) maxindex=nstuetz;
weights = sum = 0.0;
for (i=minindex;i<maxindex;i++) {
diff = y + (Real) (index-i);
a = exp( -lambda * diff * diff);
weights += a;
sum += a * stuetz[i]; // or a/stuetz[i]
}
return (Real) (weights/sum); // and then sum/weights
}
//// NOTE : `*p' may not be changed by any of the functions!
SimulationType methodnugget(int spacedim, bool grid) {return Nugget;}
Real testCov(Real *x,Real *p, int effectivedim){
Real y;
y = fabs(*x);
if (y==0) return 1.0;
return (y * M_E < 1) ? (1.0 + 1.0 / log(y)) : 0;
}
SimulationType methodtest(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return SpectralTBM;
case 2 : return SpectralTBM;
case 3 : return TBM3;
default: return Direct;
}
}
}
/* exponential model */
Real exponential(Real *x,Real *p, int effectivedim){
return exp(-fabs( *x));
}
Real Scaleexponential(Real *p,int scaling){ return MINUSINVLOG005; }
Real TBM2exponential(Real *x, Real *p, int effectivedim)
{
double y;
if (*x==0.0) {return 1.0;}
y = fabs(*x);
return 1.0 - PIHALF * y * I0mL0(y);
}
Real TBM3exponential(Real *x, Real *p, int effectivedim){
register Real y;
y = fabs( *x);
return (1.0-y)*exp(-y);
}
Real Dexponential(Real *x, Real *p, int effectivedim){
return - exp(-fabs( *x));
}
Real spectralexponential(Real *p ) { /* see Yaglom ! */
Real register y;
y = 1.0 - UNIFORM_RANDOM;
return sqrt(1.0 / (y * y) - 1.0);
}
// spectral measure, see Lantuejoul, Geostatistical Simulation (2002)
SimulationType methodexponential(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return SpectralTBM;
case 2 : return SpectralTBM;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangeexponential(int dim, int *index, Real* range){
*index = -1;
}
/* derivative of exponential: qexponential */
Real qexponential(Real *x,Real *p, int effectivedim){
register Real y;
y = exp(-fabs( *x));
return y * (2.0 - p[KAPPA] * y) / (2.0 - p[KAPPA]);
}
Real Scaleqexponential(Real *p,int scaling){
return -1.0 /
log( (1.0 - sqrt(1.0 - p[KAPPA] * (2.0 - p[KAPPA]) * 0.05)) / p[KAPPA]);
}
int checkqexponential(Real *param, int timespacedim, SimulationType method){
if ((param[KAPPA]<0.0) || (param[KAPPA]>1.0)) {
strcpy(ERRORSTRING_OK,"kappa in [0,1]");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
return 0;
}
Real Dqexponential(Real *x,Real *p, int effectivedim){
register Real y;
y = exp(-fabs( *x));
return y * (p[KAPPA] * y - 1.0) * 2.0 / (2.0 - p[KAPPA]);
}
SimulationType methodqexponential(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return Direct;
case 2 : return Direct;
default: return Direct;
}
}
}
void rangeqexponential(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
Real r[4] = {0, 1, 0, 1};
memcpy(range, r, sizeof(Real) * 4);
}
/* derivative of exponential: damped cosine */
Real dampedcosine(Real *x, Real*p, int effectivedim){
register double y;
y = fabs( *x);
return exp(-y * p[KAPPA]) * cos(y);
}
Real Scaledampedcosine(Real *p,int scaling){
if (scaling==NATSCALE_EXACT) return 0.0;
return MINUSINVLOG005;
}
Real TBM3dampedcosine(Real *x, Real *p, int effectivedim){
register double y;
y = fabs( *x);
return exp(-p[KAPPA]*y) * ((1.0 - p[KAPPA] * y) * cos(y) - y * sin(y));
}
Real Ddampedcosine(Real *x, Real *p, int effectivedim){
register double y;
y = fabs( *x);
return - exp(-p[KAPPA]*y) * (p[KAPPA] * cos(y) + sin(y));
}
int checkdampedcosine(Real *param, int timespacedim, SimulationType method){
if ((timespacedim==3) || ((method==TBM3) && (timespacedim<=3))){
if (param[KAPPA]<1.73205080756889) {
strcpy(ERRORSTRING_OK,
"kappa >= sqrt(3) for 3-dimensional simulations and techniques");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
} else if (timespacedim==2){
if (param[KAPPA]<1.0) {
strcpy(ERRORSTRING_OK,"kappa >= 1.0 for 2 dimensions");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
} else if (timespacedim==1){
if (param[KAPPA]<0.0) {
strcpy(ERRORSTRING_OK,"kappa >= 0.0 for 1 dimension");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
} else { // timespacedim ==4
return ERRORNOTPROGRAMMED;
}
return 0;
}
SimulationType methoddampedcosine(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangedampedcosine(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
range[1] = RF_INF;
range[3] = 10.0;
switch(dim) {
case 1 : range[2] = (range[0] = 0.0) + 1e-10; break;
case 2 : range[2] = (range[0] = 1.0) + 1e-10; break;
case 3 : range[2] = (range[0] = RF_M_SQRT_3) + 1e-10; break;
default: *index = -2; range[0] = range[2] = RF_INF;
}
}
/* Class of "spherical" models */
/* circular model */
Real circular(Real *x, Real *p, int effectivedim)
{
Real y;
if ((y=fabs( *x))>1.0) return 0.0;
return ( 1.0 - (2.0 * (y * sqrt(1.0- y * y) + asin(y))) / PI);
}
Real Scalecircular(Real *p,int scaling) {return 1.138509531721630274603;}
// spectral measure, see Lantue !!
SimulationType methodcircular(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return Forbidden;
}
} else {
switch(spacedim) {
case 1 : return Direct;
case 2 : return Direct;
default: return Forbidden;
}
}
}
void rangecircular(int dim, int *index, Real* range){
if (dim<=2) *index=-1; else *index=-2;
}
/* spherical model */
Real spherical(Real *x, Real *p, int effectivedim){
register Real y;
if ((y= fabs( *x))>1.0) {return 0.0;}
return (1.0+y*0.5*(y*y-3.0));
}
Real Scalespherical(Real *p,int scaling){ return 1.23243208931941;}
Real TBM2spherical(Real *x, Real *p, int effectivedim){
register Real y, y2;
y=fabs( *x); y2=y*y;
if (effectivedim <= 2) {
if (y>1.0)
{return (1.0- 0.75 * y * ((2.0 - y2) * asin(1.0/y) +
sqrt(y2 -1.0)));}
return (1.0 - 0.375 * PI * y * (2.0 - y2));
} else {
assert(false);
}
}
Real TBM3spherical(Real *x, Real *p, int effectivedim){
register Real y;
if ((y=fabs( *x))>1.0) {return 0.0;}
return (1.0 + (-3.0 + 2.0 * y * y) * y);
}
Real Dspherical(Real *x, Real *p, int effectivedim){
register Real y;
if ((y=fabs( *x))>1.0) {return 0.0;}
return 1.5 * (y * y - 1.0);
}
// spectral measure, see Lantue !!
SimulationType methodspherical(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
case 3 : return CircEmbed;
default: return Forbidden;
}
} else {
switch(spacedim) {
case 1 : return TBM2;
case 2 : return TBM2;
case 3 : return TBM3;
default: return Forbidden;
}
}
}
void rangespherical(int dim, int *index, Real* range){
if (dim<=3) *index=-1; else *index=-2;
}
/* power model */
Real power(Real *x, Real *p, int effectivedim){
register Real y;
if ((y=fabs( *x))>1.0) {return 0.0;}
return pow(1.0 - y, p[KAPPA]);
}
Real Scalepower(Real *p,int scaling){
return 1.0 / (1.0 - pow(0.05,1/p[KAPPA]));
}
Real TBM2power(Real *x, Real *p, int effectivedim){
// only kappa=2 up to now !
register Real y;
y=fabs( *x);
if (y>1.0)
{return (1.0 - 2.0 * y *(asin(1.0/y) - y + sqrt(y*y-1.0) ));}
return (1.0 - y * (PI - 2.0 * y));
}
Real TBM3power(Real *x, Real *p, int effectivedim){
register Real y;
if ((y=fabs( *x))>1.0) {return 0.0;}
return (1.0 - y - y * p[KAPPA]) *
pow(1.0 - y, p[KAPPA]-1.0);
}
Real Dpower(Real *x, Real *p, int effectivedim){
register Real y;
if ((y=fabs( *x))>1.0) {return 0.0;}
return - p[KAPPA] * pow(1.0 - y, p[KAPPA]-1.0);
}
int checkpower(Real *param, int timespacedim, SimulationType method) {
switch (method) {
case TBM2 :
if (param[KAPPA]!=2) { // not that kappa==2 if fine also for
// 3 dimensions!
strcpy(ERRORSTRING_OK,"kappa=2 in case of TBM2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
break;
case TBM3 :
if (param[KAPPA]<2) {
strcpy(ERRORSTRING_OK,"kappa>=2");
sprintf(ERRORSTRING_WRONG,"method=TBM3 and kappa=%f",
param[KAPPA]);
return ERRORCOVFAILED;
}
if (timespacedim<=3) break; //!!
default:
if (param[KAPPA] < 0.5 * (1.0 + timespacedim)) {
strcpy(ERRORSTRING_OK,"kappa >= (dim+1)/2");
sprintf(ERRORSTRING_WRONG,"kappa=%f dim=%d",param[KAPPA],
timespacedim);
return ERRORCOVFAILED;
}
}
return 0;
}
SimulationType methodpower(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
// range definition:
// 0: min, theory, 1:max, theory
// 2: min, practically 3:max, practically
void rangepower(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
range[0] = range[2] = 0.5 * ((Real) dim + 1);
range[1] = RF_INF;
range[3] = 20.0;
}
/* stable model */
Real stable(Real *x,Real *p, int effectivedim){
if (*x==0.0) {return 1.0;}
return (exp(-pow(fabs( *x),p[KAPPA])));
}
Real Scalestable(Real *p,int scaling){return pow(MINUSINVLOG005,1/p[KAPPA]);}
Real TBM3stable(Real *x, Real *p, int effectivedim){
register double y;
y = pow(fabs( *x),p[KAPPA]);
return exp(-y) * (1 - p[KAPPA] * y);
}
Real Dstable(Real *x, Real *p, int effectivedim){
register double y,z;
if ( (z = fabs( *x)) == 0.0) return 0.0; /* WRONG VALUE, but multiplied
with 0 anyway*/
y = pow(z,p[KAPPA] - 1.0);
return - p[KAPPA] * exp(-y * z);
}
int checkstable(Real *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;
}
return 0;
}
SimulationType methodstable(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangestable(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
memcpy(range, range_stable, sizeof(Real) * 4);
}
/* Whittle-Matern or Whittle or Besset */
Real WhittleMatern(Real *x, Real *p, int effectivedim)
// check calling functions, like hyperbolic and gneiting if any changings !!
{
static double kappa=RF_INF;
static double loggamma;
register double y;
if ( *x==0.0) {return 1.0;}
y = fabs( *x);
if (kappa!=p[KAPPA]) {
kappa=p[KAPPA];
loggamma = lgammafn(kappa);
}
return 2.0 *
exp(kappa*log(0.5*y)-loggamma + log(bessel_k(y,kappa,2.0))-y);
}
Real ScaleWhittleMatern(Real *p,int scaling){
// it is working reasonably well if kappa is in [0.001,100]
// happy to get any better suggestion!!
static int nstuetz = 19;
static Real stuetz[]=
{1.41062516176753e-14, 4.41556861847671e-12, 2.22633601732610e-06,
1.58113643548649e-03, 4.22181082102606e-02, 2.25024764696152e-01,
5.70478218148777e-01, 1.03102016706644e+00, 1.57836638352906e+00,
2.21866372852304e+00, 2.99573229151620e+00, 3.99852231863082e+00,
5.36837527567695e+00, 7.30561120838150e+00, 1.00809957038601e+01,
1.40580075785156e+01, 1.97332533513488e+01, 2.78005149402352e+01,
3.92400265713477e+01};
static int stuetzorigin = 11;
if (scaling==NATSCALE_EXACT) return 0.0;
return interpolate(log(p[KAPPA]) * INVLOG2, stuetz, nstuetz, stuetzorigin,
1.5, 5);
}
Real TBM2WhittleMatern(Real *x, Real *p, int effectivedim)
{
if (p[KAPPA]==0.5) return TBM2exponential(x, p, effectivedim);
assert(false);
}
Real TBM3WhittleMatern(Real *x, Real *p, int effectivedim)
// check calling functions, like hyperbolic and gneiting if any changings !!
{
static double kappa=RF_INF;
static double loggamma;
register double y,loghalfy;
if ( *x==0.0) {return 1.0;}
y = fabs( *x);
loghalfy = log(0.5*y);
if (kappa!=p[KAPPA]) {
kappa=p[KAPPA];
loggamma = lgammafn(kappa);
}
return
(2.0 * exp(kappa * loghalfy -loggamma + log(bessel_k(y, kappa, 2.0)) - y)
-(4.0 * exp((kappa + 1.0) * loghalfy - loggamma +
log(bessel_k(y, kappa - 1.0, 2.0)) - y)));
}
#define LOG05 -0.69314718055994528623
Real DWhittleMatern(Real *x, Real *p, int effectivedim)
// check calling functions, like hyperbolic and gneiting if any changings !!
{
static double kappa=RF_INF;
static double loggamma;
register double y;
if ( *x==0.0) {return 1.0;}
y = fabs( *x);
if (kappa!=p[KAPPA]) {
kappa=p[KAPPA];
loggamma = lgammafn(kappa);
}
return -2.0 *exp(kappa * log(0.5*y) - loggamma +
log(bessel_k(y, kappa-1.0, 2.0)) - y);
}
Real spectralWhittleMatern(Real *p ) { /* see Yaglom ! */
return sqrt(pow(1.0 - UNIFORM_RANDOM,-1.0/p[KAPPA]) - 1.0);
}
int checkWhittleMatern(Real *param, int timespacedim, SimulationType method)
{
static Real spectrallimit=0.17;
switch(method) {
case TBM2 :
if (param[KAPPA]!=0.5) return ERRORNOTPROGRAMMED;
break;
case SpectralTBM :
if ((param[KAPPA]>=spectrallimit)) return 0;
sprintf(ERRORSTRING_OK,
"%f<=kappa (the numerical errors are too big for 0<kappa<%f)",
spectrallimit,spectrallimit);
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
default :
if ((param[KAPPA]>0)) return 0;
strcpy(ERRORSTRING_OK,"0<kappa");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType methodWhittleMatern(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return SpectralTBM;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return SpectralTBM;
case 2 : return SpectralTBM;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangeWhittleMatern(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
memcpy(range, range_whittle, sizeof(Real) * 4);
}
Real hyperbolic(Real *x, Real*p, int effectivedim){
static double kappa = RF_INF;
static double lambda= RF_INF;
static double delta = RF_INF;
static double deltasq;
static double kappadelta;
static double logconst;
double kappay;
Real y;
if ( *x==0.0) {return 1.0;}
if (p[KAPPAIII]==0) { // whittle matern
y = *x * p[KAPPAI];
return WhittleMatern(&y, p, effectivedim);
}
if (p[KAPPAI]==0) { //cauchy => KAPPAII < 0 !
y = *x/p[KAPPAIII];
/* note change in sign as KAPPAII<0 */
return pow(1+y*y,p[KAPPAII]);
}
y=fabs( *x);
if ((p[KAPPAI]!=kappa) || (p[KAPPAII]!=lambda) || (p[KAPPAIII]!=delta)) {
kappa=p[KAPPAI];
lambda= p[KAPPAII];
delta = p[KAPPAIII];
deltasq = delta * delta;
kappadelta = kappa * delta;
logconst = kappadelta - log(bessel_k(kappadelta,lambda,2.0))
- lambda * log(delta);
}
y=sqrt(deltasq + y*y);
kappay = kappa * y;
return
exp(logconst + lambda * log(y)+log(bessel_k(kappay,lambda,2.0))-kappay);
}
Real TBM3hyperbolic(Real *x, Real*p, int effectivedim)
{
static double kappa = RF_INF;
static double lambda= RF_INF;
static double delta = RF_INF;
static double deltasq;
static double kappadelta;
static double logconst;
Real y;
double ysq,s,kappas,logs;
if ( *x==0.0) {return 1.0;}
if (p[KAPPAIII]==0) { // whittle matern
y = *x * p[KAPPAI];
return TBM3WhittleMatern(&y, p, effectivedim);
}
if (p[KAPPAI]==0) { //cauchy
register Real y,ha;
y= *x/p[KAPPAIII];
ha=y*y;
/* note change in sign as KAPPAII<0 */
return (1.0+ (1.0+2.0 * p[KAPPAII]) * ha) *
pow(1.0+ha,p[KAPPAII]-1.0);
}
y=fabs( *x);
if ((p[KAPPAI]!=kappa) || (p[KAPPAII]!=lambda) || (p[KAPPAIII]!=delta)) {
kappa=p[KAPPAI];
lambda= p[KAPPAII];
delta = p[KAPPAIII];
deltasq = delta * delta;
kappadelta = kappa * delta;
logconst = kappadelta - log(bessel_k(kappadelta,lambda,2.0))
- lambda * log(delta);
}
ysq = y * y;
s=sqrt(deltasq + ysq);
kappas = kappa * s;
logs = log(s);
return
( exp(logconst + lambda * logs +log(bessel_k(kappas,lambda,2.0))-kappas)
- ysq*kappa*exp(logconst + (lambda-1.0)*logs
+log(bessel_k(kappas,lambda-1.0,2.0))-kappas)
);
}
Real Dhyperbolic(Real *x, Real*p, int effectivedim)
{
static double kappa = RF_INF;
static double lambda= RF_INF;
static double delta = RF_INF;
static double deltasq;
static double kappadelta;
static double logconst;
Real y;
double s,kappas,logs;
if ( *x==0.0) {return 1.0;}
if (p[KAPPAIII]==0) { // whittle matern
y = *x * p[KAPPAI];
return DWhittleMatern(&y, p, effectivedim);
}
if (p[KAPPAI]==0) { //cauchy
register Real y,ha;
y= *x/p[KAPPAIII];
ha=y*y;
/* note change in sign as KAPPAII<0 */
return 2.0 * p[KAPPAII] * fabs(y) *
pow(1.0+ha,p[KAPPAII]-1.0);
}
y=fabs( *x);
if ((p[KAPPAI]!=kappa) || (p[KAPPAII]!=lambda) || (p[KAPPAIII]!=delta)) {
kappa=p[KAPPAI];
lambda= p[KAPPAII];
delta = p[KAPPAIII];
deltasq = delta * delta;
kappadelta = kappa * delta;
logconst = kappadelta - log(bessel_k(kappadelta,lambda,2.0))
- lambda * log(delta);
}
s=sqrt(deltasq + y * y);
kappas = kappa * s;
logs = log(s);
return
(
- y * kappa*exp(logconst + (lambda-1.0)*logs
+log(bessel_k(kappas,lambda-1.0,2.0))-kappas)
);
}
int checkhyperbolic(Real *param, int timespacedim, SimulationType method){
if (param[KAPPAII]>0) {
if ((param[KAPPAIII]<0) || (param[KAPPAI]<=0)) {
strcpy(ERRORSTRING_OK,"kappa1>0 and kappa3>=0 if kappa2>0");
sprintf(ERRORSTRING_WRONG,"kappa1=%f and kappa3=%f for kappa2=%f",
param[KAPPAI],param[KAPPAIII],param[KAPPAII]);
return ERRORCOVFAILED;
}
} else if (param[KAPPAII]<0) {
if ((param[KAPPAIII]<=0) || (param[KAPPAI]<0)) {
strcpy(ERRORSTRING_OK,"kappa1>=0 and kappa3>0 if kappa2<0");
sprintf(ERRORSTRING_WRONG,"kappa1=%f and kappa3=%f for kappa2=%f",
param[KAPPAI],param[KAPPAIII],param[KAPPAII]);
return ERRORCOVFAILED;
}
} else { // param[KAPPAII]==0.0
if ((param[KAPPAIII]<=0) || (param[KAPPAI]<=0)) {
strcpy(ERRORSTRING_OK,"kappa1>0 and kappa3>0 if kappa2=0");
sprintf(ERRORSTRING_WRONG,"kappa1=%f and kappa3=%f for kappa2=%f",
param[KAPPAI],param[KAPPAIII],param[KAPPAII]);
return ERRORCOVFAILED;
}
}
return 0;
}
SimulationType methodhyperbolic(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangehyperbolic(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
Real r[12] = {0, RF_INF, 0.000001, 10.0,
RF_NEGINF, RF_INF, -20.0, 20.0,
0, RF_INF, 0.000001, 10.0
};
*index = -1;
memcpy(range, r, sizeof(Real) * 12);
}
/* Gneiting's functions */
// alternative to Gaussian */
#define Sqrt2TenD47 0.30089650263257344820 /* approcx 0.3 ?? */
#define NumericalScale 0.301187465825
Real Gneiting(Real *x, Real *p, int effectivedim){
register Real y,oneMy8;
if ((y=fabs(*x * NumericalScale))>1.0) {return(0.0);}
oneMy8 = 1.0-y; oneMy8*=oneMy8; oneMy8*=oneMy8; oneMy8*=oneMy8;
return ((1.0+y * ( 8.0 + y * (25.0 + 32.0 *y)))*oneMy8);
}
Real ScaleGneiting(Real *p,int scaling) {return 0.5854160193;}
Real TBM3Gneiting(Real *x, Real *p, int effectivedim){
register Real y,oneMy7;
if ((y=fabs( *x*Sqrt2TenD47))>1.0) {return 0.0;}
oneMy7 = 1.0-y; oneMy7*=oneMy7; oneMy7 *= oneMy7 * oneMy7 * (1.0-y);
return
(1.0 + y * (7.0 - y * (5.0 + y * (147.0 + 384.0 * y))))* oneMy7;
}
Real DGneiting(Real *x, Real *p, int effectivedim){
register Real y,oneMy7;
if ((y=fabs( *x*Sqrt2TenD47))>1.0) {return 0.0;}
oneMy7 = 1.0-y; oneMy7*=oneMy7; oneMy7 *= oneMy7 * oneMy7 * (1.0-y);
return
(-y) * ( 22.0 + y * (154.0 + y * 352.0)) * oneMy7;
}
SimulationType methodGneiting(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangeGneiting(int dim, int *index, Real* range){
if (dim<=3) *index=-1; else *index=-2;
}
Real genGneiting(Real *x, Real *p, int effectivedim)
{
register Real y, s;
if ((y=fabs( *x))>1.0) return 0.0;
s = p[KAPPAII] + p[KAPPAI];
switch ((int) p[KAPPAI]) {
case 1:
return pow ((1.0-y), s) * (1.0 + s*y);
case 2:
return pow ((1.0-y), s) *
(1.0 + y * (s + y * (s*s-1.0)*0.3333333333333333));
case 3:
return pow ((1.0-y), s) *
(1.0 + y * (s
+ y * ( 0.2 * (2.0*s*s - 3.0)
+ y * (s*s-4.0)*s*0.0666666666666666666))) ;
default : assert(false);
}
}
Real TBM3genGneiting(Real *x, Real *p, int effectivedim)
{
register Real y, s;
if ((y=fabs( *x))>1.0) {return 0.0;}
s = p[KAPPAII] + p[KAPPAI];
switch ((int) p[KAPPAI]) {
case 1:
return pow((1.0 - y), p[KAPPAII]) *
(1.0 + y * (p[KAPPAII] - y * s * (p[KAPPAII] + 3.0)));
case 2:
return pow((1.0 - y), (p[KAPPAII]+1.0)) *
(1.0 + y * ((p[KAPPAII]+1.0)
- y *(1.0 + 2.0 * s
+ y * (s * s - 1.0) * 0.33333333333333 *
(p[KAPPAII] + 5.0))));
case 3:
return pow((1.0 - y), (p[KAPPAII]+2.0)) *
(1.0 + y * ((p[KAPPAII]+2.0)
+ y * (0.2 * (-9.0 + s * (-10.0 + s))
+ y * 0.0666666666666666 *
((27.0 - s*(7.0 + s*(18.0 + 2.0*s)))
- y * (s*s - 4.0) * s * (p[KAPPAII] + 7.0)))));
default : assert(false);
}
}
Real DgenGneiting(Real *x, Real *p, int effectivedim)
{
register Real y, s;
if ((y=fabs( *x))>1.0) {return 0.0;}
s = p[KAPPAII] + p[KAPPA];
switch ((int) p[KAPPAI]) {
case 1:
return - pow(1.0 - y, p[KAPPAII]) *
y * (p[KAPPAII] + 1.0) * (p[KAPPAII] + 2.0) ;
case 2:
return - pow(1.0 - y, s - 1.0) *
y * (0.333333333333333333 * s * s + 0.66666666666666666667 + s +
y * 0.333333333333333333 * (s * s - 1.0) * (p[KAPPAII] + 4.0) );
case 3:
return - pow(1.0 - y, s - 1.0) *
0.2 * y * (6 + s * (5.0 + s) +
y * (-6 + s * (1.0 + s * (4.0 + s)) +
- 0.33333333333333333 * y * s * (-12.0 + s * (-4.0 + s * (3.0 + s)))
));
default : assert(false);
}
}
int checkgenGneiting(Real *param, int timespacedim, SimulationType method)
{
if ((param[KAPPAI] != (double) ((int)param[KAPPAI])) ||
(param[KAPPAI]<0)) {
strcpy(ERRORSTRING_OK,"positive integer kappa");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAI]);
return ERRORCOVFAILED;
}
if ((param[KAPPAI] < 1.0) || (param[KAPPAI] > 3.0))
return ERRORNOTPROGRAMMED;
if ((method==TBM3) &&
(param[KAPPAII] < 2.0 + param[KAPPAI])) {
strcpy(ERRORSTRING_OK,"kappa2 <= 2.0 + kappa1");
sprintf(ERRORSTRING_WRONG, "method=TBM3, kappa1=%f and kappa2=%f ",
param[KAPPAI], param[KAPPAII]);
return ERRORCOVFAILED;
}
if (param[KAPPAII] < ( 0.5 * timespacedim +
param[KAPPAI] + 0.5)) {
strcpy(ERRORSTRING_OK, "kappa2 >= (dim + 2*kappa1 + 1)/2 ");
sprintf(ERRORSTRING_WRONG, "kappa1=%f, kappa2=%f, and dim=%d ",
param[KAPPAI], param[KAPPAII], timespacedim);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType methodgenGneiting(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangegenGneiting(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
range[0] = range[2] = range[1] = range[3] = *index;
range[4] = range[6] = 0.5 * (Real) (dim + 2 * *index + 1);
range[5] = RF_INF;
range[7] = 20.0;
if ((++(*index)) > 3) *index=-1;
if (dim>3) *index=-2;
}
/* Gausian model */
Real Gauss(Real *x, Real*p, int effectivedim) {
return exp(- *x * *x);
}
Real ScaleGauss(Real *p,int scaling) {return SQRTINVLOG005;}
Real TBM3Gauss(Real *x, Real*p, int effectivedim) {
register Real y;
y = *x * *x;
return (1-2.0*y)*exp(-y);
}
Real DGauss(Real *x, Real*p, int effectivedim) {
register Real y;
y = fabs( *x);
return - 2.0 * y * exp(- y * y);
}
Real spectralGauss(Real *p ) {
return 2.0 * sqrt(-log(1.0 - UNIFORM_RANDOM));
}
SimulationType methodGauss(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return SpectralTBM;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return SpectralTBM;
case 2 : return SpectralTBM;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangeGauss(int dim, int *index, Real* range){
*index=-1;
}
/* Cauchy models */
Real Cauchy(Real *x, Real *p, int effectivedim){
return pow(1+ *x * *x,-p[KAPPA]);
}
Real ScaleCauchy(Real *p,int scaling) {
switch(scaling) {
case NATSCALE_EXACT: case NATSCALE_APPROX:
return 1.0/sqrt(pow(0.05,-1/p[KAPPA])-1.0);
break;
case NATSCALE_MLE:
/// this values should be changed, because of
/// possibly long tails!
return 1.0/sqrt(pow(0.05,-1/p[KAPPA])-1.0);
break;
default: assert(false);
}
}
Real TBM2Cauchy(Real *x, Real *p, int effectivedim){
register Real y2, lpy2;
y2 = *x * *x;
lpy2=1.0+y2;
switch ((int) (p[KAPPA]*2.0 + 0.1)) {
case 1 : return 1.0 / lpy2;
case 3 : return (1.0 - y2)/ (lpy2*lpy2);
case 5 : return
(1.0-y2*(2.0+0.333333333333333333333*y2))/(lpy2*lpy2*lpy2);
case 7 : lpy2 *= lpy2;
return (1.0- y2*(3.0+y2*(1.0+0.2*y2)))/(lpy2 * lpy2);
default : assert(false);
}
}
Real TBM3Cauchy(Real *x, Real *p, int effectivedim){
register Real ha;
ha= *x * *x;
return (1.0+ (1.0-2.0*p[KAPPA])*ha) * pow(1.0+ha,-p[KAPPA]-1.0);
}
Real DCauchy(Real *x, Real *p, int effectivedim){
register Real y;
y = fabs( *x);
return (-2.0 * p[KAPPA] * y) * pow(1.0 + y * y, -p[KAPPA]-1.0);
}
int checkCauchy(Real *param, int timespacedim, SimulationType method){
switch (method) {
case TBM2 :
if (timespacedim>2) return ERRORCOVFAILED;
if ((param[KAPPA]!=0.5) && (param[KAPPA]!=1.5) &&
(param[KAPPA]!=2.5) && (param[KAPPA]!=3.5)) {
strcpy(ERRORSTRING_OK,"kappa in {0.5, 1.5, 2.5 ,3.5} (TBM2)");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
default :
if (param[KAPPA]<=0) {
strcpy(ERRORSTRING_OK,"0<kappa");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
}
return 0;
}
SimulationType methodCauchy(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return TBM3;
default: return TBM3;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangeCauchy(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
memcpy(range, range_cauchy, sizeof(Real) * 4);
}
Real generalisedCauchy(Real *x, Real *p, int effectivedim){
return pow(1.0 + pow(fabs(*x), p[KAPPAI]), -p[KAPPAII]/p[KAPPAI]);
}
Real ScalegeneralisedCauchy(Real *p,int scaling) {
switch(scaling) {
case NATSCALE_EXACT: case NATSCALE_APPROX:
return pow(pow(0.05,-p[KAPPAI]/p[KAPPAII])-1.0,-1.0/p[KAPPAI]);
break;
case NATSCALE_MLE:
// should be changed! (long tails!)
return pow(pow(0.05,-p[KAPPAI]/p[KAPPAII])-1.0,-1.0/p[KAPPAI]);
break;
default: assert(false);
}
}
Real TBM3generalisedCauchy(Real *x, Real *p, int effectivedim){
register Real ha;
ha=pow(fabs( *x),p[KAPPAI]);
return
(1.0+ (1.0-p[KAPPAII])*ha) * pow(1.0+ha,-p[KAPPAII]/p[KAPPAI]-1.0);
}
Real DgeneralisedCauchy(Real *x, Real *p, int effectivedim){
register Real ha,y;
if ((y = fabs(*x))==0.0) return 0.0; // WRONG VALUE, but multiplied with
// zero anyway
ha=pow(y,p[KAPPAI] - 1.0);
return - p[KAPPAII] * ha * pow(1.0 + ha * y,-p[KAPPAII]/p[KAPPAI]-1.0);
}
int checkgeneralisedCauchy(Real *param, int timespacedim, SimulationType method){
if ((param[KAPPAI]<=0) || (param[KAPPAI]>2.0)) {
strcpy(ERRORSTRING_OK,"0<kappa1<=2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAI]);
return ERRORCOVFAILED;
if (param[KAPPAII]>0) return 0;
}
if (param[KAPPAII]<=0) {
strcpy(ERRORSTRING_OK,"0<kappa2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAII]);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType methodgeneralisedCauchy(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return TBM3;
default: return TBM3;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangegeneralisedCauchy(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
memcpy(range, range_genCauchy, sizeof(Real) * 8);
}
Real Cauchytbm(Real *x, Real *p, int effectivedim){
register Real ha;
if ( *x==0) {return 1.0;}
ha=pow(fabs( *x),p[KAPPAI]);
return
(1.0+ (1.0-p[KAPPAII]/p[KAPPAIII])*ha) * pow(1.0+ha,-p[KAPPAII]/p[KAPPAI]-1.0);
}
Real TBM3Cauchytbm(Real *x, Real *p, int effectivedim){
register Real bg,ha;
ha=pow(fabs( *x),p[KAPPAI]);
bg=p[KAPPAII]/p[KAPPAIII];
return
(1 + ha * (1-bg*(1+p[KAPPAI])+(1-p[KAPPAII]) * (1+(1-bg)*ha)))*
pow(1+ha,-p[KAPPAII]/p[KAPPAI]-2.0);
}
Real DCauchytbm(Real *x, Real *p, int effectivedim){
register Real y,ha;
if ((y = fabs(*x)) == 0.0) return 0.0; // WRONG VALUE, but multiplied
// zero anyway
ha=pow(y,p[KAPPAI] - 1.0);
return
p[KAPPAII] * ha * (-1.0 - p[KAPPAI]/p[KAPPAIII] +
ha * y * (p[KAPPAII]/p[KAPPAIII] - 1.0)) *
pow(1.0 + ha * y,-p[KAPPAII]/p[KAPPAI]-2.0);
}
int checkCauchytbm(Real *param, int timespacedim, SimulationType method){
if ((method==TBM3) && (3.0 > param[KAPPAIII])) {
strcpy(ERRORSTRING_OK,"3 <= kappa3 and method=TBM3");
sprintf(ERRORSTRING_WRONG,
"method=TBM3 and kappa3=%f ",param[KAPPAIII]);
return ERRORCOVFAILED;
}
if (timespacedim > param[KAPPAIII]) {
strcpy(ERRORSTRING_OK,"dim <= kappa3");
sprintf(ERRORSTRING_WRONG,
"kappa3=%f < dim=%d ",param[KAPPAIII], timespacedim);
return ERRORCOVFAILED;
}
// theory: should check whether p[KAPPAIII] is an integer.
// practically, it goes through without this restriction !!
if ((param[KAPPAI]<=0) || (param[KAPPAI]>2.0)) {
strcpy(ERRORSTRING_OK,"0<kappa1<=2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAI]);
return ERRORCOVFAILED;
}
if (param[KAPPAII]<=0) {
strcpy(ERRORSTRING_OK,"0<kappa2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAII]);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType methodCauchytbm(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return TBM3;
default: return TBM3;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangeCauchytbm(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
memcpy(range, range_genCauchy, sizeof(Real) * 8);
range[8] = range[10] = (double) dim;
range[9] = RF_INF;
range[11] = (double) dim + 10.0;
}
/* Bessel function */
Real Bessel(Real *x,Real *p, int effectivedim){
static double kappa=RF_INF;
static double gamma;
Real y;
if ( *x==0.0) {return 1.0;}
y = fabs( *x);
if (kappa!=p[KAPPA]) {
kappa=p[KAPPA];
gamma = gammafn(kappa+1.0);
}
return gamma * pow(0.5 * y,-kappa) * bessel_j(y,kappa);
}
Real spectralBessel(Real *p ) {
return p[KAPPA]==0 ? 1.0 :
(sqrt(1-pow(UNIFORM_RANDOM,1.0/p[KAPPA])));
}
int checkBessel(Real *param, int timespacedim, SimulationType method){
// Whenever TBM3Bessel exists, add further check against too small kappa!
if (param[KAPPA]<0.5*((double) timespacedim) - 1.0) {
strcpy(ERRORSTRING_OK,"kappa >= dim/2 - 1");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
if (method==SpectralTBM && param[KAPPA]<0.0) {
strcpy(ERRORSTRING_OK,"spectral TBM and kappa>=0");
sprintf(ERRORSTRING_WRONG,"%f", param[KAPPA]);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType methodBessel(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return SpectralTBM;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return SpectralTBM;
case 2 : return SpectralTBM;
case 3 : return Direct;
default: return Direct;
}
}
}
void rangeBessel(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
range[2] = 0.0001 + (range[0] = 0.5 * ((Real) dim - 2.0));
range[1] = RF_INF;
range[3] = 10.0;
}
Real wave(Real *x, Real *p, int effectivedim) {
if ( *x==0.0) { return 1.0;}
return sin( *x)/ *x;
}
Real Scalewave(Real *p,int scaling) {return 0.302320850755833;}
Real spectralwave(Real *p ) {
Real x; x=UNIFORM_RANDOM; return sqrt(1- x * x);
}
SimulationType methodwave(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return SpectralTBM;
case 2 : return SpectralTBM;
case 3 : return Direct;
default: return Direct;
}
}
}
void rangewave(int dim, int *index, Real* range){
if(dim<=3) *index=-1; else *index=-2;
}
Real cubic(Real *x, Real *p, int effectivedim)
{ ///
Real y, y2;
if ((y=fabs(*x)) >= 1.0) {return 0.0;}
y2 = y * y;
return (1.0 + (((0.75 * y2 - 3.5) * y2 + 8.75) * y - 7) * y2);
}
Real TBM3cubic(Real *x, Real *p, int effectivedim)
{ ///
Real y, y2;
if ((y=fabs(*x)) >= 1.0) {return 0.0;}
y2 = y * y;
return (1.0 + y2 * (-21.0 + y * (35.0 + y2 * (-21.0 + 6.0 * y2))));
}
Real Dcubic(Real *x, Real *p, int effectivedim)
{ ///
Real y,y2;
if ((y=fabs(*x)) >= 1.0) {return 0.0;}
y2 = y * y;
return y * (-14.0 + y * (26.25 + y2 * (-17.5 + 5.25 * y2)));
}
Real Scalecubic(Real *p,int scaling) {return 1.44855683156829;}
SimulationType methodcubic(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangecubic(int dim, int *index, Real* range){
if(dim<=3) *index=-1; else *index=-2;
}
Real penta(Real *x, Real *p, int effectivedim)
{ ///
Real y,y2;
if ((y=fabs( *x))>=1.0) { return 0.0; }
y2=y * y;
return
(1.0 + y2 * (-7.333333333333333
+ y2 * (33.0 +
y * (-38.5 +
y2 * (16.5 +
y2 * (-5.5 +
y2 * 0.833333333333333))))));
}
Real TBM3penta(Real *x, Real *p, int effectivedim)
{ ///
Real y,y2;
if ((y=fabs( *x))>=1.0) { return 0.0; }
y2 = y * y;
return
(1.0 + y2 * (-22.0
+ y2 * (165.0
+ y * (-231.0
+ y2 * (132.0
+ y2 * (-55.0
+ 10.0 * y2))))));
}
Real Dpenta(Real *x, Real *p, int effectivedim)
{ ///
Real y,y2;
if ((y=fabs( *x))>=1.0) { return 0.0; }
y2 = y * y;
return
y * (-14.66666666666666667 +
y2 * (132.0 +
y * (-192.5 +
y2 * (115.5 +
y2 * (-49.5 +
y2 * 9.16666666666666667)))));
}
Real Scalepenta(Real *p,int scaling) {return 1.6552838957365;}
SimulationType methodpenta(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangepenta(int dim, int *index, Real* range){
if(dim<=3) *index=-1; else *index=-2;
}
/* Tilmann Gneiting's space time models */
/* following function is unused! -- it has been nsst3*/
Real spacetime3(Real *x,Real *p, int effectivedim){
// turning bands modified spacetime1 covariance function
Real y, z, invsqrtpsi;
int i;
assert(effectivedim==2);
return spacetime1(x,p,effectivedim) + fabs(x[0]) * Dspacetime1(x,p,effectivedim);
}
Real InvSqrtPsi(Real x, Real a, Real b, int c) {
Real y;
y = pow(fabs(x), a);
switch(c) {
case 1 :
return pow(1.0 + y, - 0.5 * b);
case 2 :
return sqrt( (y + 1.0) / (y/b + 1.0) );
case 3 :
return sqrt(-log(b) / log(y + 1.0 / b));
default: assert(false);
}
}
/* Tilmann Gneiting's space time models, part I */
Real spacetime1(Real *x, Real *p, int effectivedim){
Real y, z, invsqrtpsi;
int i;
assert(effectivedim==2);
invsqrtpsi = InvSqrtPsi(x[1], p[KAPPAIII], p[KAPPAIV], (int)p[KAPPAV]);
y = x[0] * invsqrtpsi;
switch((int) p[KAPPAII]) {
case 1 : z = stable(&y, p, 1); break;
case 2 : z = WhittleMatern(&y, p, 1); break;
case 3 : z = Cauchy(&y, p, 1); break;
default: assert(false);
}
return pow(invsqrtpsi, p[KAPPAVI]) * z;
}
Real TBM2spacetime1(Real *x, Real *p, int effectivedim){
Real y, z, invsqrtpsi;
int i;
assert(effectivedim==2);
invsqrtpsi = InvSqrtPsi(x[1], p[KAPPAIII], p[KAPPAIV], (int)p[KAPPAV]);
y = x[0] * invsqrtpsi;
switch((int) p[KAPPAII]) {
case 1 : assert(false); break;
case 2 : z = TBM2WhittleMatern(&y, p, 1); break;
case 3 : z = TBM2Cauchy(&y, p, 1); break;
default: assert(false);
}
return pow(invsqrtpsi, p[KAPPAVI]) * z;
}
Real TBM3spacetime1(Real *x, Real *p, int effectivedim){
assert(false); // should never happen
return spacetime1(x, p, effectivedim) +
Dspacetime1(x, p, effectivedim) * fabs(x[0]);
}
Real Dspacetime1(Real *x, Real *p, int effectivedim){
Real y, z, invsqrtpsi;
int i;
assert(effectivedim==2);
invsqrtpsi = InvSqrtPsi(x[1], p[KAPPAIII], p[KAPPAIV], (int)p[KAPPAV]);
y = x[0] * invsqrtpsi;
switch((int) p[KAPPAII]) {
case 1 : z = Dstable(&y, p, 1); break;
case 2 : z = DWhittleMatern(&y, p, 1); break;
case 3 : z = DCauchy(&y, p, 1); break;
default: assert(false);
}
return pow(invsqrtpsi, p[KAPPAVI] + 1.0) * z;
}
int checkspacetime1(Real *param, int timespacedim, SimulationType method) {
int error;
// 1 : stable
// 2 : whittle
// 3 : cauchy (not generalised)
// first parameter(s) for phi
// then choice of phi; then two parameters for psi, then choice of psi
if (timespacedim<=1) {
strcpy(ERRORSTRING_OK,"total dim>=2");
sprintf(ERRORSTRING_WRONG,"%d", timespacedim);
return ERRORCOVFAILED;
}
if (param[KAPPAII] != (double)((int) param[KAPPAII])) {
strcpy(ERRORSTRING_OK,"kappa3 an integer");
sprintf(ERRORSTRING_WRONG,"%f", param[KAPPAII]);
return ERRORCOVFAILED;
}
switch((int) param[KAPPAII]) {
case 1 :
if (method==TBM2) {
strcpy(ERRORSTRING_OK,"TBM2 and kappa2=2,3");
sprintf(ERRORSTRING_WRONG,"%f", param[KAPPAII]);
return ERRORCOVFAILED;
}
error=checkstable(param, timespacedim, method);
break;
case 2 : error=checkWhittleMatern(param, timespacedim, method); break;
case 3 : error=checkCauchy(param, timespacedim, method); break;
default :
strcpy(ERRORSTRING_OK,"kappa2=1,2,3");
sprintf(ERRORSTRING_WRONG,"%d",(int) param[KAPPAIII]);
return ERRORCOVFAILED;
}
if (error) return error;
if ((param[KAPPAIII]<=0) || (param[KAPPAIII]>2)) {
strcpy(ERRORSTRING_OK,"0<kappa3<=2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAIII]);
return ERRORCOVFAILED;
}
if ((param[KAPPAIV]<=0) || (param[KAPPAIV]>1)) {
strcpy(ERRORSTRING_OK,"0<kappa4<=1");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAIV]);
return ERRORCOVFAILED;
}
if (param[KAPPAV] != (double)((int) param[KAPPAV])) {
strcpy(ERRORSTRING_OK,"kappa5 an integer");
sprintf(ERRORSTRING_WRONG,"%f", param[KAPPAV]);
return ERRORCOVFAILED;
}
switch((int) param[KAPPAV]) {
case 1 : case 2: case 3: break;
default :
strcpy(ERRORSTRING_OK,"kappa5=1,2,3");
sprintf(ERRORSTRING_WRONG,"%d",(int) param[KAPPAV]);
return ERRORCOVFAILED;
}
if (timespacedim-1 > param[KAPPAVI]) {
strcpy(ERRORSTRING_OK,"kappa6>=dim-1");
sprintf(ERRORSTRING_WRONG,"%f for spatial dim=%d",
param[KAPPAVI],timespacedim-1);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType methodspacetime1(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangespacetime1(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
int b;
Real *r;
if ((*index<=0) || (*index>9)) { // see also last line
int i; for (i=0; i<24; i++) range[i]= RF_NAN; *index=-1; return;
}
switch ((*index-1) % 3) {
case 0 : r=range_stable; break;
case 1 : r=range_whittle; break;
case 2 : r=range_cauchy; break;
default: assert(false);
}
memcpy(range, r, sizeof(Real) * 4);
range[4] = range[6] =
range[5] = range[7] = (Real) (1 + (*index - 1) % 3);
range[8] = 0;
range[9] = range[11] = 2;
range[10] = RANGE_EPSILON;
range[12] = 0;
range[13] = 1;
range[14] = RANGE_EPSILON;
range[15] = 1.0 - RANGE_EPSILON;
range[16] = range[17] =
range[18] = range[19] = (Real) (1 + (*index - 1) / 3);
range[20] = range[22] = (Real) dim;
range[21] = RF_INF;
range[23] = (Real) dim + 10.0;
if ( (++(*index)) > 9) *index=-1;
}
/* Tilmann Gneiting's space time models, part II*/
Real spacetime2(Real *x,Real *p, int effectivedim){
Real y, z, invsqrtpsi;
int i;
assert(effectivedim==2);
invsqrtpsi = InvSqrtPsi(x[1], p[KAPPAIV], p[KAPPAV], (int)p[KAPPAVI]);
y = x[0] * invsqrtpsi;
switch((int) p[KAPPAIII]) {
case 1 : z = generalisedCauchy(&y, p, 1); break;
default: assert(false);
}
return pow(invsqrtpsi, p[KAPPAVII]) * z;
}
Real TBM3spacetime2(Real *x, Real *p, int effectivedim){
assert(false); // should never happen
return spacetime2(x, p, effectivedim) +
Dspacetime2(x, p, effectivedim) * fabs(x[0]);
}
Real Dspacetime2(Real *x, Real *p, int effectivedim){
Real y, z, invsqrtpsi;
int i;
assert(effectivedim==2);
invsqrtpsi = InvSqrtPsi(x[1], p[KAPPAIV], p[KAPPAV], (int)p[KAPPAVI]);
y = x[0] * invsqrtpsi;
switch((int) p[KAPPAIII]) {
case 1 : z = DgeneralisedCauchy(&y, p, 1); break;
default: assert(false);
}
return pow(invsqrtpsi, p[KAPPAV] + 1.0) * z;
}
int checkspacetime2(Real *param, int timespacedim, SimulationType method) {
int error;
// 1 : generalisedcauchy
// first parameter(s) for phi
// then choice of phi; then two parameters for psi, then choice of psi
if (timespacedim<=1) {
strcpy(ERRORSTRING_OK,"total dim>=2");
sprintf(ERRORSTRING_WRONG,"%d", timespacedim);
return ERRORCOVFAILED;
}
if (param[KAPPAIII] != (double)((int) param[KAPPAIII])) {
strcpy(ERRORSTRING_OK,"kappa3 an integer");
sprintf(ERRORSTRING_WRONG,"%f", param[KAPPAIII]);
return ERRORCOVFAILED;
}
switch((int) param[KAPPAIII]) {
case 1 : error=checkgeneralisedCauchy(param, timespacedim, method); break;
default :
strcpy(ERRORSTRING_OK,"kappa3=1,2,3");
sprintf(ERRORSTRING_WRONG,"%d",(int) param[KAPPAIII]);
return ERRORCOVFAILED;
}
if (error) return error;
if ((param[KAPPAIV]<=0) || (param[KAPPAIV]>2)) {
strcpy(ERRORSTRING_OK,"0<kappa4<=2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAIV]);
return ERRORCOVFAILED;
}
if ((param[KAPPAV]<=0) || (param[KAPPAV]>1)) {
strcpy(ERRORSTRING_OK,"0<kappa5<=1");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAV]);
return ERRORCOVFAILED;
}
if (param[KAPPAVI] != (double)((int) param[KAPPAVI])) {
strcpy(ERRORSTRING_OK,"kappa6 an integer");
sprintf(ERRORSTRING_WRONG,"%f", param[KAPPAVI]);
return ERRORCOVFAILED;
}
switch((int) param[KAPPAVI]) {
case 1 : case 2: case 3: break;
default :
strcpy(ERRORSTRING_OK,"kappa6=1,2,3");
sprintf(ERRORSTRING_WRONG,"%d",(int) param[KAPPAVI]);
return ERRORCOVFAILED;
}
if (timespacedim>param[KAPPAVII]) {
strcpy(ERRORSTRING_OK,"kappa7>=dim");
sprintf(ERRORSTRING_WRONG,"%d",(int) param[KAPPAVII]);
}
return 0;
}
SimulationType methodspacetime2(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 : return TBM3;
case 2 : return TBM3;
case 3 : return TBM3;
default: return Direct;
}
}
}
void rangespacetime2(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
int b;
Real *r;
if ((*index<=0) || (*index>3)) { // see also last line
int i; for (i=0; i<28; i++) range[i]= RF_NAN; *index=-1; return;
}
switch ((*index-1) % 1) {
case 0 : r = range_genCauchy; break;
default: assert(false);
}
memcpy(range, r, sizeof(Real) * 8);
range[8] = range[10] =
range[9] = range[11] = (Real) (1 + (*index - 1) % 1);
range[12] = 0;
range[13] = range[15] = 2;
range[14] = RANGE_EPSILON;
range[16] = 0;
range[17] = 1;
range[18] = RANGE_EPSILON;
range[19] = 1.0 - RANGE_EPSILON;
range[20] = range[21] =
range[22] = range[23] = (Real) (1 + (*index - 1) / 1);
range[24] = range[26] = (Real) dim;
range[25] = RF_INF;
range[27] = (Real) dim + 10.0;
if ( (++(*index)) > 3) *index= -1;
}
/* locally defined functions */
// Brownian motion
Real fractalBrownian(Real *x, Real *p, int effectivdim) {
return pow(*x,p[KAPPA]);
}
Real twodimfractalBrownianlocal(Real *x, Real *p, int effectivdim) {
register Real y;
if (p[KAPPA]<=1.5) {
if ((y=fabs(*x))>1.0) return 0.0;
return 1.0-pow(y,p[KAPPA])+0.5*p[KAPPA]*(y*y-1.0);
} else {
if ((y=fabs(*x))>2.0) return 0.0;
if (y>1.0) {
register Real z;
z = 2.0-y;
return p[KAPPA] * (2.0-p[KAPPA]) / 18.0 * z * z * z / y;
}
return
1.0 - 0.16666666666666666666666666667 * p[KAPPA] * (1+p[KAPPA])
- pow(y,p[KAPPA]) + p[KAPPA] * (5.0 + 2.0 * p[KAPPA])/ 18.0 * y *y;
}
}
// also corrects the variance and other parameters in param !
Real twodimfractalBrownianS(Real *p) {
if (p[KAPPA]<=1.5) return 0.5*p[KAPPA];
else return p[KAPPA]*(5+2.0*p[KAPPA])/18.0;
}
int check2dfractalBrownian(Real *param, int timespacedim, SimulationType method){
if ((timespacedim>2) && (method!=Nothing)) {
strcpy(ERRORSTRING_OK,"total dim<=2");
sprintf(ERRORSTRING_WRONG,"%d",timespacedim);
return ERRORCOVFAILED;
}
if ((param[KAPPA]<=0.0) || (param[KAPPA]>=2.0)) {
strcpy(ERRORSTRING_OK,"kappa in (0,2)");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
if ((method!=Nothing) && (method!=CircEmbedLocal)) {
strcpy(ERRORSTRING_OK,"method=local CE");
sprintf(ERRORSTRING_WRONG,"%s",METHODNAMES[method]);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType method2dfractalBrownian(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbedLocal;
case 2 : return CircEmbedLocal;
default: return Forbidden;
}
} else {
switch(spacedim) {
case 1 : return Direct;
case 2 : return Direct;
default: return Forbidden;
}
}
}
void range2dfractalBrownian(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
if (dim<=2) *index = -1; else *index=-2;
Real r[4] = {0, 2, RANGE_EPSILON, 2};
memcpy(range, r, sizeof(Real) * 4);
}
Real threedimfractalBrownianlocal(Real *x, Real *p, int effectivedim) {
register Real y;
if (p[KAPPA] <= 1.0) {
if ((y=fabs(*x)) > 1.0) return 0.0;
return 1.0 - pow(y, p[KAPPA]) + 0.5 * p[KAPPA] * (y * y - 1.0);
} else {
if ((y=fabs(*x)) > 2.0) return 0.0;
if (y > 1.0) {
register Real z;
z = 2.0 - y;
return p[KAPPA] * (2.0 - p[KAPPA]) / 18.0 * z * z * z / y;
}
return
1.0 - 0.16666666666666666666666666667 * p[KAPPA] * (1+p[KAPPA])
- pow(y,p[KAPPA]) + p[KAPPA] * (5.0 + 2.0 * p[KAPPA])/ 18.0 * y *y;
}
}
Real threedimfractalBrownianS(Real *p) {
if (p[KAPPA] <= 1.0) return 0.5 * p[KAPPA];
else return p[KAPPA] * (5 + 2.0 * p[KAPPA]) / 18.0;
}
int check3dfractalBrownian(Real *param, int timespacedim, SimulationType method){
if ((timespacedim>3) && (method!=Nothing)) {
strcpy(ERRORSTRING_OK,"total dim<=3");
sprintf(ERRORSTRING_WRONG,"%d",timespacedim);
return ERRORCOVFAILED;
}
if ((param[KAPPA]<=0.0) || (param[KAPPA]>=2.0)) {
strcpy(ERRORSTRING_OK,"kappa in (0,2)");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
if ((method!=Nothing) && (method!=CircEmbedLocal)) {
strcpy(ERRORSTRING_OK,"method=local CE");
sprintf(ERRORSTRING_WRONG,"%s",METHODNAMES[method]);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType method3dfractalBrownian(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbedLocal;
case 2 : return CircEmbedLocal;
case 3 : return CircEmbedLocal;
default: return Forbidden;
}
} else {
switch(spacedim) {
case 1 : return Direct;
case 2 : return Direct;
case 3 : return Direct;
default: return Forbidden;
}
}
}
void range3dfractalBrownian(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
if (dim<=3) *index = -1; else *index=-2;
Real r[4] = {0, 2, RANGE_EPSILON, 2};
memcpy(range, r, sizeof(Real) * 4);
}
Real fractGauss(Real *x, Real *p, int effectivedim){
register Real y;
if ((y = fabs(*x))==0.0) return 1.0;
return 0.5 * (pow(fabs(y + 1.0), p[KAPPA])
- 2.0 * pow(y, p[KAPPA])
+ pow(fabs(y - 1.0), p[KAPPA])
);
}
int checkfractGauss(Real *param, int timespacedim, SimulationType method) {
if ((timespacedim!=1) && (method!=Nothing)) {
strcpy(ERRORSTRING_OK,"dim=1");
sprintf(ERRORSTRING_WRONG,"%d",timespacedim);
return ERRORCOVFAILED;
}
if ((param[KAPPA]<=0) || (param[KAPPA]>2.0)) {
strcpy(ERRORSTRING_OK,"0<kappa<=2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPA]);
return ERRORCOVFAILED;
}
return 0;
}
SimulationType methodfractGauss(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
default: return Forbidden;
}
} else {
switch(spacedim) {
case 1 : return Direct;
default: return Forbidden;
}
}
}
void rangefractGauss(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
if (dim==1) *index = -1; else *index=-2;
Real r[4] = {0, 2, RANGE_EPSILON, 2};
memcpy(range, r, sizeof(Real) * 4);
}
/* local-global distinguisher */
Real lgd1(Real *x, Real*p, int effectivedim) {
Real y;
if ((y=fabs(*x)) < 1) return 1.0 - p[KAPPAII] / (p[KAPPAI] + p[KAPPAII])
* exp(p[KAPPAI] * log(y));
else return p[KAPPAI] / (p[KAPPAI] + p[KAPPAII])
* exp( -p[KAPPAII] * log(y));
}
Real Scalelgd1(Real *p,int scaling) {
if (19 * p[KAPPAI] < p[KAPPAII])
return exp( log(0.95 * (p[KAPPAI] + p[KAPPAII]) / p[KAPPAII]) / p[KAPPAI]);
else return exp(log(0.05 * (p[KAPPAI] + p[KAPPAII]) / p[KAPPAI])/p[KAPPAII]);
}
/*
Real TBM3(Real *x, Real*p, int effectivedim) {
register Real y;
}
Real D(Real *x, Real*p, int effectivedim) {
register Real y;
}
*/
SimulationType methodlgd1(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return CircEmbed;
case 2 : return CircEmbed;
default: return CircEmbed;
}
} else {
switch(spacedim) {
case 1 :
case 2 :
case 3 :
default: return Direct;
}
}
}
int checklgd1(Real *param, int timespacedim, SimulationType method) {
if ((timespacedim>2) && (method!=Nothing)) {
strcpy(ERRORSTRING_OK,"dim<=2");
sprintf(ERRORSTRING_WRONG,"%d",timespacedim);
return ERRORCOVFAILED;
}
if (((param[KAPPAI]<=0) || (param[KAPPAI]>1.0)) && (timespacedim==1)) {
strcpy(ERRORSTRING_OK,"0<kappa1<=1, dim=1");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAI]);
return ERRORCOVFAILED;
}
if (((param[KAPPAI]<=0) || (param[KAPPAI]>0.5)) && (timespacedim==2)) {
strcpy(ERRORSTRING_OK,"0<kappa1<=1/2, dim=2");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAI]);
return ERRORCOVFAILED;
}
if (param[KAPPAII]<=0) {
strcpy(ERRORSTRING_OK,"kappa2>0");
sprintf(ERRORSTRING_WRONG,"%f",param[KAPPAII]);
return ERRORCOVFAILED;
}
return 0;
}
static Real range_lgd1[8] = {0.0, 1.0, 0.01, 1.0,
0.0, RF_INF, 0.01, 20.0};
void rangelgd1(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
if (dim<=2) *index = -1; else *index=-2;
memcpy(range, range_lgd1, sizeof(Real) * 8);
if (dim==2) range[1]=range[3]=0.5;
}
/* FD model */
Real FD(Real *x,Real *p, int effectivedim){
static Real dold=RF_INF;
static Real kold, sk;
Real y, d, k, skP1;
d = p[KAPPAI];
y = fabs(*x);
k = (Real) trunc((double) y);
if (dold!=d || kold > k) {
sk = 1;
kold = 0.0;
}
// sign (-1)^k is (kold+d), 16.11.03, checked.
for (; kold<k; kold+=1.0) sk = sk * (kold + d) / (kold + 1.0 - d);
dold = d;
kold = k;
if (k == y) return sk;
skP1 = sk * (kold + d) / (kold + 1.0 - d);
return sk + (y - k) * (skP1 - sk);
}
SimulationType methodFD(int spacedim, bool grid){
return (grid ? CircEmbed : Direct);
}
int checkFD(Real *param, int timespacedim, SimulationType method) {
if (param[KAPPAI] < -0.5 || param[KAPPAI] >= 0.5) {
strcpy(ERRORSTRING_OK,"-0.5 <= kappa < 0.5");
sprintf(ERRORSTRING_WRONG, "%f",param[KAPPAI]);
return ERRORCOVFAILED;
}
if (timespacedim>1) {
strcpy(ERRORSTRING_OK, "dim=1");
sprintf(ERRORSTRING_WRONG, "%d", timespacedim);
return ERRORCOVFAILED;
}
return 0;
}
void rangeFD(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
if (dim==1) *index = -1; else *index=-2;
range[2] = (range[0] = -0.5) + RANGE_EPSILON;
range[4] = (range[1] = 0.5) - RANGE_EPSILON;
}
// ---------------------------------------------------------------------
/*
SimulationType method(int spacedim, bool grid){
if (grid) {
switch(spacedim) {
case 1 : return ;
case 2 : return ;
default: return ;
}
} else {
switch(spacedim) {
case 1 : return ;
case 2 : return ;
case 3 : return ;
default: return Direct;
}
}
}
void range(int dim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
Real r[] = {};
memcpy(range, r, sizeof(Real) * );
}
*/
Computing file changes ...
