https://github.com/cran/RandomFields
Tip revision: e994a4415e67fa60cbfd3f208aaab20872521c0b authored by Martin Schlather on 14 February 2019, 21:02:19 UTC
version 3.3
version 3.3
Tip revision: e994a44
Gneiting.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
Gneiting's space-time covariance models and related models
Copyright (C) 2006 -- 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 <Rmath.h>
#include <R_ext/Lapack.h>
#include <R_ext/Linpack.h>
#include "questions.h"
#include "operator.h"
#include "shape.h"
#include "Processes.h"
#include "startGetNset.h"
#define LOG2 M_LN2
#define AVESTP_MINEIGEN 2
#define AVESTP_LOGDET 3
#define AVESTP_V 4
#define AVESTP_LOGV 5
#define AVESTP_LOGMIXDENS 6
#define TOTALAVESTP (AVESTP_LOGMIXDENS + 1)
#define AVE_A 0
#define AVE_Z 1
#define AVE_SPACETIME 2
#define AVE_PHI 0
#define AVE_GAUSS 1
void kappa_ave(int i, model *cov, int *nr, int *nc){
bool
spacetime = (bool) (PisNULL(AVE_SPACETIME) || P0INT(AVE_SPACETIME));
int dim = OWNLOGDIM(0);
if (spacetime) dim--;
*nr = (i==AVE_A || i==AVE_Z) ? dim : 1;
*nc = (i==AVE_A) ? dim : i < DefList[COVNR].kappas ? 1 : -1;
}
void ave(double *h, model *cov, double *v) {
// f = uAu +zu;
bool
spacetime = (bool) (PisNULL(AVE_SPACETIME) || P0INT(AVE_SPACETIME));
model *next = cov->sub[0];
int i,j,k,d,
dim = OWNLOGDIM(0);
if (spacetime) dim--;
double Ah[AveMaxDim], Eplus2B[AveMaxDim],
dummy,
hh,
*A = P(AVE_A),
*z = P(AVE_Z),
c = spacetime ? h[dim] : 0.0; // Sockelwert fuer c
hh = 0.0;
for (k=d=0; d<dim; d++) {
for (dummy = 0.0, j=0; j<dim; j++) dummy += h[j] * A[k++];
Ah[d] = dummy;
c += z[d] * h[d];
hh += h[d] * h[d];
}
for (j=d=0; d<dim; d++) {
for (i=0; i<dim; i++) {
Eplus2B[j++] = 2.0 * Ah[d] * Ah[i];
}
Eplus2B[j - dim + d] += 1.0;
}
double det, AhEinvAh;
Ext_XCinvXdet(Eplus2B, dim, Ah, 1, &AhEinvAh, &det, false, NULL);
double y = SQRT(0.5 * hh + c * c * (1.0 - 2.0 * AhEinvAh));
COV(&y, next, v);
*v /= SQRT(det);
}
int checkave(model *cov) {
ASSERT_UNREDUCED;
model *next = cov->sub[0];
bool
spacetime = (bool) (PisNULL(AVE_SPACETIME) || P0INT(AVE_SPACETIME));
int i, j, err,
dim = OWNLOGDIM(0),
spdim = spacetime ? dim - 1 : dim;
double
*A = P(AVE_A);
char dim_info[2][4] = {"d", "d-1"};
if (OWNTOTALXDIM < 2) SERR("The spatial dimension must be at least 2.");
if (dim > AveMaxDim)
SERR2("For technical reasons max. dimension for ave is %d. Got %d.",
AveMaxDim, dim);
if (cov->ncol[AVE_A] != spdim || cov->nrow[AVE_A] != spdim)
SERR5("A not %.50sx%.50s matrix, but %dx%d (dim=%d)", dim_info[spacetime],
dim_info[spacetime], cov->ncol[AVE_A], cov->nrow[AVE_A], spdim);
if (cov->ncol[AVE_Z] != 1 || cov->nrow[AVE_Z] != spdim)
SERR1("z not (%.50s)-dim vector", dim_info[spacetime]);
for (i=0; i<spdim; i++)
for (j=i+1; j<spdim; j++)
if (A[i + j * spdim] != A[j + i * spdim]) {
A[j + i * spdim] = A[i + j * spdim];
warning("A is not symmetric -- lower part used");
}
// naechste Zeile stimmt nicht mit Bernoulli ueberein
kdefault(cov, AVE_SPACETIME, true);
if ((err = checkkappas(cov)) != NOERROR) RETURN_ERR(err);
if ((err = CHECK(next, dim, 1, PosDefType, XONLY, ISOTROPIC,
SCALAR, cov->frame // VariogramType changed 20.7.14 wg spectral
)) != NOERROR) RETURN_ERR(err);
if (!isNormalMixture(next->monotone)) RETURN_ERR(ERRORNORMALMIXTURE);
if (DefList[NEXTNR].spectral == NULL) RETURN_ERR(ERRORSPECTRAL); // nicht gatter
// updatepref(cov, next); ## gute idee?
if (next->pref[SpectralTBM] == PREF_NONE)
cov->pref[RandomCoin] = cov->pref[Average] = PREF_NONE;
// kein setbackward??
// no setbackard
RETURN_NOERROR;
}
void rangeave(model VARIABLE_IS_NOT_USED *cov, range_type* ra){
int i;
for (i=0; i<=1; i++) {
ra->min[i] = RF_NEGINF;
ra->max[i] = RF_INF;
ra->pmin[i] = - 10.0;
ra->pmax[i] = 10.0;
ra->openmin[i] = true;
ra->openmax[i] = true;
}
ra->min[2] = 0;
ra->max[2] = 1;
ra->pmin[2] = 0;
ra->pmax[2] = 1;
ra->openmin[2] = false;
ra->openmax[2] = false;
}
void sd_avestp(model *cov, gen_storage VARIABLE_IS_NOT_USED *S, int dim, double *sd){
///// window_info *w = &(S->window);
int d;
double b, alphamin, x2, InvSqrt2a, EmA,
*q = cov->q;
// see article/GEOSTATS/simuspacetime/simuspacetime2008/simuspacetime.tex
// for the reasoning of these calculations
BUG; // auf keinen fall parallel tauglich
q[AVESTP_LOGV] = LOG(q[AVESTP_V]);
for (x2=0.0, d=0; d<dim; d++) {
double lensimu = RF_NA; //// w->max[d] - w->min[d];
x2 += lensimu * lensimu;
}
// x2 *= 0.25;
b = 3.0 * q[AVESTP_V] * x2 / dim;
alphamin = (4.0 + 4.0 * b - 2.0 * SQRT(4.0 * b * b + 8.0 * b + 1.0)) / 3.0;
InvSqrt2a = 1.0 / SQRT(2.0 * alphamin * 6.0 * q[AVESTP_V]);
*sd = InvSqrt2a;
EmA = 1.0 - alphamin;
cov->mpp.maxheights[0] = EXP(-0.5 * LOG(EmA) - 0.25 * LOG(alphamin) + b / EmA -
2 * x2); // proportional zum dritten Moment !
/*
double radius =
SQRT((-9 // so e^{-9} as threshold
- 0.25 * dim * (q[AVESTP_LOGV] - 1.14473) // log pi
- 0.25 * q[AVESTP_LOGDET]
//+ 0.5 * cov_a->logdens
- q[AVESTP_LOGMIXDENS]
) / ( - q[AVESTP_MINEIGEN] * q[AVESTP_V]) ); // ???
assert(radius > 0);
*/
}
int structAve(model *cov, model **newmodel) {
model *shape; //, *gaussmix=NULL;
int err;
ASSERT_NEWMODEL_NOT_NULL;
ERR("'ave' currently does not work"); // what to do with the next lines??
if ((err = covcpy(newmodel, cov)) != NOERROR) RETURN_ERR(err);
shape = *newmodel;
SET_NR(shape, SHAPEAVE);
addModel(shape, AVE_GAUSS, GAUSS);
ERR("'ave' currently does not work");
// what to do ERR1("'ave' currently does not work (% ld)", (long int) gaussmix); // what to do with the next lines??
/*
gaussmix = shape->sub[AVE_GAUSS];
gaussmix->ts dim = 1;
gaussmix->frame = GaussMethodType;
gaussmix->method = SpectralTBM;
*/
RETURN_NOERROR;
}
void logshapeave(double *x, model *cov, double *v, double *Sign) {
// nur stationaer
bool
spacetime = (bool) (PisNULL(AVE_SPACETIME) || P0INT(AVE_SPACETIME));
int d, j, k,
dim = OWNLOGDIM(0);
if (spacetime) dim--;
double f, dummy, r2,
*A = P(AVE_A),
*z = P(AVE_Z),
t = spacetime ? x[dim] : 0.0,
*q = cov->q;
f = r2 = 0.0;
for (k=d=0; d<dim; d++) {
r2 += x[d] * x[d];
for (dummy = z[d], j=0; j<dim; j++) dummy += x[j] * A[k++];
f += dummy * x[d];
}
#ifdef SCHLATHERS_MACHINE
warning("is exponent of V correct?");
#endif
v[0] = 0.25 * dim * q[AVESTP_LOGV] // V^{d/4 oder d/2} !!!!!!!!!!!
- 0.5 * (LOG2 - dim * M_LN_SQRT_PId2) - r2
// LOG2 : wegen spectral simulation;
// zweiter Term fuer logg
//+ DefList[phi->nr].logmixdens(x, q[AVESTP_LOGV], phi); /* g */// nicht gatternr
;
Sign[0] = 1.0;
double phase = q[AVERAGE_YPHASE] + q[AVERAGE_YFREQ] * (f - t); // Y
Sign[1] = phase > 0.0 ? 1.0 : phase < 0.0 ? -1.0 : 0.0;
v[1] = LOG(FABS(phase));
}
int check_shapeave(model *cov) {
if (cov->sub[AVE_GAUSS] == NULL)
SERR1("both submodels must be set to '%.50s'", DefList[GAUSS].nick);
cov->mpp.maxheights[0] = RF_NA;
return checkave(cov); // !! not next
}
int init_shapeave(model *cov, gen_storage *s) {
ASSERT_GAUSS_METHOD(Average);
model
*gaussmix = cov->sub[AVE_GAUSS];
double sd,
*q = cov->q;
bool
spacetime = (bool) (PisNULL(AVE_SPACETIME) || P0INT(AVE_SPACETIME));
int err = NOERROR,
dim = OWNLOGDIM(0);
if (spacetime) dim--;
q[AVESTP_V] = 0.0;
q[AVESTP_MINEIGEN] = 1.0;
q[AVESTP_LOGDET] = 0.0;
sd_avestp(cov, s, dim, &sd); // sd->gauss
assert(VDIM0 == 1); assert(VDIM0 == VDIM1);
if (cov->mpp.moments >= 0) {
cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0;
if (cov->mpp.moments >= 1) {
if ((err = INIT(gaussmix, cov->mpp.moments, s)) != NOERROR) RETURN_ERR(err);
if (cov->mpp.moments >= 2) {
cov->mpp.mM[2] = 1.0;
}
}
}
RETURN_ERR(err);
}
void do_shapeave(model *cov, gen_storage *S) {
// Simulation of V; sopee Bernoulli Sec. 4.2
model
*aveGAUSS = cov->sub[AVE_GAUSS],
*phi = cov->sub[AVE_PHI];
double spec_ret[StpMaxDim], sd,
*q = cov->q;
bool
spacetime = (bool) (PisNULL(AVE_SPACETIME) || P0INT(AVE_SPACETIME));
int
dim = OWNLOGDIM(0);
if (spacetime) dim--;
DefList[MODELNR(phi)].drawmix(phi, q + AVESTP_V); // nicht gatternr
sd_avestp(cov, S, dim, &sd); // sd->gauss
BUG;
SPECTRAL(aveGAUSS, S, spec_ret); // nicht gatternr
q[AVERAGE_YFREQ] = *spec_ret * q[AVESTP_V];
q[AVERAGE_YPHASE] = TWOPI * UNIFORM_RANDOM;
BUG; // what to do with the next line?
}
/* coxgauss, cmp with nsst1 !! */
// C = 2 (C + 4 M H M), H = h h^t
// a = t - h M h - zh
// EXP(- 0.5 * (h *h + 2 a^2 - mu C mu)) // stimmen die Vorzeichen??
// mu = h - 2 a M h
/* cox, cmp with nsst1 !! */
// coxisham
#define COX_MU 0
#define COX_D 1
#define COX_BETA 2
void GetEu2Dinv(model *cov, double *x, int dim,
double *det, double *Eu2Dinv,
double *newxsq, double *newx, double *z) {
double
y[CoxMaxDim],
*V = P(COX_MU),
*D= P(COX_D),
beta = P0(COX_BETA),
t = x[dim],
t2 = POW(FABS(t), beta); // standard t^2
int d,
dimP1 = dim + 1,
dimsq = dim * dim;
for (d=0; d<dim; d++) {
y[d] = x[d] - t * V[d];
}
for (d=0; d<dimsq; d++) {
Eu2Dinv[d] = t2 * D[d];
}
for (d=0; d<dimsq; d+=dimP1) Eu2Dinv[d] += 1.0; // D + E
solve_storage *PT = NULL;
if (z !=NULL) PT->result = z;
int err = Ext_XCinvXdet(Eu2Dinv, dim, y, 1, newxsq, det, false, PT);
if (z != NULL) PT->result = NULL;
if (err != NOERROR) ERR("error occuredin 'GetEu2Dinv'");
*newx = SQRT(*newxsq);
}
void cpyUf(double *Eu2Dinv, double factor, int dim, int logicaldim, double *v) {
// Eu2Dinv has dimension dim^2; v dimension logicaldim^2
// Eu2Dinv is copied into the upper left part of v and
// multiplied by factor
int d, i, k, j,
logicaldimsq = logicaldim * logicaldim;
for (i=0; i<logicaldimsq; v[i++] = 0.0);
for (i=0; i<dim; i++) {
for (d=i * logicaldim, k=i * dim, j=0; j<=i; j++)
v[d++] = Eu2Dinv[k++] * factor;
for (k += dim - 1; j<dim; j++, k+=dim) {
v[d++] = Eu2Dinv[k] * factor;
}
}
}
void addzzT(double *v, double factor, double *z, int dim, int logicaldim) {
// z has dimension dim; v dimension logicaldim^2
// zzT is copied into the upper left part of v after being
// multiplied by factor
int i,j,k;
for (i=0; i<dim; i++) {
k = i * logicaldim;
for (j=0; j<dim; j++) {
v[k++] += z[i] * z[j] * factor;
}
}
}
void kappa_cox(int i, model *cov, int *nr, int *nc){
switch (i) {
case COX_MU :
*nc = 1;
*nr = OWNLOGDIM(0) - 1;
break;
case COX_D :
*nc = *nr = OWNLOGDIM(0) - 1;
break;
case COX_BETA :
*nc = *nr = 1;
break;
default: *nc = *nr = OUT_OF_RANGE;
}
}
void cox(double *x, model *cov, double *v) {
model *next = cov->sub[0];
int dim = OWNLOGDIM(0) - 1,
dimsq = dim * dim;
double det, newx, newxsq;
TALLOC_XX1(Eu2Dinv, dimsq);
GetEu2Dinv(cov, x, dim, &det, Eu2Dinv, &newxsq, &newx, NULL);
END_TALLOC_XX1;
// printf("cox det=%10g %10g\n", det, newx);
COV(&newx, next, v);
*v /= SQRT(det);
}
void coxhess(double *x, model *cov, double *v) {
model *next = cov->sub[0];
int logicaldim = OWNLOGDIM(0),
dim = logicaldim - 1,
dimsq = dim * dim;
double z[CoxMaxDim], det, newx, newxsq, phiD, phiD2;
TALLOC_XX1(Eu2Dinv, dimsq);
GetEu2Dinv(cov, x, dim, &det, Eu2Dinv, &newxsq, &newx, z);
Abl2(&newx, next, &phiD2);
if (newxsq == 0.0) {
cpyUf(Eu2Dinv, phiD2 / SQRT(det), dim, logicaldim, v);
} else {
Abl1(&newx, next, &phiD);
cpyUf(Eu2Dinv, phiD / (SQRT(det) * newx), dim, logicaldim, v);
addzzT(v, (phiD2 - phiD/newx) / (SQRT(det) * newxsq), z, dim, logicaldim);
}
END_TALLOC_XX1;
}
void coxnabla(double *x, model *cov, double *v) {
model *next = cov->sub[0];
int d,
logicaldim = OWNLOGDIM(0),
dim = logicaldim - 1,
dimsq=dim * dim;
double z[CoxMaxDim], det, newx, newxsq, phiD, factor;
TALLOC_XX1(Eu2Dinv, dimsq);
GetEu2Dinv(cov, x, dim, &det, Eu2Dinv, &newxsq, &newx, z);
END_TALLOC_XX1;
if (newxsq == 0.0) {
for (d=0; d<=dim; d++) v[d] = 0.0;
} else {
newx = SQRT(newxsq);
Abl1(&newx, next, &phiD);
factor = phiD / (det * newx);
for (d=0; d<dim; d++) {
v[d] = factor * z[d];
}
for (d=0; d<logicaldim; v[d++]=0.0);
}
}
int checkcox(model *cov) {
model *next = cov->sub[0];
int err, i,
logicaldim = OWNLOGDIM(0),
dim = logicaldim - 1,
dimsq = dim * dim;
if (OWNTOTALXDIM < 2) SERR("The space-time dimension must be at least 2.");
if (cov->ncol[COX_MU] != 1 || cov->nrow[COX_MU] != dim) {
if (cov->ncol[COX_MU] == dim && cov->nrow[COX_MU] == 1) {
cov->nrow[COX_MU] = dim;
cov->ncol[COX_MU] = 1;
} else {
SERR3("mu is not given or not a vector of dimen. %d (nrow=%d ncol=%d)",
dim, cov->nrow[COX_MU], cov->ncol[COX_MU]);
}
}
// is matrix positive definite?
if (PisNULL(COX_D)) {
PALLOC(COX_D, dim, dim);
for (i=0; i<dimsq; i++) P(COX_D)[i] = 1.0;
} else {
if (!Ext_is_positive_definite(P(COX_D), dim))
SERR("D is not (strictly) positive definite");
}
kdefault(cov, COX_BETA, 2.0);
if ((err = checkkappas(cov)) != NOERROR) RETURN_ERR(err);
if ((err = CHECK(next, dim, 1, PosDefType, XONLY, ISOTROPIC,
SCALAR, cov->frame // VariogramType changed 20.7.14 wg spectral
)) != NOERROR)
RETURN_ERR(err);
if (logicaldim != 3) cov->pref[SpectralTBM] = PREF_NONE;
if (!isNormalMixture(next->monotone)) RETURN_ERR(ERRORNORMALMIXTURE);
if (DefList[NEXTNR].spectral == NULL) RETURN_ERR(ERRORSPECTRAL); // nicht gatternr
// no setbackard
updatepref(cov, next);
if (P0(COX_BETA) != 2.0) cov->pref[SpectralTBM] = 0;
cov->hess = true;
if (cov->Ssolve == NULL) {EXT_NEW_STORAGE(solve);}
EXTRA_STORAGE;
RETURN_NOERROR;
}
void rangecox(model VARIABLE_IS_NOT_USED *cov, range_type* range){
range->min[COX_MU] = RF_NEGINF;
range->max[COX_MU] = RF_INF;
range->pmin[COX_MU] = - 100.0;
range->pmax[COX_MU] = +100.0;
range->openmin[COX_MU] = true;
range->openmax[COX_MU] = true;
range->min[COX_D] = RF_NEGINF;
range->max[COX_D] = RF_INF;
range->pmin[COX_D] = - 100.0;
range->pmax[COX_D] = +100.0;
range->openmin[COX_D] = false;
range->openmax[COX_D] = false;
range->min[COX_BETA] = 0.0;
range->max[COX_BETA] = 2.0;
range->pmin[COX_BETA] = 0.1;
range->pmax[COX_BETA] = 2.0;
range->openmin[COX_BETA] = true;
range->openmax[COX_BETA] = false;
}
int initcox(model *cov, gen_storage *s) {
ASSERT_GAUSS_METHOD(SpectralTBM);
model *next = cov->sub[0];
return INIT(next, 0, s);
}
void spectralcox(model *cov, gen_storage *s, double *e) {
model *next = cov->sub[0];
int d,
logicaldim = OWNLOGDIM(0),
dim = logicaldim - 1;
double t, v[CoxMaxDim],
*V = P(COX_MU),
rho= P0(COX_D);
SPECTRAL(next, s, e); // nicht gatternr
v[0] = rnorm(0.0, INVSQRTTWO);
v[1] = rho * v[0] + SQRT(1.0 - rho * rho) * rnorm(0.0, INVSQRTTWO);
for (t = 0.0, d=0; d<dim; d++) {
t += (v[d] + V[d]) * e[d];
}
e[dim] = -t;
}
// GenPS (generalisation of paciore and stein)
// Q = (x-y) Sy (Sx + Sy)^{-1} Sx (x-y) (weicht etwas von Stein ab)
#define STP_XI2 0
#define STP_PHI 1
#define STP_GAUSS 3
void kappa_stp(int i, model *cov, int *nr, int *nc){
*nc = (i == STP_S || i == STP_M) ? OWNLOGDIM(0) : 1;
*nr = i < DefList[COVNR].kappas ? OWNLOGDIM(0) : -1;
}
void stp(double *x, double *y, model *cov, double *v) {
int d, j, k, err,
dim = OWNLOGDIM(0),
dimsq = dim * dim;
double h[StpMaxDim],
Mh[StpMaxDim], hSx[StpMaxDim],
Syh[StpMaxDim], xi2x, xi2y,
detA, hMh, cxy, zh, Q, Amux[StpMaxDim], Amuy[StpMaxDim],
*Sc = P(STP_S),
*M = P(STP_M),
*z = P(STP_Z);
model *phi = cov->sub[STP_PHI],
*Sf = cov->kappasub[STP_S],
*xi2 =cov->sub[STP_XI2];
TALLOC_X1(Sx, dimsq);
TALLOC_X2(Sy, dimsq);
TALLOC_X3(A, dimsq);
if (Sf != NULL) {
FCTN(x, Sf, Sx); // symmetric, pos definite !!
FCTN(y, Sf, Sy);
//
// if (false) {
// int ii;
// printf("x=%10g %10g y=%10g %10g\n", x[0], x[1], y[0], y[1]);
//for (ii=0; ii<4; ii++) printf("%10g ", Sx[ii]);
// printf("\n");
// for (ii=0; ii<4; ii++) printf("%10g ", Sy[ii]);
// printf("\n");
// }
} else {
int bytes = sizeof(double) * dimsq;
MEMCOPY(Sx, Sc, bytes);
MEMCOPY(Sy, Sc, bytes);
}
if (xi2 != NULL) {
FCTN(x, xi2, &xi2x);
FCTN(y, xi2, &xi2y);
} else {
xi2x = xi2y = 0.0;
}
for (k=0, d=0; d<dim; d++) {
h[d] = x[d] - y[d];
}
zh = hMh = 0.0;
for (k=0, d=0; d<dim; d++) {
Mh[d] = hSx[d] = Syh[d] = 0.0;
for (j=0; j<dim; j++, k++) {
Mh[d] += h[j] * M[k];
hSx[d] += h[j] * Sx[k];
Syh[d] += h[j] * Sy[k]; // uses that S is symmetric
}
zh += z[d] * h[d];
hMh += Mh[d] * h[d];
}
cxy = xi2x - xi2y - zh;
for (k=d=0; d<dim; d++) {
for (j=0; j<dim; j++, k++) {
A[k] = Sx[k] + Sy[k] + 4.0 * Mh[d] * Mh[j];
assert(R_FINITE(A[k]));
}
Amux[d] = hSx[d] + 2.0 * (hMh + cxy) * Mh[d]; // uses that M is symmetric
Amuy[d] = Syh[d] + 2.0 * (hMh - cxy) * Mh[d];
}
double
dx = Ext_detPosDef(Sx, dim),
dy = Ext_detPosDef(Sy, dim);
END_TALLOC_X1;
END_TALLOC_X2;
// det_UpperInv(A, &detA, dim); // here
double muxAmuy; // changed 5.4.2018, 3.2.0.28
err = Ext_XCinvYdet(A, dim, true, Amux, Amuy, 1, &muxAmuy, &detA, false,NULL);
END_TALLOC_X3;
if (err != NOERROR) ERR("error occuredin 'GetEu2Dinv'");
Q = cxy * cxy - hMh * hMh + muxAmuy;
if (Q < 0.0) {
PRINTF("x=%10g,%10g y=%10g,%10g detA=%10g\n",
x[0], x[1], y[0], y[1], detA);
PRINTF("cxy=%4f hMh=%10g Amux=%10g Amuy=%10g\ndim=%d h=%10g,%10g hSx=%10g,%10g, Q=%10g\n",
cxy, hMh, Amux[0], Amuy[0],
dim, h[0], h[1], hSx[0], hSx[1], Q);
BUG;
}
Q = SQRT(Q);
aux_covfct auxcf;
if ((auxcf = DefList[MODELNR(phi)].aux_cov) != NULL) {
BUG;
auxcf(x, y, Q, phi, v);
} else
FCTN(&Q, phi, v);
*v *= POW(2.0, 0.5 * double(dim)) * POW(dx * dy / (detA * detA), 0.25);
}
int checkstp(model *cov){
ASSERT_UNREDUCED;
model
*phi = cov->sub[STP_PHI],
*Sf = cov->kappasub[STP_S],
*xi2 =cov->sub[STP_XI2];
int err,
dim = OWNLOGDIM(0);
ASSERT_CARTESIAN;
if (dim > StpMaxDim)
SERR2("For technical reasons max. dimension for ave is %d. Requested is %d",
StpMaxDim, GATTERXDIM(0));
if (PisNULL(STP_S) && Sf==NULL) { // Sc
if ((cov->px[STP_S] = EinheitsMatrix(dim)) == NULL)
RETURN_ERR(ERRORMEMORYALLOCATION);
cov->ncol[STP_S] = cov->nrow[STP_S] = dim;
}
if (PisNULL(STP_M)) { // M
if ((cov->px[STP_M] = EinheitsMatrix(dim)) == NULL)
RETURN_ERR(ERRORMEMORYALLOCATION);
cov->ncol[STP_M] = cov->nrow[STP_M] = dim;
}
if (PisNULL(STP_Z)) { // z
PALLOC(STP_Z, dim, 1);
}
if ((err = CHECK(phi, dim, 1, PosDefType, XONLY, ISOTROPIC,
SCALAR, cov->frame // VariogramType changed 20.7.14 wg spectral
)) != NOERROR)
RETURN_ERR(err);
if (!isNormalMixture(phi->monotone)) RETURN_ERR(ERRORNORMALMIXTURE);
cov->pref[Average] = 5;
if (Sf != NULL) {
if ((err = CHECK(Sf, dim, dim, ShapeType, XONLY, CARTESIAN_COORD,
dim, cov->frame // VariogramType changed 20.7.14 wg spectral
)) != NOERROR)
RETURN_ERR(err);
}
if (xi2 != NULL) {
if ((err = CHECK(xi2, dim, dim, ShapeType, XONLY, CARTESIAN_COORD,
SCALAR, cov->frame // VariogramType changed 20.7.14 wg spectral
)) != NOERROR)
RETURN_ERR(err);
}
// kein setbackward??
EXTRA_STORAGE;
cov->mpp.maxheights[0] = RF_NA;
RETURN_NOERROR;
}
void rangestp(model VARIABLE_IS_NOT_USED *cov, range_type* range){
int i;
for (i=0; i<=2; i++) { /* S, M, z */
range->min[i] = RF_NEGINF;
range->max[i]= RF_INF;
range->pmin[i] = - 1e10;
range->pmax[i] = 1e10;
range->openmin[i] = true;
range->openmax[i] = true;
}
}
int structStp(model *cov, model **newmodel) {
model *shape;
int err;
ASSERT_NEWMODEL_NOT_NULL;
if ((err = covcpy(newmodel, cov)) != NOERROR) RETURN_ERR(err);
shape = *newmodel;
SET_NR(shape, SHAPESTP);
addModel(shape, STP_GAUSS, GAUSS);
ERR("'stp' currently does not work"); // what to do with the next line??
/*
shape->sub[STP_GAUSS]->ts dim = 1;
*/
RETURN_NOERROR;
}
void logshapestp(double *x, double *u, model *cov, double *v, double *Sign){
// kann um ca. Faktor 2 beschleunigt werden, wenn
// Sx , log det U, Hx fuer alle x abgespeichert werden
// und die Werte dann aus dem Speicher gelesen werden
// jedoch sehr Speicherintensiv. MEMCOPY braucht man auch nicht
model
*Sf = cov->kappasub[STP_S],
*xi2 =cov->sub[STP_XI2];
int j, k, d,
// dim= cov->xdimprev,
dim= OWNXDIM(0),
dimsq = dim * dim,
bytes = sizeof(double) * dimsq;
double h[StpMaxDim], hSxh, hSx, xi, Mhd,
*Sc = P(STP_S),
*M = P(STP_M),
*z = P(STP_Z),
*q = cov->q;
assert(q != NULL);
TALLOC_X1(Sx, dimsq);
if (Sf == NULL) {
MEMCOPY(Sx, Sc, bytes);
} else {
FCTN(x, Sf, Sx); // symmetric, pos definite !!
}
if (xi2 == NULL) {
xi = 0.0;
} else {
FCTN(x, xi2, &xi);
}
for (k=0, d=0; d<dim; d++) {
h[d] = u[d] - x[d];
}
hSxh = 0.0;
for (k=0, d=0; d<dim; d++) {
Mhd = hSx = 0.0;
for (j=0; j<dim; j++) {
Mhd += h[j] * M[k];
hSx += h[j] * Sx[k++];
}
xi += Mhd * h[d] + z[d] * h[d];
hSxh += hSx * h[d];
}
double exponent =
0.25 * dim * (// M_LN2 + ??? !!! Rechnung!!!
q[AVESTP_LOGV] - 2.0 * M_LN_SQRT_PI) // (2V/pi)^{d/4}
+ 0.25 * LOG(Ext_detPosDef(Sx, dim)) // Sx ^1/4
- q[AVESTP_V] * hSxh // EXP(-V(U-x) S (U-x))
// + DefList[phi->nr].logmixdens(x, q[AVESTP_LOGV], phi) // g //nicht gatternr
// - 0.5 * cov_a->logdens // f
; // 1 / SQRT(f)
if (!(exponent < 5.0) && PL >= PL_DETAILS) {
if (!(exponent < 6.0)) // could be NA, too
PRINTF("\n%10g log Det U=%10g %10g expon=%10g",
0.25 * dim * (// M_LN2 + ??? !!! Rechnung!!!
q[AVESTP_LOGV] - 2.0 * M_LN_SQRT_PI) // (2V/pi)^{d/4}
, 0.25 * LOG(Ext_detPosDef(Sx, dim)) /// Sx ^1/4
, -q[AVESTP_V]* hSxh // EXP(-V(U-x) S (U-x))
// , DefList[phi->nr].logmixdens(x, q[AVESTP_LOGV], phi)// g
//, - 0.5 * cov_a->logdens // f
, exponent);
else PRINTF("!");
};
assert(EXP(exponent) < 10000000.0);
double cos_value = COS(q[AVERAGE_YPHASE] + q[AVERAGE_YFREQ] * xi);
*v = exponent + LOG(FABS(cos_value)) ; // Y
*Sign = cos_value > 0.0 ? 1.0 : cos_value < 0.0 ? -1.0 : 0.0;
END_TALLOC_X1;
}
int check_shapestp(model *cov) {
if (cov->sub[AVE_GAUSS] == NULL)
SERR1("both submodels must be set to '%.50s'", DefList[GAUSS].nick);
EXTRA_STORAGE;
return checkstp(cov); // !! not next
}
int init_shapestp(model *cov, gen_storage *s) {
ASSERT_GAUSS_METHOD(Average);
model
*Sf = cov->kappasub[STP_S],
*gaussmix = cov->sub[STP_GAUSS];
double sd,
*q = cov->q;
int err = NOERROR;
if (Sf != NULL) {
double minmax[2];
assert(DefList[MODELNR(Sf)].minmaxeigenvalue != NULL);
DefList[MODELNR(Sf)].minmaxeigenvalue(Sf, minmax);
if (minmax[0] <= 0.0) ERR("neg eigenvalue in shape function of 'stp'");
q[AVESTP_MINEIGEN] = minmax[0];
q[AVESTP_LOGDET] = (double) //cov->xdimprev
OWNXDIM(0) * LOG(minmax[1]);
} else {
#define dummyN (5 * StpMaxDim)
double value[StpMaxDim], ivalue[StpMaxDim], dummy[dummyN], det,
min = RF_INF;
int i, Ferr,
dim = OWNLOGDIM(0),
ndummy = dummyN;
F77_NAME(dgeev)("No", "No", &dim, P(STP_S), // SVD
&dim, value, ivalue, NULL, &dim, NULL, &dim,
dummy, &ndummy, &Ferr);
if (Ferr != 0) RETURN_ERR(ERRORSVD);
det = 1.0;
for (i=0; i<dim; i++) {
double v = FABS(value[i]);
det *= v;
if (min > v) min = v;
}
q[AVESTP_MINEIGEN] = min;
q[AVESTP_LOGDET] = LOG(det);
}
q[AVESTP_LOGV] = 0.0;
q[AVESTP_LOGMIXDENS] = 0.0;
sd_avestp(cov, s, OWNLOGDIM(0), &sd); // sd->gauss
assert(VDIM0 == 1); assert(VDIM0 == VDIM1);
if (cov->mpp.moments >= 0) {
cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; //// ??? notwendig
if (cov->mpp.moments >= 1) {
if ((err = INIT(gaussmix, 2, s)) != NOERROR) RETURN_ERR(err);
if (cov->mpp.moments >= 2) cov->mpp.mM[2] = 1.0;
}
}
RETURN_ERR(err);
}
void do_shapestp(model *cov, gen_storage *s) {
// Simulation of V; see Bernoulli Sec. 4.2
model
*stpGAUSS = cov->sub[STP_GAUSS],
*phi = cov->sub[STP_PHI];
double spec_ret[StpMaxDim], sd,
*q = cov->q;
DefList[MODELNR(phi)].drawmix(phi, &(q[AVESTP_V]));
sd_avestp(cov, s, OWNLOGDIM(0), &sd); // sd->gauss
BUG;
SPECTRAL(stpGAUSS, s, spec_ret); // nicht gatternr
q[AVERAGE_YFREQ] = *spec_ret * SQRT(q[AVESTP_V]);
q[AVERAGE_YPHASE] = TWOPI * UNIFORM_RANDOM;
BUG; /// what to do with the next line?
// info->logdens = DefList[phi->nr].logmixdens(ZERO, q[AVESTP_LOGV], phi);
//info->radius = RF_INF;
// info-sd s.o.
}
/* nsst */
/* Tilmann Gneiting's space time models, part I */
void nsst(double *x, model *cov, double *v) {
model *subphi = cov->sub[0];
model *subpsi = cov->sub[1];
double v1, v2, psi, y;
COV(ZERO(subpsi), subpsi, &v1);
COV(x + 1, subpsi, &v2);
psi = SQRT(1.0 + v1 - v2); // C0 : C(0) oder 0 // Cx : C(x) oder -gamma(x)
y = x[0] / psi;
COV(&y, subphi, v);
*v *= POW(psi, -P0(NSST_DELTA));
}
void Dnsst(double *x, model *cov, double *v) {
model *subphi = cov->sub[0];
model *subpsi = cov->sub[1];
double v1, v2, psi, y;
COV(ZERO(subpsi), subpsi, &v1);
COV(x + 1, subpsi, &v2);
psi = SQRT(1.0 + v1 - v2); // C0 : C(0) oder 0 // Cx : C(x) oder -gamma(x)
y = x[0] / psi;
Abl1(&y, subphi, v);
*v *= POW(psi, -P0(NSST_DELTA) - 1.0);
}
void TBM2nsst(double *x, model *cov, double *v) {
model *subphi = cov->sub[0];
model *subpsi = cov->sub[1];
double v1, v2, psi, y;
assert(false);
COV(ZERO(subpsi), subpsi, &v1);
COV(x + 1, subpsi, &v2);
psi = SQRT(1.0 + v1 - v2); // C0 : C(0) oder 0 // Cx : C(x) oder -gamma(x)
y = x[0] / psi;
TBM2CALL(&y, subphi, v);
*v *= POW(psi, -P0(NSST_DELTA));
}
int checknsst(model *cov) {
model *subphi = cov->sub[0];
model *subpsi = cov->sub[1];
int err;
if (OWNXDIM(0) != 2) SERR("reduced dimension must be 2");
if ((err = checkkappas(cov)) != NOERROR) RETURN_ERR(err);
cov->finiterange = falsch;
if ((err = CHECK(subphi, OWNLOGDIM(0), 1, PosDefType, XONLY, ISOTROPIC,
SCALAR, cov->frame // VariogramType changed 20.7.14 wg spectral
)) != NOERROR)
RETURN_ERR(err);
if (!isNormalMixture(subphi->monotone)) return(ERRORNORMALMIXTURE);
setbackward(cov, subphi);
assert(cov->finiterange != wahr); //
if ((err = CHECK(subpsi, 1, 1, VariogramType, XONLY, ISOTROPIC,
SCALAR, cov->frame // VariogramType changed 20.7.14 wg spectral
)) != NOERROR)
RETURN_ERR(err);
// kein setbackward(cov, subpsi);
RETURN_NOERROR;
}
void rangensst(model *cov, range_type* range){
range->min[NSST_DELTA] = OWNLOGDIM(0) - 1;
range->max[NSST_DELTA] = RF_INF;
range->pmin[NSST_DELTA] = range->min[NSST_DELTA];
range->pmax[NSST_DELTA] = 10.0;
range->openmin[NSST_DELTA] = false;
range->openmax[NSST_DELTA] = true;
}
/* gennsst */
/* generalisation in schlather, bernoulli 2010 */
#define GENNSST_INTERN_A 0
void kappa_gennsst_intern(int i, model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc){
*nc = *nr = i == 0 ? 0 : -1;
}
void gennsst_intern(double *x, model *cov, double *v) {
model *next = cov->sub[0];
double det, z,
*A = P(GENNSST_INTERN_A);
int dim = OWNLOGDIM(0);
// printf("\nx=%10g %10g %10g \ny=%10g %10g %10g\nz=%10g %10g %10g\n %10g %10g %10g\n %10g %10g %10g\n", x[0], x[1], x[2], y[0],y[1],y[2], A[0], A[1], A[2], A[3], A[4], A[5], A[6], A[7], A[8]);
// PMI(cov->calling->sub[1]);
// det_UpperInv(A, &det, dim);
if (Ext_XCinvXdet(A, dim, x, 1, &z, &det, false, NULL) != NOERROR) {
*v = RF_NAN;
return;
}
z = SQRT(z);
COV(&z, next, v);
*v /= SQRT(det);
}
int checkgennsst_intern(model *cov) {
// internes modell muss zwischengeschaltet werden
// um Earth->Cartesian->isotropic hinzubekommen
// und ^- hier noch mit A zu multiplizieren
model *next = cov->sub[0];
int err,
dim = OWNXDIM(0);
assert(OWNLOGDIM(0) == OWNXDIM(0));
if ((err = CHECK(next, OWNLOGDIM(0), 1, PosDefType,
XONLY, ISOTROPIC, SCALAR, cov->frame)) != NOERROR)
RETURN_ERR(err);
if (!isNormalMixture(next->monotone)) RETURN_ERR(ERRORNORMALMIXTURE);
if (PisNULL(GENNSST_INTERN_A)) {
PALLOC(GENNSST_INTERN_A, dim, dim);
} else if (cov->nrow[GENNSST_INTERN_A] != dim) {
PFREE(GENNSST_INTERN_A);
PALLOC(GENNSST_INTERN_A, dim, dim);
}
cov->finiterange = falsch;
setbackward(cov, next);
assert(cov->finiterange != wahr);
VDIM0 = VDIM1 = 1;
EXTRA_STORAGE;
RETURN_NOERROR;
}
void range_gennsst_intern(model VARIABLE_IS_NOT_USED *cov,
range_type* range){
range->min[GENNSST_INTERN_A] = RF_NEGINF;
range->max[GENNSST_INTERN_A] = RF_INF;
range->pmin[GENNSST_INTERN_A] = -100.0;
range->pmax[GENNSST_INTERN_A] = +100.0;
range->openmin[GENNSST_INTERN_A] = true;
range->openmax[GENNSST_INTERN_A] = true;
}
#define GENNSST_DIM_U 0
void gennsst(double *x, model *cov, double *v) {
model *subphi = cov->key;
model *subpsi = cov->sub[1];
int
totaldim = ANYDIM,
udim = P0INT(GENNSST_DIM_U),
hdim = totaldim - udim,
hdimSq = hdim * hdim;
COV(x + hdim, subpsi, PARAM(subphi, GENNSST_INTERN_A));
if (isnowVariogram(subpsi)) {
TALLOC_X1(z, hdimSq);
double *zero = ZERO(subpsi),
*p = PARAM(subphi, GENNSST_INTERN_A);
COV(zero, subpsi, z);
for (int i=0; i<hdimSq; i++) p[i] = z[i] - p[i];
END_TALLOC_X1;
} else if (!equalsnowNegDef(subpsi)) BUG;
COV(x, subphi, v);
if (ISNAN(*v)) ERR("error occuredin 'GetEu2Dinv'");
}
void nonstatgennsst(double *x,double *y, model *cov, double *v) {
model *subphi = cov->key;
model *subpsi = cov->sub[1];
int
totaldim = ANYDIM,
udim = P0INT(GENNSST_DIM_U),
hdim = totaldim - udim,
hdimSq = hdim * hdim;
NONSTATCOV(x + hdim, y + hdim, subpsi, PARAM(subphi, GENNSST_INTERN_A));
TALLOC_X1(z, hdimSq);
if (isnowVariogram(subpsi)) {
double *zero = ZERO(subpsi),
*p = PARAM(subphi, GENNSST_INTERN_A);
NONSTATCOV(zero, zero, subpsi, z);
for (int i=0; i<hdimSq; i++) p[i] = z[i] - p[i];
} else if (!equalsnowNegDef(subpsi)) BUG;
for (int d=0; d<hdim; d++) z[d] = x[d] - y[d];
COV(z, subphi, v);
END_TALLOC_X1;
if (ISNAN(*v)) ERR("error occuredin 'GetEu2Dinv'");
}
#define GENNSST_DIM_U 0
int checkgennsst(model *cov) {
model *subphi = cov->sub[0];
model *subpsi = cov->sub[1];
#define ncs 3
int err;
kdefault(cov, GENNSST_DIM_U, 1);
if (OWNLOGDIM(0)!=OWNXDIM(0)) SERR("logical and physical dimension differ");
int udim = P0INT(GENNSST_DIM_U),
hdim = OWNLOGDIM(0) - udim;
if (cov->q == NULL) {
QALLOC(1);
cov->q[0]=NOERROR;
}
if (isSpherical(OWNISO(0)))
return cov->q[0] != NOERROR ? (int) cov->q[0] : ERRORFAILED;
if (cov->key == NULL) {
if ((err = covcpy(&(cov->key), subphi)) != NOERROR) RETURN_ERR(err);
addModel(&(cov->key), GENNSST_INTERN);
}
if ((cov->q[0] = err = CHECK(cov->key, hdim, hdim, PosDefType,
XONLY, SYMMETRIC, // viel zu schwach, geht aber
// z.Zt. nicht anders
// hinsichtlich der Koordinatentransformation
// -- wird das interne modell abgefangen
SCALAR, cov->frame)) != NOERROR) {
RETURN_ERR(err);
}
if (!isNormalMixture(cov->key->sub[0])) SERR("'phi' not a normal mixture.");
if (TOTALXDIM(SYSOF(cov->key)) != hdim)
SERR("given dim does not match dimension required for 'phi'");
// MUSSS ZWINGEND ALS ZWEITES KOMMEN
domain_type type = XONLY;
while (true) {
if ((err = CHECK(subpsi, udim, udim, NegDefType, type,
OWNISO(0), hdim, cov->frame)) == NOERROR) break;
if (type == XONLY) type = KERNEL; else RETURN_ERR(err);
}
if ((!equalsSpaceIsotropic(OWNISO(0)) ||
MODELNR(subpsi) != MATRIX || PisNULL(M_M) ||
subpsi->kappasub[M_M] != NULL || subpsi->nsub > 1 ||
subpsi->sub[0]->vdim[0] != 1)
&& !equalsSymmetric(OWNISO(0))) RETURN_ERR(ERRORWRONGISO);
cov->finiterange = falsch;
setbackward(cov, cov->key);
assert(cov->finiterange != wahr); //
VDIM0 = VDIM1 = 1;
COV_DELETE(cov->sub + 0);
if ((err = covcpy(cov->sub + 0, cov->key->sub[0])) != NOERROR) BUG;
SET_CALLING(cov->sub[0], cov);
// kein setbackward(cov, subpsi);
EXTRA_STORAGE;
RETURN_NOERROR;
}
bool allowedDgennsst(model *cov) {
allowedD(cov->sub[1]);
MEMCOPY(cov->allowedD, cov->sub[1]->allowedD, sizeof(allowedD_type));
return false;
}
bool allowedIgennsst(model *cov) {
bool *I = cov->allowedI;
model *subpsi = cov->sub[1];
for (int i=(int) FIRST_ISOUSER; i<=(int) LAST_ISOUSER; I[i++] = false);
I[(int) SYMMETRIC] = true;
I[(int) DOUBLEISOTROPIC] = MODELNR(subpsi) == MATRIX && PisNULL(M_M) &&
subpsi->kappasub[M_M] == NULL && subpsi->nsub > 1;
return false;
}
void rangegennsst(model VARIABLE_IS_NOT_USED *cov, range_type* ra){
ra->min[GENNSST_DIM_U] = 1;
ra->max[GENNSST_DIM_U] = OWNLOGDIM(0) - 1;
ra->pmin[GENNSST_DIM_U] = 1;
ra->pmax[GENNSST_DIM_U] = ra->max[GENNSST_DIM_U];
ra->openmin[GENNSST_DIM_U] = false;
ra->openmax[GENNSST_DIM_U] = false;
}
