swh:1:snp:8c443ee102a7f3a9fe0e834b4caec2f704d4d442
Tip revision: 874b82c3a0dd003dfee489f8e4cecf60c4149d18 authored by Martin Schlather on 14 July 2004, 00:00:00 UTC
version 1.1.13
version 1.1.13
Tip revision: 874b82c
MPPFcts.cc
/*
Authors
Martin Schlather, martin.schlather@cu.lu
*********** this function is still under construction *********
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.
*/
/******************************************************************************
*
* if deterministic then
* functions are always(!) normed to max(f)==1, i.e. for any parameter
* combination.
* ****** or shall they be normed anyhow ??
*
* the values are finally corrected in do_addmpp by
* sqrt(p[VARIANCE] / integral(f^2)) in case of Gaussian random fields
* (and integral(f) for the mean)
* 1.0 / integral(f_+) in case of max-stable processes
*
**************************************************************************** */
/*
IF NEW FUNCTION IS DEFINED MAKE SURE THAT NO MORE THAN 5 OR 6
ADDITIONALLY #defineD CONSTANTS ARE USED -- otherwise increase
c[MAXCOV][6] in RFsimu.h
*/
// basic shapes are: sphere (2 or 3 dimensions), cone, and Gaussian curve
#include <math.h>
#include <assert.h>
#include "RFsimu.h"
#include "MPPFcts.h"
static int mpp_dim;
static Real mpp_x[MAXDIM];
static Real spherical_Rsq, spherical_adddist;
Real spherical_model(Real *y)
{
Real distance;
int i;
distance = spherical_adddist;
for (i=0; i<mpp_dim; i++) {
register Real dummy;
dummy = y[i]-mpp_x[i];
distance += dummy * dummy;
}
if (distance>=spherical_Rsq) return 0;
return 1.0;
}
static Real cone_Rsq, /* square of radius */
cone_rsq, /* square of radius */
cone_height_socle,
cone_a, cone_b_socle; /* y = ax + b, for the slope */
Real cone_model(Real *y)
{
Real distance;
int i;
distance = 0.0;
for (i=0; i<mpp_dim; i++) {
register Real dummy;
dummy = y[i]-mpp_x[i];
distance += dummy * dummy;
}
if (distance>=cone_Rsq) return 0;
if (distance<=cone_rsq) return cone_height_socle;
return cone_b_socle + cone_a * sqrt(distance);
}
static Real gauss_invscale, gauss_effectiveradiussq;
Real gauss_model(Real *y)
{
Real distance;
int i;
distance = 0.0;
for (i=0; i<mpp_dim; i++) {
register Real dummy;
dummy = y[i]-mpp_x[i];
distance += dummy * dummy;
}
if (distance > gauss_effectiveradiussq) return 0;
return exp(-gauss_invscale * distance);
}
#define SPHERICALRsq 1
#define SPHERICALadddist 2
void generalspherical_init(mpp_storage* s, int v, int balldim)
{
int d;
Real gamma;
assert(s!=NULL);
s -> effectiveRadius[v] = 0.5 * s->param[v][SCALE];
if (s->addradius[v]<=0.0) {s->addradius[v] = s->effectiveRadius[v];}
s -> effectivearea [v]= 1.0;
for (d=0; d<s->timespacedim; d++)
s->effectivearea[v] *= (s->length[d] + 2.0 * s->effectiveRadius[v]);
for (d=s->timespacedim; d<balldim; d++) {
s->effectivearea[v] *= (2.0 * s->effectiveRadius[v]);
}
gamma = gammafn(1.0+0.5*balldim);
s->c[v][SPHERICALRsq] = s->effectiveRadius[v] * s->effectiveRadius[v];
s->integral[v] = s->integralsq[v] = s->integralpos[v] =
pow(SQRTPI * s->effectiveRadius[v], balldim) / gamma;
s->maxheight[v] = 1.0;
}
void circular_init(mpp_storage *s, int v)
{
int TRUESPACE=2;
generalspherical_init(s, v, TRUESPACE);
}
void circularMpp(mpp_storage *s, int v, Real *min, Real *max,
mppmodel *model )
{
int TRUESPACE = 2;
int d;
assert(s!=NULL);
mpp_dim = s->timespacedim;
spherical_Rsq = s->c[v][SPHERICALRsq];
*model = spherical_model;
for (d=0; d<s->timespacedim; d++) {
mpp_x[d]= s->min[d] - s->addradius[v] +
(s->length[d] + 2.0 * s->addradius[v]) * UNIFORM_RANDOM;
min[d] = mpp_x[d] - s->effectiveRadius[v];
max[d] = mpp_x[d] + s->effectiveRadius[v];
}
spherical_adddist = 0.0;
for (d=s->timespacedim; d<TRUESPACE; d++) {
Real dummy;
dummy = s->effectiveRadius[v] * UNIFORM_RANDOM;
spherical_adddist += dummy * dummy;
}
}
void spherical_init(mpp_storage *s, int v)
{
int TRUESPACE = 3;
generalspherical_init(s, v, TRUESPACE);
}
void sphericalMpp(mpp_storage *s, int v, Real *min, Real *max,
mppmodel *model )
{
int TRUESPACE = 3;
int d;
mpp_dim = s->timespacedim;
spherical_Rsq = s->c[v][SPHERICALRsq];
*model = spherical_model;
for (d=0; d<s->timespacedim; d++) {
mpp_x[d]= s->min[d] - s->addradius[v] +
(s->length[d] + 2.0 * s->addradius[v]) * UNIFORM_RANDOM;
min[d] = mpp_x[d] - s->effectiveRadius[v];
max[d] = mpp_x[d] + s->effectiveRadius[v];
}
spherical_adddist = 0.0;
for (d=s->timespacedim; d<TRUESPACE; d++) {
Real dummy;
dummy = s->effectiveRadius[v] * UNIFORM_RANDOM;
spherical_adddist += dummy * dummy;
}
}
#define CONErsq 1
#define CONERsq 2
#define CONEHEIGHTSOCLE 3
#define CONEA 4
#define CONEBSOCLE 5
void cone_init(mpp_storage *s, int v)
{
// ignoring: p[MEAN], p[NUGGET]
// note : sqrt(p[VARIANCE]) -> height
// p[SCALE] -> R
// p[KAPPAI] = r/R \in [0,1]
// p[KAPPAII]= socle (height)
// p[KAPPAIII]=height of cone (without socle)
//
// standard cone has radius 1/2, height of (potential) top: 1
//
int d;
Real cr,cb,CR,socle,height;
assert(s!=NULL);
socle = s->param[v][KAPPAII];
height= s->param[v][KAPPAIII];
cr= 0.5 * s->param[v][KAPPAI] * s->param[v][SCALE];
s->c[v][CONErsq] = cr * cr;
CR = 0.5 * s->param[v][SCALE];
s->c[v][CONERsq] = CR * CR;
s->c[v][CONEHEIGHTSOCLE] = socle + height;
s->c[v][CONEA] = /* a = h / (r - R) */
height / ((s->param[v][KAPPAI] - 1.0) * 0.5* s->param[v][SCALE]);
cb = height / (1.0 - s->param[v][KAPPAI]); /* b = h R / (R - r) */
s->c[v][CONEBSOCLE] = socle + cb;
switch (s->timespacedim) {
case 2 :
//
s->integral[v] = s->integralpos[v] =
PI * (s->c[v][CONERsq] * socle +
0.333333333333333333 * height *
(s->c[v][CONErsq] + s->c[v][CONERsq] + cr * CR));
//
s->integralsq[v]
= PI * (s->c[v][CONERsq] * socle * socle
+ s->c[v][CONErsq] * height * height
+ 0.5 * s->c[v][CONEA] * s->c[v][CONEA] *
(s->c[v][CONERsq] * s->c[v][CONERsq] -
s->c[v][CONErsq] * s->c[v][CONErsq])
+ 1.33333333333333333 * s->c[v][CONEA] * cb *
(s->c[v][CONERsq] * CR - s->c[v][CONErsq] * cr)
+ cb * cb * (s->c[v][CONERsq] - s->c[v][CONErsq])
);
//
break;
// case 3 : s -> integral = 4 PI / 3 * (R^3 * socle + r^3 * height)
// + 4 PI int_r^R r^2 a x + b) dx
// s -> integralsq = 4 PI / 3 * (R^3 * socle *socle
// + r^3 * height * height)
// + 4 PI int_r^R r^2 (a x + b)^2 dx
default : assert(false);
}
// here a "mean" effectiveRadius is defined + true radius for simulation;
// this is important
s->effectiveRadius[v] = CR;
if (s->addradius[v]<=0.0) {s->addradius[v] = s->effectiveRadius[v];}
// the calculation of the effective area can be very complicated
// if randomised scale is used with upper end point of the
// distribution ==oo.
// Then either the simulation is pretty approximate, and probably
// very bad (if a fixed, small border is used),
// or the simulation is very ineffective (if a large fixed border is used)
// or the formula might be very complicated (if there are adapted borders,
// and if the distribution density of the scale is weighted by about
// (length[1] + 2 scale[1]) *...*(length[dim] + 2 scale[dim])
//
s->effectivearea[v] = 1.0;
for (d=0; d<s->timespacedim; d++)
s->effectivearea[v]*= (s->length[d] + 2.0 * s->addradius[v]);
//
s->maxheight[v] = height + socle; //==s->c[CONEHEIGHTSOCLE]
}
// note. The cone_* values have to be reset all the time,
// since the function does not know whether there have
// been other simulations of random fields, as cone_init is called
// only at the very fist time.
// Furthermore, the below construction will match with the
// general construction that randomises parameters (e.g. scale mixtures)
void cone(mpp_storage *s, int v, Real *min, Real *max,
mppmodel *model )
{
int d;
mpp_dim = s->timespacedim;
cone_Rsq = s->c[v][CONERsq];
cone_rsq = s->c[v][CONErsq];
// note that c[CONEHEIGHT]==1 for Gaussian random fields !!
// However, for max-stable random fields s->height will have
// different values !!
cone_height_socle = s->c[v][CONEHEIGHTSOCLE];
cone_a = s->c[v][CONEA];
cone_b_socle = s->c[v][CONEBSOCLE];
*model = cone_model;
// logically the current effectiveRadius should be set here
// no need as here constant scale, so
// s -> effectiveRadius = 0.5 * s->param[v][SCALE];
// can be used.
// last: random point
for (d=0; d<mpp_dim; d++) {
mpp_x[d]= s->min[d] - s->addradius[v] +
(s->length[d] + 2.0 * s->addradius[v]) * UNIFORM_RANDOM;
min[d] = mpp_x[d] - s->effectiveRadius[v];
max[d] = mpp_x[d] + s->effectiveRadius[v];
}
}
int checkcone(Real *param, int timespacedim, SimulationType method) {
switch (method) {
case AdditiveMpp:
if (timespacedim!=2) {
strcpy(ERRORSTRING_OK,"2 dimensional space");
sprintf(ERRORSTRING_WRONG,"dim=%d",timespacedim);
return ERRORCOVFAILED;
}
break;
case Nothing :
break;
default :
strcpy(ERRORSTRING_OK,"method=\"add. MPP (random coins)\"");
strcpy(ERRORSTRING_WRONG,METHODNAMES[method]);
return ERRORCOVFAILED;
}
if ((param[KAPPAI] >= 1.0) ||
(param[KAPPAI] < 0.0) ||
(param[KAPPAII] < 0.0) ||
(param[KAPPAIII]< 0.0) ||
((param[KAPPAIII] + param[KAPPAII])==0.0)) {
strcpy(ERRORSTRING_OK,
"all kappa_i non-negative, kappa1<1, and kappa2+kappa3>0.0");
sprintf(ERRORSTRING_WRONG,"c(%f,%f,%f)",param[KAPPAI],
param[KAPPAII],param[KAPPAIII]);
return ERRORCOVFAILED;
}
return NOERROR;
}
SimulationType methodcone(int spacedim, bool grid){
return AdditiveMpp;
}
#define RANGE_EPSILON 1E-20
void rangecone(int spatialdim, int *index, Real* range){
// 2 x length(param) x {theor, pract }
*index = -1;
Real r[12] = {0, 1, 0, 1.0-RANGE_EPSILON,
0, RF_INF, RANGE_EPSILON, 10.0,
0, RF_INF, RANGE_EPSILON, 10.0};
memcpy(range, r, sizeof(Real) * 12);
}
#define GAUSSINVSCALE 0
#define GAUSSRADIUSSQ 1
Real gaussInt(int d, int xi, Real sigma, Real R) {
// int_{b(0,R) e^{-a r^2} dr = d b_d int_0^R r^{d-1} e^{-a r^2} dr
// where a = 2.0 * xi / sigma^2
// here : 2.0 is a factor so that the mpp function leads to the
// gaussian covariance model exp(-x^2)
// xi : 1 : integral ()
// 2 : integral ()^2
Real a;
a = 2.0 * xi / (sigma * sigma);
switch (d) {
case 1 :
Real sqrta;
sqrta = sqrt(a);
return SQRTPI / sqrta * (2.0 * pnorm(SQRT2 * sqrta * R, 0, 1, 1, 0) - 1.0);
case 2 :
return PI / a * (1.0 - exp(- a * R * R));
case 3 :
Real piDa;
piDa = PI / a;
return -2.0 * piDa * R * exp(- a * R * R )
+ piDa * sqrt(piDa) * (2.0 * pnorm(SQRT2 * sqrt(a) * R, 0, 1,1,0) - 1.0);
default : assert(false);
}
}
void gaussmpp_init(mpp_storage *s, int v)
{
int d;
s -> c[v][GAUSSINVSCALE] = 2.0 * s->param[v][INVSCALE] *
s->param[v][INVSCALE];
// e^{-2 x^2} !
s -> c[v][GAUSSRADIUSSQ] = - 0.5 * log(MPP_APPROXZERO)
* s->param[v][SCALE] * s->param[v][SCALE];
s -> effectiveRadius[v] = sqrt(s -> c[v][GAUSSRADIUSSQ]);
if (s->addradius[v]<=0.0) {s->addradius[v]=s->effectiveRadius[v];}
s->integral[v] = gaussInt(s->timespacedim, 1, s->param[v][SCALE],
s->effectiveRadius[v]);
s->integralsq[v] = gaussInt(s->timespacedim, 2, s->param[v][SCALE],
s->effectiveRadius[v]);
s->integralpos[v] = s->integral[v];
s->effectivearea[v] = 1.0;
for (d=0; d<s->timespacedim; d++)
s->effectivearea[v] *= (s->length[d] + 2.0 * s->addradius[v]);
s->maxheight[v]= 1.0;
}
void gaussmpp(mpp_storage *s, int v,Real *min, Real *max,
mppmodel *model )
{
int d;
gauss_invscale = s->c[v][GAUSSINVSCALE];
gauss_effectiveradiussq = s -> c[v][GAUSSRADIUSSQ];
*model = gauss_model;
mpp_dim = s->timespacedim;
for (d=0; d<mpp_dim; d++) {
mpp_x[d]= s->min[d] - s->addradius[v] +
(s->length[d] + 2.0 * s->addradius[v]) * UNIFORM_RANDOM;
min[d] = mpp_x[d] - s->effectiveRadius[v];
max[d] = mpp_x[d] + s->effectiveRadius[v];
}
}