https://github.com/cran/RandomFields
Tip revision: e10243fbd4eb0cbeaf518e67fbc5b8ad44889954 authored by Martin Schlather on 12 December 2019, 13:40:13 UTC
version 3.3.7
version 3.3.7
Tip revision: e10243f
operator.extremes.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
Definition of correlation functions and derivatives of hypermodels
Copyright (C) 2005 -- 2017 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 3
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.
*/
#include <R_ext/Applic.h>
#include "questions.h"
#include "operator.h"
#include "shape.h"
//////////////////////////////////////////////////////////////////////
// extremal gaussian
void extrgauss(double *x, model *cov, double *v) {
// BrownResnick to binary Gaussian
model *next = cov->sub[0];
double var, z;
COV(ZERO(next), next, &var);
COV(x, next, &z);
*v = 1.0 - SQRT(0.5 * (1.0 - z / var));
}
int check_extrgauss(model *cov) {
// to do extend to multivariate
model
*next = cov->sub[0];
double v;
int i,
vdim = VDIM0,
err = NOERROR;
if (VDIM0 != VDIM1) BUG;
if ((err = CHECK_PASSTYPE(next, PosDefType)) !=NOERROR) RETURN_ERR(err);
// if ((err = CHECK(next, cov->tsdim, cov->xdimprev, PosDefType,
// OWNDOM(0), OWNISO(0),
// SUBMODEL_DEP, cov->frame)) != NOERROR) RETURN_ERR(err);
setbackward(cov, next);
for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = 1.0;
COV(ZERO(next), next, &v);
if (v != 1.0) SERR("only correlation functions allowed");
RETURN_NOERROR;
}
//////////////////////////////////////////////////////////////////////
// Brown Resnick
#define BR_FACTOR 0.25 // KSH: '/ 4'
#define BR_SEMI_FACTOR (2 * BR_FACTOR) // hier Semi-Variogram
void brownresnick(double *x, model *cov, double *v) {
// BrownResnick to binary Gaussian
model *next = cov->sub[0];
double z;
COV(ZERO(next), next, &z);
COV(x, next, v);
*v = 2.0 * pnorm(SQRT((z - *v) * BR_SEMI_FACTOR), 0.0, 1.0, false, false);
}
void Dbrownresnick(double *x, model *cov, double *v) {
// b = BR_SEMI_FACTOR
// gamma(h) = C(0) - C(h)
// varphi(SQRT(gamma * b)) * (b C') / SQRT(gamma * b)
// varphi density standard normal
// BrownResnick to binary Gaussian
model *next = cov->sub[0];
double z, abl, s;
//
assert((hasMaxStableFrame(cov) || hasGaussMethodFrame(cov) ||
hasAnyEvaluationFrame(cov)) && cov->taylorN > 1);
if (cov->taylor[1][TaylorPow] == 0.0) { // const??
*v = 0.0;
return;
}
if (*x != 0.0) {
COV(ZERO(next), next, &z);
COV(x, next, v);
assert(DefList[NEXTNR].D != NULL);
assert(DefList[NEXTNR].D != ErrD);
Abl1(x, next, &abl); // Vorsicht: abl = -\gamma'
// printf("x=%10g abl=%10g %10g z=%10g %10g\n", *x, abl, BR_SEMI_FACTOR, z, *v);
// PMI0(next);
// crash();
abl *= BR_SEMI_FACTOR;
s = SQRT((z - *v) * BR_SEMI_FACTOR); // SQRT(c * gamma)
*v = dnorm(s, 0.0, 1.0, false) * abl / s; // -\varphi * \gamma' / sqrt\gamma
// =-2 \varphi * (-C') / (2 sqrt\gamma)
// mit C= c_0 -\gamma
//printf("v=%10g %10g %10g\n", *v, abl, s);
assert(*v <= 0);
} else {
if (cov->taylor[1][TaylorPow] < 1.0) {
*v = RF_NEGINF;
} else if (cov->taylor[1][TaylorPow] == 1.0) {
*v = FABS(cov->taylor[1][TaylorConst]);
assert(*v > 0.0);
} else BUG;
}
// 2 * c * varphi(s) * gamma' / s
}
void DDbrownresnick(double *x, model *cov, double *v) {
// D = \varphi * b C' / SQRT(\gamma b)
// b = BR_SEMI_FACTOR
// gamma(h) = C(0) - C(h)
// varphi(SQRT(gamma * b)) [(b C')^2 / [2 sqrt (gamma * b)]
// +(b C'') / sqrt (gamma * b)
// +1/2 * (b C')^2 / sqrt (gamma * b)^3 ]
// varphi density standard normal
model *next = cov->sub[0];
double z, abl, abl2, s0, s;
if (cov->taylor[1][TaylorPow] == 0.0) {
*v = 0.0;
return;
}
if (*x != 0.0) {
COV(ZERO(next), next, &z);
COV(x, next, v);
Abl1(x, next, &abl);
Abl2(x, next, &abl2);
s0 = (z - *v) * BR_SEMI_FACTOR;
s = SQRT(s0); // SQRT(c * gamma)
abl *= BR_SEMI_FACTOR;
abl2 *= BR_SEMI_FACTOR;
*v = dnorm(s, 0.0, 1.0, false) / s * (abl2 + abl * abl * 0.5 * (1/s0 + 1));
assert(*v >= 0);
} else {
*v = cov->taylor[1][TaylorPow]==1 ? 0.0 : RF_INF;
}
}
void D3brownresnick(double *x, model *cov, double *v) {
// D = \varphi * b C' / SQRT(\gamma b)
// b = BR_SEMI_FACTOR
// gamma(h) = C(0) - C(h)
// varphi(SQRT(gamma * b)) [(b C')^3 / [4 sqrt (gamma * b)]
// +3 (b C') (b C'')/ [2 sqrt (gamma * b)]
// + (b C''') / [sqrt (gamma * b)]
// + (b C')^2 / [2 sqrt (gamma * b)]
// +3(b C') (b C'') / [2 sqrt (gamma * b)^3]
// +3(b C')^3 / [4 sqrt (gamma * b)^5]
model *next = cov->sub[0];
double z, abl, abl2, abl3, s0, s;
if (cov->taylor[1][TaylorPow] == 0.0) {
*v = 0.0;
return;
}
if (*x != 0.0) {
COV(ZERO(next), next, &z);
COV(x, next, v);
Abl1(x, next, &abl);
Abl2(x, next, &abl2);
Abl3(x, next, &abl3);
s0 = (z - *v) * BR_SEMI_FACTOR;
s = SQRT(s0); // SQRT(c * gamma)
abl *= BR_SEMI_FACTOR;
abl2 *= BR_SEMI_FACTOR;
abl3 *= BR_SEMI_FACTOR;
*v = dnorm(s, 0.0, 1.0, false) / s *
(abl3 +
1.5 * abl2 * abl * (1/s0 + 1) +
abl * abl * abl * (0.25 + 0.5 / s0 + 0.75 / (s0 * s0)));
// printf("br x=%10g v=%10e s=%10g abl=%10g s0=%10g abl2=%10g\n", *x, *v, s, abl, s0, abl2);
// assert(*v >= 0);
} else {
*v = cov->taylor[1][TaylorPow]==1 ? 0.0 : RF_NEGINF; // da zweite Ableitung im Ursprung Pol +\infty hat.
}
}
int TaylorBrownresnick(model *cov) {
model
*next = cov->sub[0];
int idx = isnowPosDef(next); // Taylorentw. in 2 Glieder falls pos def.
assert(idx == 0);
// assert(!idx);
if (next->taylor[idx][TaylorPow] >= 2.0) {//for pos/neg def only == 2 possible
cov->full_derivs = 1;
} else cov->full_derivs = 0;
cov->rese_derivs = MIN(3, next->rese_derivs);
// else if (next->taylor[idx][TaylorPow] == 2) {
// if (next->taylorN < 2 + idx) cov->rese_derivs = 0;
// else if (cov->rese_derivs > 2) cov->rese_derivs = 2;
// }
if (next->taylorN >= 1 + idx && next->taylor[idx][TaylorConst] < 0.0) {
// 2 \Phi(SQRT(gamma * br_f)) =
//1+ 2 * phi(0) * SQRT(gamma * br_f) - 2 phi(0) / 6 *SQRT(gamma*br_f)^3
//sei gamma(x) = c x^alpha + d x^beta, beta > alpha. Dann gilt
// also 2 \Phi(SQRT(gamma * br_f)) ~
// 1 + 2 * phi(0) * SQRT((c x^alpha + d x^beta)* br_f)
// - 2 * phi(0) / 6 * SQRT((c x^alpha + d x^beta) *br_f)^3
// = 1 + 2 * phi(0)[ (c f x^a)^{1/2} (1 + 0.5 d / c * x^{b-a})
// (c f x^a)^{3/2} (1 + 1.5 d / c * x^{b-a}) ]
// ~ 1 + 2 * phi(0) (c f x^a)^{1/2}
// + phi(0)(c f)^{1/2} d / c * x^{b - a/2}) (*)
// + 2 * phi(0) (c f x^a)^{3/2} (**)
//
//
// da a <= 2 und b > a gilt:
// 1. Abltg ~ phi(0) (c f)^{1/2} x^{a/2 - 1}
// + phi(0) (c f)^{1/2} d / c * (b - a/2) x^{b - a/2 - 1}
// + 3 phi(0 (c f)^{3/2} a x^{3/2 * a - 1}
// 2. Abltg ~ phi(0) (c f)^{1/2} (a/2 - 1) x^{a/2 - 2} ~ infty, a != 2
// und fuer a=2 gilt folgendes: da gamma neg def ist e^{-\gamma} pos def
// und mit sasvari folgt aus der Taylor entwicklung dass (a=2) < b <= 4).
// Andererseits dominiert x{3/2 a} den Term x^{b - a/2} falls
// 3/2 a < b - a/2 <=> b > 2a = 4.
// Also ist b=4 ein Grenzfall, wo beide (*) und (**) eine Rolle spielen.
// Ansonsten nur (*).
//
assert(next->taylor[idx][TaylorConst] <= 0.0);
cov->taylorN = 2;
cov->taylor[0][TaylorConst] = 1.0;
cov->taylor[0][TaylorPow] = 0.0;
double
next_taylor_const = - next->taylor[idx][TaylorConst],
g = SQRT(next_taylor_const * BR_SEMI_FACTOR * 0.5 / M_PI);
cov->taylor[1][TaylorConst] = - 2 * g;
cov->taylor[1][TaylorPow] = 0.5 * next->taylor[idx][TaylorPow];
if (next->taylor[idx][TaylorPow] == 2) {
if (next->taylorN >= 2 + idx) {
cov->taylorN = 3;
if (next->taylor[idx + 1][TaylorConst] != 0) {
cov->taylor[2][TaylorConst] =
g * next->taylor[idx + 1][TaylorConst] / next_taylor_const;
cov->taylor[2][TaylorPow] =
next->taylor[idx + 1][TaylorPow] - 0.5*next->taylor[idx][TaylorPow];
} else { // Spezialfall:L fractionelle BR
cov->taylor[2][TaylorConst] = 0.0;
cov->taylor[2][TaylorPow] = 4.0;
}
if (next->taylor[idx + 1][TaylorPow] == 4) {
cov->taylor[1][TaylorConst] += 2 * g * next_taylor_const * BR_SEMI_FACTOR;
}
} else cov->taylorN = 0;
}
} else cov->taylorN = 0;
if (next->tailN >= 1) {
cov->tailN = 1;
cov->tail[0][TaylorPow] = -0.5 * next->tail[0][TaylorPow];
if (next->tail[0][TaylorPow] > 0) {
assert( next->tail[0][TaylorConst] < 0.0);
double next_tail_const = - next->tail[0][TaylorConst];
cov->tail[0][TaylorConst] =
2.0 / SQRT(2.0 * M_PI * BR_SEMI_FACTOR * next_tail_const);
cov->tail[0][TaylorExpConst] = 0.5 * BR_SEMI_FACTOR * next_tail_const;
cov->tail[0][TaylorExpPow] = next->tail[0][TaylorPow];
} else {
cov->tail[0][TaylorConst] =
2.0 / SQRT(2.0 * M_PI * BR_SEMI_FACTOR * next->tail[0][TaylorConst])
* EXP(-0.5 * BR_SEMI_FACTOR * next->tail[0][TaylorConst]);
cov->tail[0][TaylorPow] = cov->tail[0][TaylorExpConst] =
cov->tail[0][TaylorExpPow] = 0.0;
}
} else cov->tailN = 0;
if (cov->taylorN < 1) cov->rese_derivs = 0;
RETURN_NOERROR;
}
int checkbrownresnick(model *cov) {
// to do extend to multivariate
model
*next = cov->sub[0];
int i, err,
vdim = VDIM0;
if (vdim != VDIM1) BUG;
if ((err = CHECK_PASSTF(next, VariogramType, SUBMODEL_DEP,
// hasMaxStableFrame(cov) ? BrMethodType :
EvaluationType))
// if ((err = CHECK(next, dim, dim, VariogramType, OWNDOM(0),
// OWNISO(0), SUBMODEL_DEP,
// hasMaxStableFrame(cov) ? Ma xStableType : EvaluationType))
!= NOERROR) {
RETURN_ERR(err);
}
setbackward(cov, next);
cov->monotone = isBernstein(next) ? GNEITING_MON :
isMonotone(next) ? MONOTONE : NOT_MONOTONE;
if ((err = TaylorBrownresnick(cov)) != NOERROR) RETURN_ERR(err);
for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = 1.0;
MEMCOPY(cov->pref, DefList[COVNR].pref, sizeof(pref_shorttype));
RETURN_NOERROR;
}
int struct_brownresnick(model *cov, model VARIABLE_IS_NOT_USED **newmodel) {
model *next = cov->sub[0];
if (hasSmithFrame(cov)) {
if (1 > next->taylorN || 1 > next->tailN)
SERR2("frame '%.50s' not possible for submodel '%.50s'",
TYPE_NAMES[cov->frame], NICK(next));
BUG;
// shape ist somit die Ableitung, falls d=1 und i.w. die
// zweifache Ableitung, falls d=3
// hier ist auch Taylor zu setztn fuer den neuen Shape
} else ILLEGAL_FRAME;
RETURN_NOERROR;
}
int init_brownresnick(model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
int err;
// model *next = cov->sub[0];
if ((err = TaylorBrownresnick(cov)) != NOERROR) RETURN_ERR(err);
RETURN_NOERROR;
}
void do_brownresnick(model *cov, gen_storage *s) {
model *next = cov->sub[0];
DO(next, s); // nicht gatternr
}
void BR2EG(double *x, model *cov, double *v) {
// BrownResnick to binary Gaussian
model *next = cov->sub[0];
double z;
COV(ZERO(next), next, &z);
COV(x, next, v);
z = 2.0 * pnorm(SQRT( (z - *v) * BR_SEMI_FACTOR), 0.0, 1.0, true, false) -1.0;
*v = 1.0 - 2.0 * z * z;
}
int check_BR2EG(model *cov) {
// to do extend to multivariate
model
*next = cov->sub[0];
double v, t, alpha;
int err, i,
vdim = VDIM0;
if (VDIM0 != VDIM1) BUG;
if ((err = CHECK_PASSTYPE(next, PosDefType)) !=NOERROR) RETURN_ERR(err);
// if ((err = CHECK(next, cov->tsdim, cov->xdimown, PosDefType,
// OWNDOM(0), OWNISO(0),
// SCALAR, cov->frame)) != NOERROR) RETURN_ERR(err);
setbackward(cov, next);
for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = 1.0;
if (next->pref[Nothing] == PREF_NONE) RETURN_ERR(ERRORPREFNONECOV);
// Erfc(x) = 2(1 - Phi(x * SQRT(2)))
// Erf(x) = 2 Phi(x * SQRT(2)) - 1
// Erf^{-1}(x) = Phi^{-1}( (y + 1) / 2 ) / SQRT(2)
// Sei c = 1-2 * Erf(SQRT(semivario / 4))^2 .
// Also c = 1 - 2 [ 2 Phi(SQRT(semivario / 2)) - 1]^2
// Umgekehrt semivario = 4 * { Erf^{-1} (sqrt[0.5 (1 - c)]) } ^2
// mit c = 0 folgt SQRT(0.5 (1-c)) = 1 / SQRT(2)
// semivario = 2 * {Phi^{-1}( (1 / SQRT(2) + 1) / 2 ) } ^2
alpha = 0.0;
COV(ZERO(next), next, &v);
t = qnorm(0.5 * (1.0 + INVSQRTTWO), 0.0, 1.0, true, false);
t *= t / (BR_SEMI_FACTOR * (1.0 - alpha)); // 1/2 wegen Erf->qnorm
if (v > t)
SERR2("variance equals %10g, but must be at most 4(Erf^{-1}(1/2))^2 = %10g",
v, t);
RETURN_NOERROR;
}
void BR2BG(double *x, model *cov, double *v) {
// BrownResnick to binary Gaussian
model *next = cov->sub[0];
double z;
COV(ZERO(next), next, &z);
COV(x, next, v);
z = 2.0 * pnorm(SQRT( (z - *v) * BR_SEMI_FACTOR), 0.0, 1.0, true, false) -1;
*v = COS(M_PI * z);
}
int check_BR2BG(model *cov) {
// to do extend to multivariate
model *next = cov->sub[0];
double v, t, alpha;
int err, i,
vdim = VDIM0;
if (VDIM0 != VDIM1) BUG;
// if ((err = CHECK(next, cov->tsdim, cov->xdimown, PosDefType,
// OWNDOM(0), OWNISO(0),
// SCALAR, cov->frame)) != NOERROR) RETURN_ERR(err);
if ((err = CHECK_PASSTF(next, PosDefType, vdim, //cov->frame
//hasMaxStableFrame(cov) ? BrMethodType :
EvaluationType
)) != NOERROR)
RETURN_ERR(err);
setbackward(cov, next);
for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = 1.0;
if (next->pref[Nothing] == PREF_NONE) RETURN_ERR(ERRORPREFNONECOV);
// Erfc(x) = 2(1 - Phi(x * SQRT(2)))
// Erf(x) = 2 Phi(x * SQRT(2)) - 1
// Erf^{-1}(x) = Phi^{-1}( (y + 1) / 2 ) / SQRT(2)
// Sei c = COS(pi * Erf(SQRT(semivario / 4))) .
// Also c = COS(pi * [2 * Phi(SQRT(semivario / 2)) - 1] )
// Umgekehrt semivario = 2 * { Phi^{-1}(0.5 * [arcCOS( c ) / pi + 1]) }^2
// mit c = 0 folgt arcCOS(c)/ pi + 1 = 3/2
// semivario = 2 * { Phi^{-1}( 3 / 4) }^2
COV(ZERO(next), next, &v);
alpha = 0.0; // to do
t = qnorm(0.75, 0.0, 1.0, false, false);
t *= t / (BR_SEMI_FACTOR * (1.0 - alpha)); // 1/2 wegen Erf->qnorm
if (v > t) {
SERR2("variance equals %10g, but must be at most 4(Erf^{-1}(1 / 2))^2 = %10g",
v, t);
}
RETURN_NOERROR;
}
void strokorb(double *x, model *cov, double *v) {
// BrownResnick to strokorb Gaussian
model *next = cov->sub[0];
int dim = OWNLOGDIM(0);
double u;
int idx = 0;
u = 2 * *x;
switch (dim) {
case 1 :
Abl1(&u, next, v);
*v = - *v;
break;
case 3 :
if (*x == 0.0) {
while (idx < next->taylorN && (next->taylor[idx][TaylorPow] ==0.0 ||
next->taylor[idx][TaylorPow] ==1.0)) idx++;
if (idx >= next->taylorN) BUG;
double p = next->taylor[idx][TaylorPow];
if (p > 3.0) BUG;
// > 3 ist mathematisch vermutlich nicht moeglich, da bei strokorb
// das submodel die entwicklung 1 - c x (oder rauher) hat
if (p == 3.0)
*v = next->taylor[idx][TaylorConst] * p * (p - 1) * POW(2, p-2)/ M_PI;
else *v = RF_INF;
} else {
Abl2(&u, next, v);
*v /= (M_PI * *x);
}
break;
default: BUG;
}
if (*v<0) {
BUG;
}
// printf("x=%10g %10g ", *x, *v);
assert(*v >= 0);
assert(*x >= 0);
// assert(false);
}
int TaylorStrokorb(model *cov) {
model *next = cov->sub[0];
int dim = OWNLOGDIM(0);
if (next->tailN < 1)
SERR2("%d members of the Taylor expansion at infinity of '%.50s', but at least order 1 required.", next->tailN, NICK(next));
cov->taylorN = cov->tailN = 1;
cov->tail[0][TaylorExpConst] = next->tail[0][TaylorExpConst];
cov->tail[0][TaylorExpPow] = next->tail[0][TaylorExpPow];
switch(dim) {
case 1 :
cov->taylor[0][TaylorConst] =
-next->taylor[1][TaylorConst] * next->taylor[1][TaylorPow];
cov->taylor[0][TaylorPow] = next->taylor[1][TaylorPow] - 1.0;
cov->taylorN = 1;
cov->tailN = 1;
if (next->tail[0][TaylorExpPow] != 0.0) {
if (next->tail[0][TaylorExpConst] == 0.0) BUG;
cov->tail[0][TaylorConst] = next->tail[0][TaylorConst]
* next->tail[0][TaylorExpConst] * next->tail[0][TaylorExpPow];
cov->tail[0][TaylorPow] =
next->tail[0][TaylorPow] + next->tail[0][TaylorExpPow] - 1;
} else {
if (next->tail[0][TaylorExpConst] != 0.0) BUG;
if (next->tail[0][TaylorPow] >= 0) {
if (next->tail[0][TaylorConst] == 0) SERR("trivial submodel");
SERR1("'%.50s' is not integrable", NICK(next));
}
cov->tail[0][TaylorConst] =
- next->tail[0][TaylorConst] * next->tail[0][TaylorPow];
cov->tail[0][TaylorPow] = next->tail[0][TaylorPow] - 1.0;
}
break;
case 3 : {
int idx = 1;
if (next->taylorN <= 1)
SERR2("%d members of the Taylor expansion of '%.50s' known, but at least 2 members required.", next->taylorN, NICK(next));
if (next->taylor[idx][TaylorPow] == 1.0) {
if (next->taylorN <= idx + 1)
SERR3("%d members of the Taylor expansion of '%.50s' known, but at least %d members required.", next->taylorN, NICK(next), idx + 1);
idx++;
} else assert(next->taylor[idx][TaylorPow] < 1.0);
cov->taylor[0][TaylorPow] = (next->taylor[idx][TaylorPow] - 2.0) - 1.0;
cov->taylor[0][TaylorConst] = next->taylor[idx][TaylorConst] / M_PI *
next->taylor[idx][TaylorPow] * (next->taylor[idx][TaylorPow] - 1) *
POW(2.0, cov->taylor[0][TaylorPow]); //
cov->taylorN = 1;
cov->tailN = 1;
if (next->tail[0][TaylorExpPow] != 0.0) {
if (next->tail[0][TaylorExpConst] == 0.0) BUG;
// Achtung ! Reihenfolge der Berechnung!
double f = next->tail[0][TaylorExpConst] * next->tail[0][TaylorExpPow];
cov->tail[0][TaylorPow] = next->tail[0][TaylorPow] +
2.0 * (next->tail[0][TaylorExpPow] - 1.0);
cov->tail[0][TaylorConst]= next->tail[0][TaylorConst] * f * f *
POW(2.0, cov->tail[0][TaylorPow]);
cov->tail[0][TaylorPow] -= 1.0;
cov->tail[0][TaylorExpConst] *= POW(2.0, cov->tail[0][TaylorExpPow]);
} else {
// Achtung ! Reihenfolge der Berechnung!
if (next->tail[0][TaylorExpConst] != 0.0) BUG;
cov->tail[0][TaylorPow] = next->tail[0][TaylorPow] - 2.0;
cov->tail[0][TaylorConst] = next->tail[0][TaylorConst] / M_PI
* next->tail[0][TaylorPow] * (next->tail[0][TaylorPow] - 1.0)
* POW(2.0, next->tail[0][TaylorPow] - 2.0);
cov->tail[0][TaylorPow] -= 1.0;
}
}
break;
default: BUG;
}
RETURN_NOERROR;
}
int checkstrokorb(model *cov) {
model *next = cov->sub[0];
int
dim = OWNLOGDIM(0),
err = NOERROR;
if ((err = CHECK_PASSTF(next, TcfType, SCALAR, EvaluationType)) != NOERROR)
RETURN_ERR(err);
// if ((err = CHECK(next, cov->tsdim, cov->xdimprev, TcfType,
// OWNDOM(0), OWNISO(0),
// SCALAR, EvaluationType)) != NOERROR) RETURN_ERR(err);
if (next->randomkappa) RETURN_ERR(ERRORRANDOMKAPPA);
if (!isGneiting(next))
SERR("member of the Gneiting-Schaback class as submodel needed");
switch(dim) {
case 1:
if (next->rese_derivs < 1) {
//APMI(cov);
SERR("submodel must be once differentiable");
}
break;
case 3:
if (next->rese_derivs < 2) SERR("submodel must be twice differentiable");
break;
default:
SERR("only dimensions 1 and 3 are allowed");
}
if (false && !(hasMaxStableFrame(cov) || hasRandomFrame(cov))) {
//APMI(cov); crash(); crash();
SERR1("'%.50s' may be used only as a shape function with max-stable field simulation", NICK(cov));
}
if ((err = TaylorStrokorb(cov)) != NOERROR) RETURN_ERR(err);
setbackward(cov, next);
RETURN_NOERROR;
}
int init_strokorb(model *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
int err;
if (hasSmithFrame(cov) || hasRandomFrame(cov)) {
//model *next = cov->sub[0];
// double
// radius = P(TRUNC_RADIUS)[0]; // default -1
// if ((err = INIT(next, 0, s)) != NOERROR) RETURN_ERR(err);
//if (radius>=0 && radius < cov->mpp.refradius) cov->mpp.refradius = radius;
// Eplus, M2 are assumed to be still precise !!
assert(VDIM0 == 1 && VDIM1 == 1);
cov->mpp.maxheights[0] = RF_NA;
}
else ILLEGAL_FRAME;
/*
if ((err = CHECK(cov->sub[0], cov->tsdim, cov->xdimprev,
ShapeType, //
OWNDOM(0), OWNISO(0),
SCALAR,
M axStableType// otherwise do will fail
)) != NOERROR) RETURN_ERR(err);
*/
cov->mpp.maxheights[0] = 1.0; // all with be covered by pgs
if (cov->mpp.moments >= 1) {
cov->mpp.mM[1] = cov->mpp.mMplus[1] = 1;
}
if ((err = TaylorStrokorb(cov)) != NOERROR) RETURN_ERR(err);
RETURN_NOERROR;
}
void do_strokorb(model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
BUG;
}
////////////////// strokorbBall
#define STROKORBBALL_DIM 0 // Inner
int checkstrokorbBall(model *cov) {
model
//*key = cov->key,
*next = cov->sub[0];
int
dim = OWNLOGDIM(0),
err = NOERROR;
assert(cov->key == NULL);
if ((err = CHECK_PASSTF(next, TcfType, SCALAR, EvaluationType)) !=NOERROR) RETURN_ERR(err);
// if ((err = CHECK(next, dim, cov->xdimprev, TcfType,
// OWNDOM(0), OWNISO(0),
// SCALAR, EvaluationType)) != NOERROR) RETURN_ERR(err);
if (!isGneiting(next))
SERR("member of the Gneiting-Schaback class as submodel needed");
switch(dim) {
case 1:
if (next->rese_derivs < 2) SERR("submodel must be twice differentiable");
break;
case 3:
if (next->rese_derivs < 3)
SERR("submodel must be three times differentiable");
break;
default:
SERR("only dimensions 1 and 3 are allowed");
}
if (false && !(hasMaxStableFrame(cov) || hasRandomFrame(cov)))
SERR1("'%.50s' may be used only as a shape function with max-stable field simulation", NICK(cov));
if (next->tailN < 1)
SERR2("%d members of the Taylor expansion at infinity of '%.50s' found, but at least 1 is required.", next->tailN, NICK(next));
if (next->taylorN < 2)
SERR2("%d members of the Taylor expansion of '%.50s' found, but at least 2 is required.", next->taylorN, NICK(next));
setbackward(cov, next);
RETURN_NOERROR;
}
void ScaleToVar(model *local, model *remote,
int VARIABLE_IS_NOT_USED variant) {
// ACHTUNG!! from and to sind hier vertauscht!!
assert(MODELNR(local)==POWER_DOLLAR && MODELNR(remote)==LOC);
// int dim = local->tsdim;
double scale = PARAM0(local, POWSCALE);
// scale im entfernten model setzen
PARAM(remote, LOC_SCALE)[0] = scale;
remote->randomkappa = true;
}
int struct_strokorbBall(model *cov, model **newmodel) {
int err,
dim = OWNLOGDIM(0);
ASSERT_NEWMODEL_NOT_NULL;
if (hasSmithFrame(cov)) {
addModel(newmodel, BALL, cov);
addModel(newmodel, POWER_DOLLAR);
kdefault(*newmodel, POWSCALE, 1.0);
kdefault(*newmodel, POWPOWER, -dim);
kdefault(*newmodel, POWVAR, 1.0 / VolumeBall(dim, BALL_RADIUS));
model *pts=NULL, *scale=NULL;
if ((err = covcpy(&pts, *newmodel)) != NOERROR) RETURN_ERR(err);
if (DefList[COVNR].kappas < 2) {
if ((err = covcpy(&scale, cov)) != NOERROR) RETURN_ERR(err);
// !! inverse scale gegenueber paper
SET_NR(scale, STROKORB_BALL_INNER);
kdefault(scale, STROKORBBALL_DIM, dim);
addModel(&scale, RECTANGULAR, *newmodel);
kdefault(scale, RECT_APPROX, false);
kdefault(scale, RECT_ONESIDED, true);
(*newmodel)->kappasub[POWSCALE] = scale;
} else { // for testing only
addModelKappa(*newmodel, POWSCALE, UNIF);
kdefault((*newmodel)->kappasub[POWSCALE], UNIF_MIN, P0(0));
kdefault((*newmodel)->kappasub[POWSCALE], UNIF_MAX, P0(1));
}
addModel(&pts, RECTANGULAR, *newmodel);
addModel(&pts, LOC, *newmodel);
kdefault(pts, LOC_SCALE, 1.0);
kdefault(pts, LOC_POWER, -dim);
addModelKappa(pts, LOC_SCALE, NULL_MODEL);
kdefault(pts->kappasub[LOC_SCALE], NULL_TYPE, RandomType);
addSetParam(newmodel, pts, ScaleToVar, true, 0);
addModel(newmodel, ZHOU); // to do : unif better ?!
(*newmodel)->sub[PGS_LOC] = pts;
SET_CALLING(pts, *newmodel);
} else ILLEGAL_FRAME_STRUCT;
RETURN_NOERROR;
}
void rangestrokorbball(model VARIABLE_IS_NOT_USED *cov, range_type *range){
range->min[0] = RF_NEGINF;
range->max[0] = RF_INF;
range->pmin[0] = -4.0;
range->pmax[0] = 4.0;
range->openmin[0] = false;
range->openmax[0] = false;
range->min[1] = RF_NEGINF;
range->max[1] = RF_INF;
range->pmin[1] = -4.0;
range->pmax[1] = 4.0;
range->openmin[1] = false;
range->openmax[1] = false;
}
/*
int init_strokorbBall(model *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
if (cov->frame = = M axStableType) {
//model *next = cov->sub[0];
// double
// radius = P(TRUNC_RADIUS)[0]; // default -1
// if ((err = INIT(next, 0, s)) != NOERROR) RETURN_ERR(err);
//if (radius>=0 && radius < cov->mpp.refradius) cov->mpp.refradius = radius;
// Eplus, M2 are assumed to be still precise !!
cov->mpp.maxheight = RF_NA;
}
else ILLEGAL_FRAME;
cov->mpp.maxheight = 1.0; // all with be covered by pgs
if (cov->mpp.moments >= 1) {
cov->mpp.mM[1] = cov->mpp.mMplus[1] = 1;
}
RETURN_NOERROR;
}
void do_strokorbBall(model VARIABLE_IS_NOT_USED *cov, gen_storage VA// RIABLE_IS_NOT_USED *s) {
// assert(false);
// }
// */
void strokorbBallInner(double *x, model *cov, double *v) {
// only proportional to a density !!!
model *next = cov->sub[0];
int dim = COVNR != STROKORB_BALL_INNER || PisNULL(STROKORBBALL_DIM)
? OWNLOGDIM(0) : P0INT(STROKORBBALL_DIM);
if (*x <= 0) {
*v = 0.0;
return;
}
double y = 2 * *x;
switch (dim) {
case 1 :
// \chi''(s)*s [ diameter ] -> 4\chi''(2 s)*s
Abl2(&y, next, v);
*v *= 2.0 * y;
break;
case 3 :
double w;
Abl2(&y, next, v);
Abl3(&y, next, &w);
*v = 2.0 * (*v - w * y) * y / 3.0;
break;
default: BUG;
}
if (*v<0) {
BUG;
}
// assert(*v >= 0);
}
int TaylorBall(model *cov) {
model *next = cov->sub[0];
if (next->tailN < 1 || next->taylorN < 2)
SERR1("taylor expansions of '%.50s' not programmed yet", NICK(next));
double tep = next->tail[0][TaylorExpPow],
tp = next->tail[0][TaylorPow];
cov->taylorN = cov->tailN = 1;
cov->tail[0][TaylorExpConst] = POW(2.0, tep) * next->tail[0][TaylorExpConst];
cov->tail[0][TaylorExpPow] = tep;
int idx = 1;
if (next->taylor[1][TaylorPow] == (int) next->taylor[1][TaylorPow]) {
assert(next->taylor[1][TaylorPow] == 1.0);
// 2 sollte nie auftreten, laeuft aber gleich
if (next->taylorN >= 3) idx++;
else SERR1("%.50s does not have a long enough taylor development programmed",
NICK(next));
}
double TP = next->taylor[idx][TaylorPow];
switch(P0INT(STROKORBBALL_DIM)) {
case 1 :
if (tep != 0.0) {
cov->tail[0][TaylorPow] = tp + 2.0 * (tep - 1.0) + 1.0;// !! +1, da noch ein x drauf-multipliziert wird
cov->tail[0][TaylorConst] = next->tail[0][TaylorExpConst] * tep;
cov->tail[0][TaylorConst] *= cov->tail[0][TaylorConst];
} else {
cov->tail[0][TaylorConst] = tp * (tp -1.0);
cov->tail[0][TaylorPow] = tp - 1.0; // !! nicht -2.0, da noch ein x drauf-multipliziert wird
}
cov->taylor[0][TaylorConst] = TP * (TP - 1.0);
cov->taylor[0][TaylorPow] = TP - 1.0;
break;
case 3 :
if (tep != 0.0) {
double dummy = next->tail[0][TaylorExpConst] * tep;
cov->tail[0][TaylorConst] = dummy * dummy * dummy / 3.0;
cov->tail[0][TaylorPow] = tp + 3 * tep - 1.0;
} else {
cov->tail[0][TaylorConst] = tp * (tp - 1.0) * (3.0 - tp) / 3.0;
cov->tail[0][TaylorPow] = tp - 1.0; // !! nicht -2.0, da noch ein x drauf-multipliziert wird
}
cov->taylor[0][TaylorConst] = TP * (TP - 1.0) * (3.0 - TP) / 3.0;
cov->taylor[0][TaylorPow] = TP - 2.0;
break;
default: BUG;
}
cov->tail[0][TaylorConst] *= 2.0 * next->tail[0][TaylorConst] *
POW(2.0, cov->tail[0][TaylorPow]);
cov->taylor[0][TaylorConst] *= 2.0 * next->taylor[idx][TaylorConst] *
POW(2.0, cov->taylor[0][TaylorPow]);
RETURN_NOERROR;
}
int check_strokorbBallInner(model *cov) {
model *next = cov->sub[0];
int err;
if ((err = checkkappas(cov)) != NOERROR) RETURN_ERR(err);
if (OWNLOGDIM(0) != 1) SERR("only dimension 1 allowed");
if ((err = checkstrokorbBall(cov)) != NOERROR) RETURN_ERR(err);
switch(P0INT(STROKORBBALL_DIM)) {
case 1:
if (next->rese_derivs < 2) SERR("submodel must be twice differentiable");
break;
case 3:
if (next->rese_derivs < 3)
SERR("submodel must be three times differentiable");
break;
default:
SERR("only dimensions 1 and 3 are allowed");
}
if ((err = TaylorBall(cov)) != NOERROR) RETURN_ERR(err);
RETURN_NOERROR;
}
int init_strokorbBallInner(model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
//
model *next = cov->sub[0];
int err;
if (next->randomkappa) RETURN_ERR(ERRORRANDOMKAPPA);// u.a. Taylor fehlt
assert(VDIM0 == 1 && VDIM1 == 1);
cov->mpp.maxheights[0] = 1.0;
cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1;
if (cov->mpp.moments >= 1) {
cov->mpp.mM[1] = cov->mpp.mMplus[1] = 1;
}
if ((err = TaylorBall(cov)) != NOERROR) RETURN_ERR(err);
RETURN_NOERROR;
}
void do_strokorbBallInner(model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
//model *next = cov->sub[0];
//DO(next, s); // nicht gatternr
}
void range_strokorbBallInner(model VARIABLE_IS_NOT_USED *cov, range_type *range){
range->min[STROKORBBALL_DIM] = 1;
range->max[STROKORBBALL_DIM] = 3;
range->pmin[STROKORBBALL_DIM] = 1;
range->pmax[STROKORBBALL_DIM] = 3;
range->openmin[STROKORBBALL_DIM] = false;
range->openmax[STROKORBBALL_DIM] = false;
}
void poly2unif(model *local, model *remote,
int VARIABLE_IS_NOT_USED variant) {
assert(MODELNR(local)==POLYGON && MODELNR(remote)==UNIF &&
local->Spolygon != NULL && local->Spolygon->P != NULL);
polygon *P = local->Spolygon->P;
assert(P->e != NULL);
int d,
dim = LOGDIM(SYSOF(local), 0);
assert(dim == 2);
for (d=0; d<dim; d++) {
PARAM(remote, UNIF_MIN)[d] = P->box0[d];
PARAM(remote, UNIF_MAX)[d] = P->box1[d];
}
remote->randomkappa = true;
}
void strokorbPoly(double *x, model *cov, double *v) {
// only proportional to a density !!!
model *next = cov->sub[0];
COV(x, next, v);
}
int checkstrokorbPoly(model *cov) {
model
// *key = cov->key,
*next = cov->sub[0];
int
dim = OWNLOGDIM(0),
err = NOERROR;
assert(cov->key == NULL);
if ((err= CHECK_PASSTF(next, TcfType, SCALAR, EvaluationType)) != NOERROR) RETURN_ERR(err);
// if ((err = CHECK(next, dim, cov->xdimprev, TcfType,
// OWNDOM(0), OWNISO(0),
// SCALAR, EvaluationType)) != NOERROR) RETURN_ERR(err);
if (!isGneiting(next))
SERR("member of the Gneiting-Schaback class as submodel needed");
if (dim != 2) SERR("only dimension 2 currently programmed");
if (!(hasSmithFrame(cov))) {
SERR1("'%.50s' may be used only as a shape function of a Smith field", NICK(cov));
}
setbackward(cov, next);
RETURN_NOERROR;
}
int struct_strokorbPoly(model *cov, model **newmodel) {
model *sub = cov->sub[0];
int
dim = OWNLOGDIM(0);
double var = 1.0;
model *pts=NULL, *shape=NULL;
ASSERT_NEWMODEL_NOT_NULL;
if (hasSmithFrame(cov)) {
if (SUBNR != BROWNRESNICK)
SERR1("only tcf '%.50s' allowed", DefList[BROWNRESNICK].nick);
sub = sub->sub[0];
if (isDollar(sub)) {
var = PARAM0(sub, DVAR);
sub = sub->sub[0];
}
if (SUBNR != BROWNIAN || PARAM0(sub, BROWN_ALPHA) != 1.0) {
SERR2("Numerical inverse Laplace transform has not been implemented yet. Currently, only '%.50s' with parameter %.50s=1 is a valid submodel", DefList[BROWNIAN].nick,
DefList[BROWNIAN].kappanames[BROWN_ALPHA]);
}
addModel(&pts, UNIF, NULL, true);
kdefault(pts, UNIF_NORMED, (int) false);
PARAMALLOC(pts, UNIF_MIN, dim, 1);
PARAMALLOC(pts, UNIF_MAX, dim, 1);
addModel(&shape, POLYGON, NULL, true);
addModelKappa(shape, POLYGON_BETA, ARCSQRT_DISTR);
kdefault(shape->kappasub[POLYGON_BETA], ARCSQRT_SCALE, 1.0 / var);
addSetParam(&shape, pts, poly2unif, true, 0);
addModel(newmodel, ZHOU);
kdefault(*newmodel, ZHOU_NORMED, false);
kdefault(*newmodel, ZHOU_ISOTROPIC, false);
SET_CALLING(shape, *newmodel);
SET_CALLING(pts, *newmodel);
(*newmodel)->sub[PGS_LOC] = pts;
(*newmodel)->sub[PGS_FCT] = shape;
} else ILLEGAL_FRAME_STRUCT;
RETURN_NOERROR;
}