https://github.com/cran/RandomFields
Tip revision: 6332d8d86088cebf6f828f1d29c71b8060f7e88b authored by Martin Schlather on 23 June 2016, 09:04:50 UTC
version 3.1.16
version 3.1.16
Tip revision: 6332d8d
plusmalS.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
Definition of auxiliary correlation functions
Note:
* Never use the below functions directly, but only by the functions indicated
in RFsimu.h, since there is gno error check (e.g. initialization of RANDOM)
* VARIANCE, SCALE are not used here
* definitions for the random coin method can be found in MPPFcts.cc
* definitions for genuinely anisotropic or nonsta tionary models are in
SophisticatedModel.cc; hyper models also in Hypermodel.cc
Copyright (C) 2001 -- 2003 Martin Schlather
Copyright (C) 2004 -- 2004 Yindeng Jiang & Martin Schlather
Copyright (C) 2005 -- 2015 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 <math.h>
#include <stdio.h>
#include <R_ext/Lapack.h>
#include "RF.h"
#include "Operator.h"
#include "variogramAndCo.h"
//#include <R_ext/Applic.h>
//#include <R_ext/Utils.h>
//#include <R_ext/BLAS.h>
// $
bool hasVarOnly(cov_model *cov) {
if (cov == NULL || !isDollar(cov)) BUG;
if (!PisNULL(DSCALE) && P0(DSCALE) != 1.0) return false;
if (!PisNULL(DANISO) || !PisNULL(DPROJ)) return false;
int i,
kappas = CovList[cov->nr].kappas;
for (i=0; i<kappas; i++)
if (cov->kappasub[i] != NULL) return false;
return true;
}
void kappaS(int i, cov_model *cov, int *nr, int *nc){
switch(i) {
case DVAR : case DSCALE :
*nr = *nc = 1;
break;
case DANISO :
*nr = cov->xdimown;
*nc = SIZE_NOT_DETERMINED;
break;
case DAUSER :
*nr = SIZE_NOT_DETERMINED;
*nc = cov->xdimown;
break;
case DPROJ :
*nr = SIZE_NOT_DETERMINED;
*nc = 1;
break;
default : *nr = *nc = -1;
}
}
// simple transformations as variance, scale, anisotropy matrix, etc.
void Siso(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[DOLLAR_SUB];
int i,
vdim = cov->vdim[0],
vdimSq = vdim * vdim;
double y,
*aniso=P(DANISO),
*scale =P(DSCALE),
var = P0(DVAR);
assert(cov->Sdollar->simplevar);
y = *x;
if (aniso != NULL) y = fabs(y * aniso[0]);
// printf("scale = %d %d\n", scale == NULL, aniso == NULL);
// printf("scale = %f %f \n", scale[0], *x);
if (scale != NULL)
y = scale[0]>0.0 ? y / scale[0]
: (y == 0.0 && scale[0]==0.0) ? 0.0 : RF_INF;
// letzteres darf nur passieren wenn dim = 1!!
COV(&y, next, v);
for (i=0; i<vdimSq; i++) v[i] *= var;
}
// simple transformations as variance, scale, anisotropy matrix, etc.
void logSiso(double *x, cov_model *cov, double *v, double *Sign){
cov_model *next = cov->sub[DOLLAR_SUB];
int i,
vdim = cov->vdim[0],
vdimSq = vdim * vdim;
double y,
*aniso=P(DANISO),
*scale =P(DSCALE),
logvar = log(P0(DVAR));
assert(cov->Sdollar->simplevar);
y = *x;
if (aniso != NULL) y = fabs(y * aniso[0]);
if (scale != NULL)
y = scale[0]>0.0 ? y / scale[0]
: (y == 0.0 && scale[0]==0.0) ? 0.0 : RF_INF;
LOGCOV(&y, next, v, Sign);
for (i=0; i<vdimSq; i++) v[i] += logvar;
}
void Sstat(double *x, cov_model *cov, double *v){
logSstat(x, cov, v, NULL);
}
void logSstat(double *x, cov_model *cov, double *v, double *Sign){
assert(cov->kappasub[DSCALE] == NULL &&
(cov->kappasub[DAUSER]==NULL ||
CovList[cov->kappasub[DAUSER]->nr].check==checkAngle));
cov_model *next = cov->sub[DOLLAR_SUB];
// *Aniso = cov->kappasub[DAUSER],
// *Scale = cov->kappasub[DSCALE];
double *z1 = NULL,
*scale =P(DSCALE),
*aniso=P(DANISO);
int i,
nproj = cov->nrow[DPROJ],
vdim = cov->vdim[0],
vdimSq = vdim * vdim;
if (nproj > 0) {
int *proj = PINT(DPROJ);
ALLOC_DOLLAR(z, nproj);
if (scale == NULL || scale[0] > 0.0) {
if (scale == NULL) for (i=0; i<nproj; i++) z[i] = x[proj[i] - 1];
else {
double invscale = 1.0 / scale[0];
for (i=0; i<nproj; i++) {
z[i] = invscale * x[proj[i] - 1];
}
}
} else {
// projection and aniso may not be given at the same time
for (i=0; i<nproj; i++)
z[i] = (x[proj[i] - 1] == 0 && scale[0] == 0) ? 0.0 : RF_INF;
}
z1 = z;
// } else if (Aniso != NULL) {
// int dim = Aniso->vdim[0];
// ALLOC_DOLLAR(z, dim);
// FCTN(x, Aniso, z);
// z1 = z;
// } else if (Scale != NULL) {
// int dim = Aniso->vdim[0];
// ALLOC_DOLLAR(z, dim);
// FCTN(x, Aniso, z);
// z1 = z;
} else if (aniso==NULL && (scale == NULL || scale[0] == 1.0)) {
z1 = x;
} else {
int xdimown = cov->xdimown;
double *xz;
// printf("xdimown %d\n", xdimown); BUG;
ALLOC_DOLLAR(z, xdimown);
if (aniso!=NULL) {
xA(x, aniso, cov->nrow[DANISO], cov->ncol[DANISO], z);
xz = z;
} else xz = x;
if (scale != NULL) {
if (scale[0] > 0.0) {
double invscale = 1.0 / scale[0];
for (i=0; i < xdimown; i++) z[i] = invscale * xz[i];
} else {
for (i=0; i < xdimown; i++)
z[i] = (xz[i] == 0.0 && scale[0] == 0.0) ? 0.0 : RF_INF;
}
}
z1 = z;
}
double var;
if (cov->Sdollar->simplevar) {
var = P0(DVAR);
} else {
FCTN(x, cov->kappasub[DVAR], &var);
}
if (Sign==NULL) {
COV(z1, next, v);
for (i=0; i<vdimSq; i++) v[i] *= var;
} else {
LOGCOV(z1, next, v, Sign);
var = log(var);
for (i=0; i<vdimSq; i++) v[i] += var;
}
}
void Snonstat(double *x, double *y, cov_model *cov, double *v){
logSnonstat(x, y, cov, v, NULL);
}
void logSnonstat(double *x, double *y, cov_model *cov, double *v, double *Sign){
cov_model
*next = cov->sub[DOLLAR_SUB],
*Aniso = cov->kappasub[DAUSER],
*Scale = cov->kappasub[DSCALE];
double *z1, *z2,
s1 = RF_NA, s2 = RF_NA, smeanSq=RF_NA,
*scale =P(DSCALE),
*aniso=P(DANISO);
int i,
nproj = cov->nrow[DPROJ],
vdim = cov->vdim[0],
vdimSq = vdim * vdim;
if (nproj > 0) {
int *proj = PINT(DPROJ);
ALLOC_DOLLARY(Z1, Z2, nproj);
z1 = Z1;
z2 = Z2;
if (scale==NULL || scale[0] > 0.0) {
double invscale = scale==NULL ? 1.0 : 1.0 / scale[0];
for (i=0; i<nproj; i++) {
z1[i] = invscale * x[proj[i] - 1];
z2[i] = invscale * y[proj[i] - 1];
}
} else {
double s = scale[0]; // kann auch negativ sein ...
for (i=0; i<nproj; i++) {
z1[i] = (x[proj[i] - 1] == 0.0 && s == 0.0) ? 0.0 : RF_INF;
z2[i] = (y[proj[i] - 1] == 0.0 && s == 0.0) ? 0.0 : RF_INF;
}
}
} else if (Aniso != NULL) {
int dim = Aniso->vdim[0];
ALLOC_DOLLARY(Z1, Z2, dim);
z1 = Z1;
z2 = Z2;
FCTN(x, Aniso, z1);
FCTN(y, Aniso, z2);
} else if (Scale != NULL && !isRandom(Scale)) {
int xdimown = cov->xdimown;
double s;
ALLOC_DOLLARY(Z1, Z2, xdimown);
z1 = Z1;
z2 = Z2;
FCTN(x, Scale, &s1);
FCTN(y, Scale, &s2);
if (s1 <= 0.0 || s2 <= 0.0)
ERR1("'%s' must be a positive function", KNAME(DSCALE));
smeanSq = 0.5 * (s1 * s1 + s2 * s2);
s = sqrt(smeanSq);
for (i=0; i<xdimown; i++) {
z1[i] = x[i] / s;
z2[i] = y[i] / s;
}
} else if (aniso==NULL && (scale==NULL || scale[0] == 1.0)) {
z1 = x;
z2 = y;
} else {
int xdimown = cov->xdimown;
double *xz1, *xz2;
ALLOC_DOLLARY(Z1, Z2, xdimown);
z1 = Z1;
z2 = Z2;
if (aniso != NULL) {
xA( x, y, aniso,cov->nrow[DANISO], cov->ncol[DANISO], z1, z2);
xz1 = z1;
xz2 = z2;
} else {
xz1 = x;
xz2 = y;
}
if (scale != NULL) {
double s = scale[0];
if (s > 0.0) {
double invscale = 1.0 / s;
for (i=0; i<xdimown; i++) {
z1[i] = invscale * xz1[i];
z2[i] = invscale * xz2[i];
}
} else {
for (i=0; i<nproj; i++) {
z1[i] = (xz1[i] == 0.0 && s == 0.0) ? 0.0 : RF_INF;
z2[i] = (xz2[i] == 0.0 && s == 0.0) ? 0.0 : RF_INF;
}
}
}
}
double var;
if (cov->Sdollar->simplevar) {
var = P0(DVAR);
if (Sign != NULL) var = log(var);
} else {
double w;
location_type *loc = Loc(cov);
int dummy = loc->i_row;
loc->i_row = loc->i_col;
FCTN(y, cov->kappasub[DVAR], &w);
loc->i_row = dummy;
FCTN(x, cov->kappasub[DVAR], &var);
var *= w;
var = Sign == NULL ? sqrt(var) : 0.5 * log(var);
}
if (Scale != NULL) {
if (Sign != NULL) var += 0.5 * log(s1 * s2 / smeanSq);
else var *= sqrt(s1 * s2 / smeanSq);
}
if (Sign == NULL) {
NONSTATCOV(z1, z2, next, v);
// printf("S-dim %d %f %f %f %f\n", vdimSq, v[0], v[1], v[2], v[3]);
for (i=0; i<vdimSq; i++) v[i] *= var;
// printf("--> S-dim %d %f %f %f %f\n", vdimSq, v[0], v[1], v[2], v[3]);
} else {
LOGNONSTATCOV(z1, z2, next, v, Sign);
for (i=0; i<vdimSq; i++) {
v[i] += var;
//printf("%f %f \n", v[i], var);
}
}
}
void covmatrixS(cov_model *cov, double *v) {
location_type *loc = Loc(cov);
cov_model *next = cov->sub[DOLLAR_SUB];
location_type *locnext = Loc(next);
assert(locnext != NULL);
int i, tot, totSq,
dim = loc->timespacedim,
vdim = cov->vdim[0];
assert(dim == cov->tsdim);
// printf("hier\n");
if ((!PisNULL(DSCALE) && P0(DSCALE) != 1.0) ||
!PisNULL(DANISO) || !PisNULL(DPROJ) ||
cov->kappasub[DSCALE] != NULL ||
cov->kappasub[DAUSER] != NULL ||
cov->kappasub[DPROJ] != NULL
) {
cov_model *prev = cov->calling;
if (prev == NULL || (!isInterface(prev) && !isProcess(prev))) prev = cov;
if (prev->Spgs == NULL && alloc_cov(prev, dim, vdim, vdim) != NOERROR)
ERR("memory allocation error in 'covmatrixS'");
// printf("hier B\n");
CovarianceMatrix(cov, v);
// printf("hierC\n");
return;
}
// printf("hier A\n");
if (cov->Spgs == NULL && alloc_cov(cov, dim, vdim, vdim) != NOERROR)
ERR("memory allocation error in 'covmatrixS'");
if (next->xdimprev != next->xdimown) {
BUG; // fuehrt zum richtigen Resultat, sollte aber nicht
// vorkommen!
CovarianceMatrix(cov, v);
return;
}
int next_gatter = next->gatternr,
next_xdimgatter = next->xdimgatter,
next_xdim = next->xdimprev;
next->gatternr = cov->gatternr;
next->xdimprev = cov->xdimprev;
next->xdimgatter = cov->xdimgatter;
CovList[next->nr].covmatrix(next, v);//hier wird uU next->totalpoints gesetzt
next->gatternr = next_gatter;
next->xdimgatter = next_xdimgatter;
next->xdimprev = next_xdim;
tot = cov->vdim[0] * locnext->totalpoints;
totSq = tot * tot;
if (cov->Sdollar->simplevar) {
double var = P0(DVAR);
if (var == 1.0) return;
for (i=0; i<totSq; v[i++] *= var);
} else {
BUG;
}
}
char iscovmatrixS(cov_model *cov) {
cov_model *sub = cov->sub[DOLLAR_SUB];
return (int) ((PisNULL(DSCALE) || P0(DSCALE) ==1.0) &&
PARAMisNULL(cov, DAUSER) &&
PARAMisNULL(cov, DPROJ) &&
cov->Sdollar->simplevar && // to do
PARAMisNULL(cov, DANISO)) * CovList[sub->nr].is_covariance(sub);
}
void DS(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[DOLLAR_SUB];
assert( cov->kappasub[DAUSER] == NULL && cov->kappasub[DSCALE] == NULL);
int i,
vdim = cov->vdim[0],
vdimSq = vdim * vdim,
nproj = cov->nrow[DPROJ];
double y[2], varSc,
*scale =P(DSCALE),
*aniso=P(DANISO),
spinvscale = 1.0;
assert(cov->Sdollar->simplevar);
assert(isCartesian(cov->isoown));
if (aniso != NULL) {
spinvscale *= aniso[0];
// was passiert, wenn es aniso nicht vom TypeIso ist ??
}
if (scale != NULL) spinvscale /= scale[0];
varSc = P0(DVAR) * spinvscale;
if (nproj == 0) {
y[0] = x[0] * spinvscale;
y[1] = (cov->isoown==ISOTROPIC || cov->ncol[DANISO]==1) ? 0.0
: x[1] * aniso[3]; // temporal; temporal scale
} else {
BUG;
// for (i=0; i<nproj; i++) {
// y[i] = spinvscale * x[proj[i] - 1];
}
Abl1(y, next, v);
for (i=0; i<vdimSq; i++) v[i] *= varSc;
}
void DDS(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[DOLLAR_SUB];
assert(cov->kappasub[DAUSER] == NULL);
int i,
vdim = cov->vdim[0],
vdimSq = vdim * vdim,
nproj = cov->nrow[DPROJ],
*proj = PINT(DPROJ);
double y[2], varScSq,
*scale =P(DSCALE),
*aniso=P(DANISO),
spinvscale = 1.0;
assert(isCartesian(cov->isoown));
assert(cov->Sdollar->simplevar);
if (aniso != NULL) {
spinvscale *= aniso[0];
// was passiert, wenn es aniso nicht vom TypeIso ist ??
}
if (scale != NULL) spinvscale /= scale[0];
varScSq = P0(DVAR) * spinvscale * spinvscale;
if (nproj == 0) {
y[0] = x[0] * spinvscale;
y[1] = (cov->isoown==ISOTROPIC || cov->ncol[DANISO]==1) ? 0.0
: x[1] * aniso[3]; // temporal; temporal scale
} else {
BUG;
for (i=0; i<nproj; i++) {
y[0] += x[proj[i] - 1] * x[proj[i] - 1];
}
y[0] = sqrt(y[0]) * spinvscale;
}
Abl2(y, next, v);
for (i=0; i<vdimSq; i++) v[i] *= varScSq;
}
void D3S(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[DOLLAR_SUB];
assert(cov->kappasub[DAUSER] == NULL);
int i,
vdim = cov->vdim[0],
vdimSq = vdim * vdim,
nproj = cov->nrow[DPROJ],
*proj = PINT(DPROJ);
double y[2], varScS3,
*scale =P(DSCALE),
*aniso=P(DANISO),
spinvscale = 1.0;
assert(isCartesian(cov->isoown));
assert(cov->Sdollar->simplevar);
if (aniso != NULL) {
spinvscale *= aniso[0];
// was passiert, wenn es aniso nicht vom TypeIso ist ??
}
if (scale != NULL) spinvscale /= scale[0];
varScS3 = P0(DVAR) * spinvscale * spinvscale * spinvscale;
if (nproj == 0) {
y[0] = x[0] * spinvscale;
y[1] = (cov->isoown==ISOTROPIC || cov->ncol[DANISO]==1) ? 0.0
: x[1] * aniso[3]; // temporal; temporal scale
} else {
BUG;
for (i=0; i<nproj; i++) {
y[0] += x[proj[i] - 1] * x[proj[i] - 1];
}
y[0] = sqrt(y[0]) * spinvscale;
}
Abl3(y, next, v);
for (i=0; i<vdimSq; i++) v[i] *= varScS3;
}
void D4S(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[DOLLAR_SUB];
assert(cov->kappasub[DAUSER] == NULL);
int i,
vdim = cov->vdim[0],
vdimSq = vdim * vdim,
nproj = cov->nrow[DPROJ],
*proj = PINT(DPROJ);
double y[2], varScS4,
*scale =P(DSCALE),
*aniso=P(DANISO),
spinvscale = 1.0;
assert(isCartesian(cov->isoown));
if (aniso != NULL) {
spinvscale *= aniso[0];
// was passiert, wenn es aniso nicht vom TypeIso ist ??
}
if (scale != NULL) spinvscale /= scale[0];
varScS4 = spinvscale * spinvscale;
varScS4 *= varScS4 * P0(DVAR);
if (nproj == 0) {
y[0] = x[0] * spinvscale;
y[1] = (cov->isoown==ISOTROPIC || cov->ncol[DANISO]==1) ? 0.0
: x[1] * aniso[3]; // temporal; temporal scale
} else {
BUG;
for (i=0; i<nproj; i++) {
y[0] += x[proj[i] - 1] * x[proj[i] - 1];
}
y[0] = sqrt(y[0]) * spinvscale;
}
Abl3(y, next, v);
for (i=0; i<vdimSq; i++) v[i] *= varScS4;
}
void nablahessS(double *x, cov_model *cov, double *v, bool nabla){
cov_model *next = cov->sub[DOLLAR_SUB],
*Aniso = cov->kappasub[DAUSER];
if (Aniso != NULL) BUG;
int i, endfor,
dim = cov->nrow[DANISO],// == ncol == x d i m ?
xdimown = cov->xdimown,
nproj = cov->nrow[DPROJ];
double *xy, *vw,
*scale =P(DSCALE),
*aniso=P(DANISO),
var = P0(DVAR);
if (nproj != 0) BUG;
if (dim != xdimown) BUG;
if (!cov->Sdollar->simplevar)
NotProgrammedYet("nabla not programmed for arbitrary 'var'");
if (aniso != NULL) {
ALLOC_DOLLAR(z, xdimown);
ALLOC_DOLLAR2(w, xdimown);
xA(x, aniso, xdimown, xdimown, z);
xy = z;
vw = w;
} else {
xy = x;
vw = v;
}
if (scale != NULL) {
ALLOC_DOLLAR3(y, xdimown);
assert(scale[0] > 0.0);
double spinvscale = 1.0 / scale[0];
var *= spinvscale;
if (!nabla) var *= spinvscale; // gives spinvscale^2 in case of hess
for (i=0; i<xdimown; i++) y[i] = xy[i] * spinvscale;
xy = y;
}
endfor = xdimown;
if (nabla) {
NABLA(xy, next, vw);
} else {
HESSE(xy, next, vw);
endfor *= xdimown;
}
if (aniso != NULL) {
if (nabla) Ax(aniso, vw, xdimown, xdimown, v);
else {
XCXt(aniso, vw, v, xdimown, xdimown); //todo:?reprogramm XCXt with alloc here ?
}
}
for (i=0; i<endfor; i++) v[i] *= var;
}
void nablaS(double *x, cov_model *cov, double *v){
nablahessS(x, cov, v, true);
}
void hessS(double *x, cov_model *cov, double *v){
nablahessS(x, cov, v, false);
}
void inverseS(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[DOLLAR_SUB];
int i,
idx[2] = {DAUSER, DPROJ};
double
scale;
if (cov->kappasub[DVAR] != NULL)
NotProgrammedYet("nabla not programmed for arbitrary 'var'");
for (i=0; i<2; i++) {
if (cov->kappasub[idx[i]] != NULL) {
char Msg[100];
sprintf(Msg, "inverse can only be calculated if '%s' is not an arbitrary function",
KNAME(idx[i]));
ERR(Msg);
}
}
if (cov->kappasub[DSCALE] != NULL) {
double left;
NONSTATINVERSE(ZERO, cov->kappasub[DSCALE], &left, &scale);
if (left < 0.0) ERR("scale not defined to be non-negative.");
} else scale = PisNULL(DSCALE) ? 1.0 : P0(DSCALE);
int
dim = cov->xdimown,
nproj = cov->nrow[DPROJ];
// *proj = (int *)P(DPROJ];
double y,
s = 1.0,
*aniso=P(DANISO),
var = P0(DVAR);
//PMI(cov); crash();
if (dim != 1) BUG;
if (aniso != NULL) {
if (isMiso(Type(aniso, cov->nrow[DANISO], cov->ncol[DANISO])))
s /= aniso[0];
else NotProgrammedYet(""); // to do
}
s *= scale;
if (nproj == 0) {
y= *x / var; // inversion, so variance becomes scale
} else {
BUG; //ERR("nproj is not allowed in invS");
}
if (CovList[next->nr].inverse == ErrCov) BUG;
INVERSE(&y, next, v);
for (i=0; i<dim; i++) v[i] *= s; //!
}
void nonstatinverseS(double *x, cov_model *cov, double *left, double*right,
bool log){
cov_model
*next = cov->sub[DOLLAR_SUB],
*Aniso = cov->kappasub[DAUSER],
*Scale = cov->kappasub[DSCALE];
int i,
dim = cov->tsdim,
nproj = cov->nrow[DPROJ];
// *proj = (int *)P(DPROJ];
double y,
s = 1.0,
*scale =P(DSCALE),
*aniso=P(DANISO),
var = P0(DVAR);
if (cov->kappasub[DVAR] != NULL)
NotProgrammedYet("nabla not programmed for arbitrary 'var'");
if (nproj == 0) {
y= *x / var; // inversion, so variance becomes scale
} else {
BUG; //ERR("nproj is not allowed in invS");
}
if (CovList[next->nr].nonstat_inverse == ErrInverseNonstat) BUG;
if (log) {
NONSTATLOGINVERSE(&y, next, left, right);
} else {
NONSTATINVERSE(&y, next, left, right);
}
if (aniso != NULL) {
if (isMiso(Type(aniso, cov->nrow[DANISO], cov->ncol[DANISO]))) s/=aniso[0];
else {
dollar_storage *S = cov->Sdollar;
int
ncol = cov->ncol[DANISO],
nrow = cov->nrow[DANISO],
ncnr = ncol * nrow,
bytes = ncnr * sizeof(double),
size = ncol * sizeof(double);
bool redo;
if (ncol != nrow) BUG;
if ((redo = S->save_aniso == NULL)) {
S->save_aniso = (double *) MALLOC(bytes);
S->inv_aniso = (double *) MALLOC(bytes);
}
ALLOC_DOLLAR4(LR, ncol);
double *save = S->save_aniso,
*inv = S->inv_aniso;
if (!redo) {
for (i=0; i<ncnr; i++) if ((redo = save[i] != P(DANISO)[i])) break;
}
if (redo) {
MEMCOPY(save, P(DANISO), bytes);
MEMCOPY(inv, P(DANISO), bytes);
if (Ext_invertMatrix(inv, nrow) != NOERROR)
ERR("inversion of anisotropy matrix failed");
}
MEMCOPY(LR, right, size);
xA(LR, inv, nrow, ncol, right);
MEMCOPY(LR, left, size);
xA(LR, inv, nrow, ncol, left);
}
}
if (Aniso != NULL) {
if (aniso != NULL) BUG;
if (CovList[Aniso->nr].inverse == ErrInverse)
ERR("inverse of anisotropy matrix function unknown");
int
nrow = Aniso->vdim[0],
ncol = Aniso->vdim[1],
size = nrow * sizeof(double);
if (cov->xdimown != ncol || ncol != 1)
ERR("anisotropy function not of appropriate form");
ALLOC_DOLLAR4(LR, nrow);
MEMCOPY(LR, right, size);
INVERSE(LR, Aniso, right);
MEMCOPY(LR, left, size);
INVERSE(LR, Aniso, left);
}
if (Scale != NULL && !isRandom(Scale)) {
double dummy;
COV(ZERO, Scale, &dummy);
s *= dummy;
} else if (scale != NULL) s *= scale[0];
if (s != 1.0) {
for (i=0; i<dim; i++) {
left[i] *= s; //!
right[i] *= s;
}
}
}
void nonstatinverseS(double *x, cov_model *cov, double *left, double*right) {
nonstatinverseS(x, cov, left, right, false);
}
void nonstat_loginverseS(double *x, cov_model *cov, double *left, double*right){
nonstatinverseS(x, cov, left, right, true);
}
void coinitS(cov_model *cov, localinfotype *li) {
assert(cov->Sdollar->simplevar);
cov_model *next = cov->sub[DOLLAR_SUB];
if ( CovList[next->nr].coinit == NULL)
ERR("# cannot find coinit -- please inform author");
CovList[next->nr].coinit(next, li);
}
void ieinitS(cov_model *cov, localinfotype *li) {
assert(cov->Sdollar->simplevar);
cov_model *next = cov->sub[DOLLAR_SUB];
if ( CovList[next->nr].ieinit == NULL)
ERR("# cannot find ieinit -- please inform author");
CovList[next->nr].ieinit(next, li);
}
void tbm2S(double *x, cov_model *cov, double *v){
assert(cov->Sdollar->simplevar);
cov_model *next = cov->sub[DOLLAR_SUB];
cov_fct *C = CovList + next->nr; // kein gatternr, da isotrop
double y[2], *xy,
*scale =P(DSCALE),
*aniso = P(DANISO);
assert(cov->kappasub[DAUSER] == NULL);
assert(cov->nrow[DPROJ] == 0);
if (aniso!=NULL) {
if (cov->ncol[DANISO]==2) { // turning layers
y[0] = x[0] * aniso[0]; // spatial
y[1] = x[1] * aniso[3]; // temporal
assert(aniso[1] == 0.0 && aniso[2] == 0.0);
if (y[0] == 0.0) C->cov(y, next, v);
else C->tbm2(y, next, v);
} else {
assert(cov->ncol[DANISO]==1);
if (cov->nrow[DANISO] == 1) { // turning bands
y[0] = x[0] * aniso[0]; // purely spatial
C->tbm2(y, next, v);
} else { // turning layers, dimension reduction
if (P0(DANISO) == 0.0) {
y[0] = x[1] * aniso[1]; // temporal
C->cov(y, next, v);
} else {
y[0] = x[0] * aniso[0]; // temporal
C->tbm2(y, next, v);
}
}
}
xy = y;
} else xy = x;
if (scale != NULL) {
double s = scale[0];
if (s > 0) {
double invscale = 1.0 / s;
if (cov->xdimown == 2){
y[0] = xy[0] * invscale; // spatial
y[1] = xy[1] * invscale; // temporal
if (y[0] == 0.0) C->cov(y, next, v);
else C->tbm2(y, next, v);
} else {
y[0] = xy[0] * invscale; // purely spatial
C->tbm2(y, next, v);
}
} else {
y[0] = (s < 0.0 || xy[0] != 0.0) ? RF_INF : 0.0;
if (cov->xdimown == 2)
y[1] = (s < 0.0 || xy[1] != 0.0) ? RF_INF : 0.0;
C->tbm2(y, next, v);
}
}
*v *= P0(DVAR);
}
// TODO : Aniso=Matrix: direkte implementierung in S,
// sodass nicht ueber initS gegangen werden muss, sondern
// e < - e^\top Aniso
int TaylorS(cov_model *cov) {
cov_model
*next = cov->sub[DOLLAR_SUB],
*sub = cov->key == NULL ? next : cov->key;
int i;
if (PisNULL(DPROJ) && PisNULL(DANISO)) {
double scale = PisNULL(DSCALE) ? 1.0 : P0(DSCALE);
cov->taylorN = sub->taylorN;
for (i=0; i<cov->taylorN; i++) {
cov->taylor[i][TaylorPow] = sub->taylor[i][TaylorPow];
cov->taylor[i][TaylorConst] = sub->taylor[i][TaylorConst] *
P0(DVAR) * pow(scale, -sub->taylor[i][TaylorPow]);
}
cov->tailN = sub->tailN;
for (i=0; i<cov->tailN; i++) {
cov->tail[i][TaylorPow] = sub->tail[i][TaylorPow];
cov->tail[i][TaylorExpPow] = sub->tail[i][TaylorExpPow];
cov->tail[i][TaylorConst] = sub->tail[i][TaylorConst] *
P0(DVAR) * pow(scale, -sub->tail[i][TaylorPow]);
cov->tail[i][TaylorExpConst] = sub->tail[i][TaylorExpConst] *
pow(scale, -sub->tail[i][TaylorExpPow]);
}
} else {
cov->taylorN = cov->tailN = 0;
}
return NOERROR;
}
int checkS(cov_model *cov) {
// hier kommt unerwartet ein scale == nan rein ?!!
cov_model
*next = cov->sub[DOLLAR_SUB],
*Aniso = cov->kappasub[DAUSER],
*Scale = cov->kappasub[DSCALE],
*sub = cov->key == NULL ? next : cov->key;
isotropy_type isoown = cov->isoown;
int i, err,
xdimown = cov->xdimown,
xdimNeu = xdimown,
// KOMMENTAR NICHT LOESCHEN
// *proj = PINT(DPROJ), // auf keinen Fall setzen, da Pointer unten neu
// gesetzt wird!!!!
nproj = cov->nrow[DPROJ];
// bool skipchecks = GLOBAL.general.skipchecks;
matrix_type type = TypeMany;
assert(isAnyDollar(cov));
if (!isDollarProc(cov)) cov->nr = DOLLAR; // wegen nr++ unten !
// cov->q[1] not NULL then A has been given
// if (cov->method == SpectralTBM && cov->xdimown != next->xdimprev)
// return ERRORDIM;
if ((err = checkkappas(cov, false)) != NOERROR) {
return err;
}
kdefault(cov, DVAR, 1.0);
bool angle = Aniso!=NULL && CovList[Aniso->nr].check == checkAngle;
if (angle && PisNULL(DANISO) && PisNULL(DAUSER) &&
cov->xdimown == cov->tsdim) {
int dim = cov->tsdim;
ASSERT_CARTESIAN;
if (!isCartesian(cov->isoown)) return ERRORANISO;
if ((err = CHECK(Aniso, dim, dim, ShapeType, XONLY,
CARTESIAN_COORD, SUBMODEL_DEP, cov->role)) == NOERROR) {
if (verysimple(Aniso)) {
PALLOC(DAUSER, dim, dim);
AngleMatrix(Aniso, P(DAUSER));
COV_DELETE(cov->kappasub + DAUSER);
Aniso = cov->kappasub[DAUSER] = NULL;
angle = false;
}
}
}
if (cov->kappasub[DANISO] != NULL)
SERR1("'%s' may not be a function.", KNAME(DANISO));
if (!PisNULL(DAUSER)) {
if (!isCartesian(cov->isoown)) return ERRORANISO;
// if (GLOBAL.internal.warn_Aniso) {
// PRINTF("NOTE! Starting with RandomFields 3.0, the use of '%s' is different from\nthe former '%s' insofar that '%s' is multiplied from the right by 'x' (i.e. Ax),\nwhereas '%s' had been multiplied from the left by 'x' (i.e. xA).\n", KNAME(DAUSER), KNAME(DANISO), KNAME(DANISO), KNAME(DAUSER));
// }
GLOBAL.internal.warn_Aniso = false;
// here for the first time
if (!PisNULL(DANISO)) return ERRORANISO_T;
int k,
lnrow = cov->nrow[DAUSER],
lncol = cov->ncol[DAUSER];
long j,
total = lncol * lnrow;
double
*pA = P(DAUSER);
PALLOC(DANISO, lncol, lnrow); // !! ACHTUNG col, row gekreuzt
for (i=k=0; i<lnrow; i++) {
for (j=i; j<total; j+=lnrow) P(DANISO)[k++] = pA[j];
}
PFREE(DAUSER);
}
bool simplevar = cov->kappasub[DVAR] == NULL || isRandom(cov->kappasub[DVAR]);
if (!simplevar) {
//cov->domown = KERNEL;
ptwise_type ptt = cov->ptwise_definite;
isotropy_type isonew = UpgradeToCoordinateSystem(cov->isoown);
if (isonew == ISO_MISMATCH) SERR("reduced systems not allowed");
if ((err = CHECK(cov->kappasub[DVAR], cov->tsdim, cov->xdimown,
ShapeType, // only!! -- for pos def use RMprod
XONLY, isonew,
SCALAR, ROLE_BASE)) != NOERROR) return err;
if (cov->kappasub[DVAR]->ptwise_definite != pt_posdef) {
if (cov->kappasub[DVAR]->ptwise_definite == pt_unknown) {
if (GLOBAL.internal.warn_negvar && cov->q==NULL) {
QALLOC(1);
warning("positivity of the variance in '%s' cannot be detected.",
NICK(next));
}
} else {
SERR2("positivity of '%s' required. Got '%s'", KNAME(DVAR),
POSITIVITY_NAMES[cov->kappasub[DVAR]->ptwise_definite]);
}
}
cov->ptwise_definite = ptt;
}
DOLLAR_STORAGE;
if (Aniso != NULL) {
if (!isDollarProc(cov) && !angle && cov->domown != KERNEL)
return ERRORFAILED;
if (!PisNULL(DANISO) || !PisNULL(DPROJ) || !PisNULL(DSCALE) ||
(Scale!=NULL && !isRandom(Scale)) )
SERR2("if '%s' is an arbitrary function, only '%s' may be given as additional parameter.", KNAME(DAUSER), KNAME(DVAR));
if (cov->isoown != SYMMETRIC && !isCoordinateSystem(cov->isoown)) {
return ERRORANISO;
}
cov->full_derivs = cov->rese_derivs = 0;
cov->loggiven = true;
isotropy_type isonew = UpgradeToCoordinateSystem(cov->isoown);
if (isonew == ISO_MISMATCH) SERR("reduced systems not allowed");
if ((err = CHECK(Aniso, cov->tsdim, cov->xdimown, ShapeType, XONLY,
isonew, SUBMODEL_DEP, cov->role)) != NOERROR) {
return err;
}
if (Aniso->vdim[1] != 1)
SERR4("'%s' returns a matrix of dimension %d x %d, but dimension %d x 1 is need.",
KNAME(DANISO), Aniso->vdim[0], Aniso->vdim[1], cov->xdimown);
if (cov->key==NULL) {
if ((err = CHECK(sub, Aniso->vdim[0], Aniso->vdim[0], cov->typus,
cov->domown,
cov->isoown, SUBMODEL_DEP, cov->role)) != NOERROR) {
return err;
}
}
cov->pref[Nugget] = cov->pref[RandomCoin] = cov->pref[Average] =
cov->pref[Hyperplane] = cov->pref[SpectralTBM] = cov->pref[TBM] =
PREF_NONE;
sub->pref[CircEmbed] = sub->pref[CircEmbedCutoff] =
sub->pref[CircEmbedIntrinsic] = sub->pref[Sequential] =
sub->pref[Specific] = PREF_NONE;
} else if (Scale != NULL && !isRandom(Scale)) {
if (!isDollarProc(cov) && cov->domown != KERNEL) return ERRORFAILED;
if (!PisNULL(DANISO) || !PisNULL(DPROJ) || !PisNULL(DSCALE))
SERR2("if '%s' is an arbitrary function, only '%s' may be given as additional parameter.", KNAME(DSCALE), KNAME(DVAR));
if (cov->isoown != SYMMETRIC && !isCoordinateSystem(cov->isoown)) {
return ERRORANISO;
}
cov->full_derivs = cov->rese_derivs = 0;
cov->loggiven = true;
isotropy_type isonew = UpgradeToCoordinateSystem(cov->isoown);
if (isonew == ISO_MISMATCH) SERR("reduced systems not allowed");
if ((err = CHECK(Scale, cov->tsdim, cov->xdimown, ShapeType, XONLY,
isonew, SUBMODEL_DEP, cov->role)) != NOERROR) {
return err;
}
if (Scale->vdim[1] != 1 || Scale->vdim[0] != 1)
SERR3("'%s' must be scalar, not %d x %d",
KNAME(DSCALE), Aniso->vdim[0], Aniso->vdim[1]);
if (cov->key==NULL) {
if ((err = CHECK(sub, Scale->vdim[0], Scale->vdim[0], cov->typus,
cov->domown,
cov->isoown, SUBMODEL_DEP, cov->role)) != NOERROR) {
return err;
}
if (!isNormalMixture(next))
SERR("scale function only allowed for normal mixtures.");
}
cov->pref[Nugget] = cov->pref[RandomCoin] = cov->pref[Average] =
cov->pref[Hyperplane] = cov->pref[SpectralTBM] = cov->pref[TBM] =
PREF_NONE;
sub->pref[CircEmbed] = sub->pref[CircEmbedCutoff] =
sub->pref[CircEmbedIntrinsic] = sub->pref[Sequential] =
sub->pref[Specific] = PREF_NONE;
} else if (!PisNULL(DANISO)) { // aniso given
int
nrow = cov->nrow[DANISO],
ncol = cov->ncol[DANISO];
if (nrow==0 || ncol==0) SERR("dimension of the matrix is 0");
if (!PisNULL(DPROJ)) return ERRORANISO_T;
if (xdimown < nrow) {
if (PL >= PL_ERRORS) {LPRINT("xdim=%d != nrow=%d\n", xdimown, nrow);}
SERR("#rows of anisotropy matrix does not match dim. of coordinates");
}
if (xdimown != cov->tsdim && nrow != ncol)
SERR("non-quadratic anisotropy matrices only allowed if dimension of coordinates equals spatio-temporal dimension");
type = Type(P(DANISO), nrow, ncol);
cov->full_derivs = cov->rese_derivs = 0;
cov->loggiven = true;
switch (cov->isoown) {
case ISOTROPIC :
if (cov->tsdim != 1) return ERRORANISO;
cov->full_derivs = cov->rese_derivs = 2;
break;
case SPACEISOTROPIC :
cov->full_derivs = cov->rese_derivs = isMdiag(type) ? 2 : 0;
if (nrow != 2 || !isMdiag(type)) {
SERR("spaceisotropy needs a 2x2 diagonal matrix");
}
break;
case ZEROSPACEISO :
if (!isMdiag(type)) return ERRORANISO;
break;
case VECTORISOTROPIC :
if (!isMiso(type)) return ERRORANISO;
break;
case SYMMETRIC:
break;
case PREVMODELI : BUG;
case CARTESIAN_COORD : case GNOMONIC_PROJ : case ORTHOGRAPHIC_PROJ:
if (!isProcess(cov->typus)) return ERRORANISO;
break;
default :
if (isCartesian(cov->isoown)) { BUG; }
else return ERRORANISO;
}
if (!isMdiag(type))
cov->pref[Nugget] = cov->pref[RandomCoin] = cov->pref[Average] = cov->pref[Hyperplane] = PREF_NONE;
if (cov->isoown != SPACEISOTROPIC && !isMiso(type))
cov->pref[SpectralTBM] = cov->pref[TBM] = PREF_NONE;
if (cov->key==NULL) {
if ((err = CHECK(next, ncol, ncol, cov->typus, cov->domown,
ncol==1 && !isProcess(cov->typus) ? ISOTROPIC : cov->isoown,
SUBMODEL_DEP, cov->role))
!= NOERROR) {
return err;
}
if (next->domown == cov->domown && next->isoown == cov->isoown &&
xdimown > 1) next->delflag = DEL_COV - 7;
} else {
if ((err = CHECK(cov->key, ncol, ncol, cov->typus, cov->domown,
ncol==1 &&
!isProcess(cov->typus) ? ISOTROPIC : cov->isoown,
SUBMODEL_DEP, cov->role)) != NOERROR) return err;
}
} else { // PisNULL(DANISO)
int tsdim = cov->tsdim;
if (nproj > 0) {
if (cov->Sdollar->proj < 0) { // different isoown cause different interpretations; here: restore original value
PFREE(DPROJ);
kdefault(cov, DPROJ, cov->Sdollar->proj);
nproj = 1;
}
if (P0INT(DPROJ) < 0) { // here: interprete "space" and "time"
//printf("%d\n", nproj);
if (nproj != 1) {
BUG;
SERR1("unallowed use of '%s'", KNAME(DPROJ));
}
cov->Sdollar->proj = P0INT(DPROJ);
if (Loc(cov)->Time && tsdim >= 2) {
assert(P0INT(DPROJ) == PROJ_TIME || P0INT(DPROJ) == PROJ_SPACE);
bool Time = P0INT(DPROJ) == PROJ_TIME;
PFREE(DPROJ);
if (Time) {
kdefault(cov, DPROJ, tsdim);
} else {
nproj = tsdim - 1;
PALLOC(DPROJ, nproj, 1);
for (i=1; i<=nproj; i++) PINT(DPROJ)[i-1] = i;
}
} else {
SERR1("unallowed use of '%s' or model to complicated.", KNAME(DPROJ));
}
// PMI(cov, 0);
}
if (nproj < tsdim) {
for (i = 0; i<Nothing; i++) { if (cov->pref[i] > 0) cov->pref[i] = 1; }
cov->pref[Specific] = PREF_BEST;
}
bool ok = (cov->xdimprev == cov->xdimown) ||
((cov->calling == NULL || cov->calling->isoown == EARTH_COORDS ||
cov->calling->isoown == EARTH_SYMMETRIC) &&
isCartesian(cov->isoown) && cov->xdimprev == cov->xdimown - 1);
if (!ok) {
// printf("!ok %d %d %d %d\n", cov->xdimprev, cov->xdimown, cov->isoprev == EARTH_COORDS, isCartesian(cov->isoown) );
//PMI(cov->calling);
return ERRORANISO;
}
for (i=0; i<nproj; i++) {
int j,
idx = PINT(DPROJ)[i];
if (idx <= 0)
SERR1("only positive values allowed for '%s'", KNAME(DPROJ));
if (idx > xdimown)
SERR4("%d%s value of '%s' (%d) out of range",
i + 1,
i == 0 ? "st" : i == 1 ? "nd" : i == 2 ? "rd" : "th",
KNAME(DPROJ), PINT(DPROJ)[i]);
for (j=i+1; j<nproj; j++) {
if (PINT(DPROJ)[j] == idx) { BUG;
SERR1("values of '%s' must be distinct.", KNAME(DPROJ));
}
}
}
xdimNeu = nproj;
tsdim = nproj;
switch (cov->isoown) {
case ISOTROPIC : case EARTH_ISOTROPIC : case SPHERICAL_ISOTROPIC :
if (cov->tsdim != 1) return ERRORANISO;
break;
case SPACEISOTROPIC :
if (nproj > 2) SERR("maximum length of projection vector is 2");
if (nproj == 2) {
if (PINT(DPROJ)[0] >= PINT(DPROJ)[1])
SERR1("in case of '%s' projection directions must be ordered",
ISONAMES[SPACEISOTROPIC]);
}
if (P0INT(DPROJ) == 1) {
tsdim = cov->tsdim - 1 ;
isoown = ISOTROPIC;
xdimNeu = 1;
} else {
assert(P0INT(DPROJ) == 2);
tsdim = 1;
isoown = ISOTROPIC;
xdimNeu = 1;
}
break;
case ZEROSPACEISO :
isoown = SYMMETRIC;
break;
case VECTORISOTROPIC : SERR("projection of vectorisotropic fields not programmed yet"); // to do: vdim muss auch reduziert werden ... --- das wird
// grausam !
break;
case SYMMETRIC: case CARTESIAN_COORD:
break;
case GNOMONIC_PROJ : case ORTHOGRAPHIC_PROJ :
isoown = CARTESIAN_COORD;
break;
case PREVMODELI : BUG;
break;
case SPHERICAL_SYMMETRIC : case EARTH_SYMMETRIC :
if (nproj != 2 || PINT(DPROJ)[0] != 1 || PINT(DPROJ)[1] != 2)
isoown = SYMMETRIC;
break;
case SPHERICAL_COORDS : case EARTH_COORDS :
if (nproj != 2 || PINT(DPROJ)[0] != 1 || PINT(DPROJ)[1] != 2)
isoown = CARTESIAN_COORD;
break;
default :
if (isCartesian(cov->isoown)) {BUG;}
else return ERRORANISO; // todo
}
}
// verhindern, dass die coordinaten-transformation anlaeuft,
// da aus z.B. Erd-coordinaten durch Projektion (auf die Zeitachse)
// kartesische werden
if (cov->key==NULL) {
if ((err = CHECK_NO_TRAFO(next, tsdim, xdimNeu, cov->typus, cov->domown,
isoown,
cov->vdim[0], // SUBMODEL_DEP; geaendert 20.7.14
cov->role)) != NOERROR) {
return err;
}
if (next->domown == cov->domown &&
next->isoown == cov->isoown) // darf nicht nach C-HECK als allgemeine Regel ruebergezogen werden, u.a. nicht wegen stat/nicht-stat wechsel !!
// hier geht es, da domown und isoown nur durchgegeben werden und die Werte // bereits ein Schritt weiter oben richtig/minimal gesetzt werden.
next->delflag = DEL_COV - 8;
} else {
if ((err = CHECK_NO_TRAFO(cov->key, tsdim, xdimNeu, cov->typus,
cov->domown, isoown,
SUBMODEL_DEP, cov->role)) != NOERROR) return err;
}
// PMI(cov);
} // end no aniso
if (( err = checkkappas(cov, false)) != NOERROR) {
return err;
}
setbackward(cov, sub);
if ((Aniso != NULL || (Scale != NULL && !isRandom(Scale)) ||
!PisNULL(DANISO)|| !PisNULL(DPROJ)) && cov->maxdim < cov->xdimown)
cov->maxdim = cov->xdimown;
if (!isAnyIsotropic(cov->isoown) && !isDollarProc(cov)) { // multivariate kann auch xdimNeu == 1 problematisch sein
cov->nr++;
}
if (xdimNeu > 1) {
cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] = 0;
}
// 30.10.11 kommentiert:
// cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] =
// cov->pref[TBM] = cov->pref[SpectralTBM] = 0;
if ( (PisNULL(DANISO) || isMiso(type)) && PisNULL(DPROJ)) {
cov->logspeed = sub->logspeed * P0(DVAR);
}
////////////////
if (sub->pref[Nugget] > PREF_NONE) cov->pref[Specific] = 100;
if (nproj == 0) cov->matrix_indep_of_x = sub->matrix_indep_of_x;
if ((err = TaylorS(cov)) != NOERROR) return err;
cov->Sdollar->simplevar = simplevar;
cov->Sdollar->isoown = isoown; // for struct Sproc !
if (isProcess(cov->typus)) {
MEMCOPY(cov->pref, PREF_NOTHING, sizeof(pref_shorttype));
}
if (GLOBAL.coords.coord_system == earth &&
is_all(isCartesian, CovList + next->nr) &&
GLOBAL.internal.warn_scale &&
(PisNULL(DSCALE) ||
P0(DSCALE) < (strcmp(GLOBAL.coords.newunits[0], "km")== 0 ? 10 : 6.3))) {
char msg[300];
sprintf(msg, "value of scale parameter equals '%4.2f',\nwhich is less than 100, although models defined on R^3 are used in the\ncontext of earth coordinates where larger scales are expected.\n(This warning appears only ones per session.)\n", PisNULL(DSCALE) ? 1.0 : P0(DSCALE));
GLOBAL.internal.warn_scale = false;
warning(msg);
}
return NOERROR;
}
void rangeS(cov_model *cov, range_type* range){
int i;
bool negdef = isNegDef(cov->typus);
range->min[DVAR] = negdef ? 0.0 : RF_NEGINF;
range->max[DVAR] = RF_INF;
range->pmin[DVAR] = negdef ? 0.0 : -10000;
range->pmax[DVAR] = 100000;
range->openmin[DVAR] = !negdef;
range->openmax[DVAR] = true;
range->min[DSCALE] = 0.0;
range->max[DSCALE] = RF_INF;
range->pmin[DSCALE] = 0.0001;
range->pmax[DSCALE] = 10000;
range->openmin[DSCALE] = true;
range->openmax[DSCALE] = true;
for (i=DANISO; i<= DAUSER; i++) {
range->min[i] = RF_NEGINF;
range->max[i] = RF_INF;
range->pmin[i] = -10000;
range->pmax[i] = 10000;
range->openmin[i] = true;
range->openmax[i] = true;
}
range->min[DPROJ] = -2;
range->max[DPROJ] = cov->tsdim;
range->pmin[DPROJ] = 1;
range->pmax[DPROJ] = cov->tsdim;
range->openmin[DPROJ] = false;
range->openmax[DPROJ] = false;
}
bool TypeS(Types required, cov_model *cov, int depth) {
cov_model *sub = cov->key==NULL ? cov->sub[0] : cov->key;
// print(" TypeS %s %s : %d %d\n", TYPENAMES[required], NAME(sub), isShape(required) || isTrend(required) || isProcess(required), TypeConsistency(required, sub, depth-1));
return (isShape(required) || isTrend(required) || isProcess(required))
&& TypeConsistency(required, sub, depth-1);
}
void spectralS(cov_model *cov, gen_storage *s, double *e) {
cov_model *next = cov->sub[DOLLAR_SUB];
int d,
ncol = PisNULL(DANISO) ? cov->tsdim : cov->ncol[DANISO];
double sube[MAXTBMSPDIM],
*scale =P(DSCALE),
invscale = 1.0;
SPECTRAL(next, s, sube); // nicht gatternr
// Reihenfolge nachfolgend extrem wichtig, da invscale auch bei aniso
// verwendet wird
if (scale != NULL) invscale /= scale[0];
if (!PisNULL(DANISO)) {
int
nrow = cov->nrow[DANISO];
long j, k, m,
total = ncol * nrow;
double
*A = P(DANISO);
for (d=0, k=0; d<nrow; d++, k++) {
e[d] = 0.0;
for (m=0, j=k; j<total; j+=nrow) {
e[d] += sube[m++] * A[j] * invscale;
}
}
} else {
for (d=0; d<ncol; d++) e[d] = sube[d] * invscale;
}
}
void ScaleDollarToLoc(cov_model *to, cov_model *from,
int VARIABLE_IS_NOT_USED depth) {
assert(!PARAMisNULL(to, LOC_SCALE));
assert(isDollar(from));
assert(!PARAMisNULL(from, DSCALE));
PARAM(to, LOC_SCALE)[0] = PARAM0(from, DSCALE);
}
bool ScaleOnly(cov_model *cov) {
return isDollar(cov) &&
PisNULL(DPROJ) && cov->kappasub[DPROJ] == NULL &&
PisNULL(DANISO) && cov->kappasub[DAUSER] == NULL &&
(PisNULL(DVAR) || P0(DVAR)==1.0) && cov->kappasub[DVAR] == NULL;
}
int addScales(cov_model **newmodel, double anisoScale, cov_model *scale,
double Scale) {
// cov_model *calling = (*newmodel)->calling;
if (anisoScale != 1.0) {
addModel(newmodel, LOC);
kdefault(*newmodel, LOC_SCALE, anisoScale);
}
if (scale != NULL) {
if (!isRandom(scale)) SERR("unstationary scale not allowed yet");
addModel(newmodel, LOC);
addSetDistr(newmodel, scale->calling, ScaleDollarToLoc, true, MAXINT);
} else {
if (Scale != 1.0) {
addModel(newmodel, LOC);
kdefault(*newmodel, LOC_SCALE, Scale);
}
}
return NOERROR;
}
int structS(cov_model *cov, cov_model **newmodel) {
if (cov->role == ROLE_GAUSS && isProcess(cov->typus)) {
cov->nr = DOLLAR_PROC;
return structSproc(cov, newmodel); // kein S-TRUCT(...) !!
}
cov_model *local = NULL,
*dummy = NULL,
*next = cov->sub[DOLLAR_SUB],
*Aniso = cov->kappasub[DAUSER],
*Scale = cov->kappasub[DSCALE];
int err = NOERROR;
bool generalAniso = (Aniso != NULL && CovList[Aniso->nr].check != checkAngle)
|| (Scale != NULL && !isRandom(Scale));
ASSERT_NEWMODEL_NOT_NULL;
if (cov->kappasub[DVAR] != NULL) {
GERR2("Arbitrary functions for '%s' should be replaced by multiplicative models using '%s'", KNAME(DVAR), CovList[PROD].nick);
}
if (generalAniso)
GERR1("complicated models including arbitrary functions for '%s' cannot be simulated yet", KNAME(DAUSER));
switch (cov->role) {
case ROLE_POISSON_GAUSS :
case ROLE_SMITH :
if (!next->deterministic) GERR("random shapes not programmed yet");
if (!PisNULL(DANISO)) GERR("anisotropy parameter not allowed yet");
if (!PisNULL(DPROJ)) GERR("projections not allowed yet");
if ((err = STRUCT(next, newmodel)) > NOERROR) return err;
addModel(newmodel, DOLLAR);
if (!PisNULL(DVAR)) kdefault(*newmodel, DVAR, P0(DVAR));
if (!PisNULL(DSCALE)) kdefault(*newmodel, DSCALE, P0(DSCALE));
break;
case ROLE_MAXSTABLE : // eigentlich nur von RPSmith moeglich !
if (!next->deterministic) GERR("random shapes not programmed yet");
if (!PisNULL(DPROJ)) GERR("projections not allowed yet");
// P(DVAR) hat keine Auswirkungen
if (!PisNULL(DANISO)) {
if (!isMonotone(next) || !isIsotropic(next->isoown) ||
PisNULL(DANISO) || cov->ncol[DANISO] != cov->nrow[DANISO])
GERR("anisotropy parameter only allowed for simple models up to now");
}
assert(cov->calling != NULL);
double anisoScale;
if (!PisNULL(DANISO)) {
anisoScale = 1.0 / getMinimalAbsEigenValue(P(DANISO), cov->nrow[DANISO]);
if (!R_FINITE(anisoScale))
GERR("singular anisotropy matrices not allowed");
} else anisoScale = 1.0;
if (cov->calling->nr == SMITHPROC) {
if ((err = STRUCT(next, newmodel)) == NOERROR && *newmodel != NULL) {
assert( (*newmodel)->calling == cov);
Types type =
TypeConsistency(PointShapeType, *newmodel, 0) ? PointShapeType :
TypeConsistency(RandomType, *newmodel, 0) ? RandomType :
TypeConsistency(ShapeType, *newmodel, 0) ? ShapeType :
OtherType;
double scale = PisNULL(DSCALE) ? 1.0 : P0(DSCALE);
if (type == RandomType) { // random locations given;
// so, it must be of pgs type (or standard):
if ((err = CHECK_R(*newmodel, cov->tsdim)) != NOERROR) {
goto ErrorHandling;
}
dummy = *newmodel;
if ((err=addScales(&dummy, anisoScale, Scale, scale))!=NOERROR){
goto ErrorHandling;
}
*newmodel = NULL;
assert(cov->vdim[0] == cov->vdim[1]);
if ((err = addPointShape(newmodel, cov, dummy, cov->tsdim,
cov->vdim[0])) != NOERROR) {
goto ErrorHandling;
}
ASSERT_NEWMODEL_NOT_NULL;
(*newmodel)->calling = cov;
} else { // type not RandomType
if (type == PointShapeType &&
(err = addScales((*newmodel)->sub + PGS_LOC, anisoScale, Scale,
scale)) != NOERROR) goto ErrorHandling;
if ((err = CHECK(*newmodel, cov->tsdim, cov->xdimprev, type,
cov->domprev, cov->isoprev, cov->vdim,
ROLE_MAXSTABLE))
!= NOERROR) {
goto ErrorHandling;
}
if (type==PointShapeType) {
if ((err = PointShapeLocations(*newmodel, cov)) != NOERROR)
goto ErrorHandling;
} else if (type == ShapeType) {
dummy = *newmodel;
*newmodel = NULL;
// suche nach geeigneten locationen
if ((err = addScales(&local, anisoScale, Scale, scale))!=NOERROR)
goto ErrorHandling;
if ((err = addPointShape(newmodel, dummy, NULL, local,
cov->tsdim, cov->vdim[0]))
!= NOERROR)
goto ErrorHandling;
} else BUG;
} // ! randomtype
} else { // S TRUCT does not return anything
int err2;
if ((err2 = addPointShape(newmodel, cov, NULL, cov->tsdim,
cov->vdim[0]))
!= NOERROR) {
if (err == NOERROR) err = err2;
goto ErrorHandling;
}
err = NOERROR;
}
} else { // not from RPsmith
BUG;
//
if ((err = STRUCT(next, newmodel)) > NOERROR) return err;
}
break;
case ROLE_GAUSS :
if (cov->key != NULL) COV_DELETE(&(cov->key));
if (PrevLoc(cov)->distances)
GERR("distances do not allow for more sophisticated simulation methods");
if ((err = STRUCT(next, newmodel)) > NOERROR) return err;
addModel(newmodel, DOLLAR);
if (!PisNULL(DVAR)) kdefault(*newmodel, DVAR, P0(DVAR));
if (!PisNULL(DSCALE) ) kdefault(*newmodel, DSCALE, P0(DSCALE));
if (!next->deterministic) GERR("random shapes not programmed yet");
break;
default :
GERR2("%s : changes in scale/variance not programmed yet for '%s'",
NICK(cov), ROLENAMES[cov->role]);
}
ErrorHandling:
if (dummy != NULL) COV_DELETE(&dummy);
if (local != NULL) COV_DELETE(&local);
return err;
}
int initS(cov_model *cov, gen_storage *s){
// am liebsten wuerde ich hier die Koordinaten transformieren;
// zu grosser Nachteil ist dass GetDiameter nach trafo
// grid def nicht mehr ausnutzen kann -- umgehbar?!
cov_model *next = cov->sub[DOLLAR_SUB],
*varM = cov->kappasub[DVAR],
*scaleM = cov->kappasub[DSCALE],
*anisoM = cov->kappasub[DAUSER],
*projM = cov->kappasub[DPROJ];
int
vdim = cov->vdim[0],
nm = cov->mpp.moments,
nmvdim = (nm + 1) * vdim,
err = NOERROR;
bool
angle = anisoM != NULL && CovList[anisoM->nr].check == checkAngle,
maxstable = hasExactMaxStableRole(cov);// Realisationsweise
if (hasAnyShapeRole(cov)) { // Average !! ohne maxstable selbst !!
double
var[MAXMPPVDIM],
scale = PisNULL(DSCALE) ? 1.0 : P0(DSCALE);
int i,
dim = cov->tsdim;
if (!PisNULL(DPROJ) || !PisNULL(DANISO) ||
projM!= NULL || (anisoM != NULL && (!angle || !anisoM->deterministic))){
SERR("(complicated) anisotropy ond projection not allowed yet in Poisson related models");
}
// Achtung I-NIT_RANDOM ueberschreibt mpp.* !!
if (varM != NULL) {
int nm_neu = nm == 0 && !maxstable ? 1 : nm;
if ((err = INIT_RANDOM(varM, nm_neu, s, P(DVAR))) != NOERROR) return err;
int nmP1 = varM->mpp.moments + 1;
for (i=0; i<vdim; i++) {
int idx = i * nmP1;
var[i] = maxstable ? P0(DVAR) : varM->mpp.mM[idx + 1];
}
} else for (i=0; i<vdim; var[i++] = P0(DVAR));
if (scaleM != NULL) {
if (dim + nm < 1) SERR("found dimension <= 0");
int dim_neu = maxstable ? nm : (dim + nm) < 1 ? 1 : dim + nm;
if ((err = INIT_RANDOM(scaleM, dim_neu, s, P(DSCALE)))
!= NOERROR) return err;
scale = maxstable ? P0(DSCALE) : scaleM->mpp.mM[1];
}
if ((err = INIT(next, nm, s)) != NOERROR) return err;
for (i=0; i < nmvdim; i++) {
cov->mpp.mM[i] = next->mpp.mM[i];
cov->mpp.mMplus[i] = next->mpp.mMplus[i];
}
if (varM != NULL && !maxstable) {
for (i=0; i < nmvdim; i++) {
cov->mpp.mM[i] *= varM->mpp.mM[i];
cov->mpp.mMplus[i] *= varM->mpp.mMplus[i];
}
} else {
int j, k;
double pow_var;
for (k=j=0; j<vdim; j++) {
pow_var = 1.0;
for (i=0; i <= nm; i++, pow_var *= var[j], k++) {
cov->mpp.mM[k] *= pow_var;
cov->mpp.mMplus[k] *= pow_var;
}
}
}
if (scaleM != NULL && !maxstable) {
if (scaleM->mpp.moments < dim) SERR("moments can not be calculated.");
int j, k,
nmP1 = scaleM->mpp.moments + 1;
for (k=j=0; j<vdim; j++) {
double pow_scale = scaleM->mpp.mM[dim + j * nmP1];
for (i=0; i <= nm; i++, k++) {
cov->mpp.mM[k] *= pow_scale;
cov->mpp.mMplus[k] *= pow_scale;
}
}
} else {
double pow_scale = pow(scale, dim);
for (i=0; i < nmvdim; i++) {
cov->mpp.mM[i] *= pow_scale;
cov->mpp.mMplus[i] *= pow_scale;
}
}
if (!PisNULL(DANISO)) {
if (cov->nrow[DANISO] != cov->ncol[DANISO])
SERR("only square anisotropic matrices allowed");
double invdet = fabs(1.0 / getDet(P(DANISO), cov->nrow[DANISO]));
if (!R_FINITE(invdet))
SERR("determinant of the anisotropy matrix could not be calculated.");
for (i=0; i < nmvdim; i++) {
cov->mpp.mM[i] *= invdet;
cov->mpp.mMplus[i] *= invdet;
}
}
if (anisoM != NULL) {
double invdet;
if (angle) {
int
ncol = anisoM->vdim[0],
nrow = anisoM->vdim[1];
if (nrow != ncol) SERR("only square anisotropic matrices allowed");
double
*diag = PARAM(anisoM, ANGLE_DIAG);
if (diag != NULL) {
invdet = 1.0;
for (i=0; i<ncol; i++) invdet /= diag[i];
} else {
invdet = PARAM0(anisoM, ANGLE_RATIO);
}
} else {
SERR("general anisotropy matrices not allowed yet");
}
invdet = fabs(invdet);
if (!R_FINITE(invdet))
SERR("determinant of the anisotropy matrix could not be calculated.");
for (i=0; i < nmvdim; i++) {
cov->mpp.mM[i] *= invdet;
cov->mpp.mMplus[i] *= invdet;
}
}
if (R_FINITE(cov->mpp.unnormedmass)) {
if (vdim > 1) BUG;
cov->mpp.unnormedmass = next->mpp.unnormedmass * var[0];
} else {
for (i=0; i<vdim; i++)
cov->mpp.maxheights[i] = next->mpp.maxheights[i] * var[i];
}
}
else if (cov->role == ROLE_GAUSS) {
cov_model
*key = cov->key,
*sub = key == NULL ? next : key;
assert(sub != NULL);
assert(key == NULL || ({PMI(cov);false;}));//
if ((err=INIT(sub, 0, s)) != NOERROR) return err;
}
else if (cov->role == ROLE_BASE) {
cov_model
*key = cov->key,
*sub = key == NULL ? next : key;
assert(sub != NULL);
assert(key == NULL || ({PMI(cov);false;}));//
if ((err=INIT(sub, 0, s)) != NOERROR) return err;
} else {
if ((err=INIT(next, 0, s)) != NOERROR) return err; // e.g. from MLE
// PMI(cov->calling);
// printf("%s ", NAME(cov->calling));
// SERR("initiation of scale and/or variance failed");
}
if ((err = TaylorS(cov)) != NOERROR) return err;
return NOERROR;
}
void doS(cov_model *cov, gen_storage *s){
cov_model
*varM = cov->kappasub[DVAR],
*scaleM = cov->kappasub[DSCALE];
int i,
vdim = cov->vdim[0];
if (varM != NULL && !varM->deterministic && !isRandom(varM)) {
assert(!PisNULL(DVAR));
DO(varM, s);
}
if (scaleM != NULL && !scaleM->deterministic && !isRandom(scaleM)) {
assert(!PisNULL(DSCALE));
DO(scaleM, s);
}
if (hasAnyShapeRole(cov)) {
cov_model *next = cov->sub[DOLLAR_SUB];
DO(next, s);// nicht gatternr
for (i=0; i<vdim; i++)
cov->mpp.maxheights[i] = next->mpp.maxheights[i] * P0(DVAR);
return;
}
else if (cov->role == ROLE_GAUSS) {
double
*res = cov->rf,
sd = sqrt(P0(DVAR));
int
totalpoints = Gettotalpoints(cov);
assert(res != NULL);
if (cov->key == NULL) BUG;
if (sd != 1.0) for (i=0; i<totalpoints; i++) res[i] *= (double) sd;
return;
}
BUG;
}
int checkplusmal(cov_model *cov) {
cov_model *sub;
int i, j, err,
vdim[2] = {1, 1},
dim = cov->tsdim,
xdim = cov->xdimown,
role = cov->role;
bool plus = CovList[cov->nr].check == checkplus,
trend = isTrend(cov->typus);
// bool linearmodel = !plus && trend && cov->sub[0]->nr == CONST;
Types covtype = cov->typus;
domain_type covdom = trend ? XONLY : cov->domown;
int trendiso = UpgradeToCoordinateSystem(cov->isoown);
if (trendiso == ISO_MISMATCH) trendiso = cov->isoown;
assert(trendiso != ISO_MISMATCH);
int coviso = trend ? trendiso : cov->isoown;
assert(cov->Splus == NULL);
// printf("\n\n %s %s\n", TYPENAMES[cov->typus], ISONAMES[cov->isoprev]);
int variants = 1 +
(int) (!trend && (cov->calling == NULL || !isShape(cov->calling)));
// PMI(cov->calling);
cov->matrix_indep_of_x = true;
for (i=0; i<cov->nsub; i++) {
Types type = covtype;
domain_type dom = covdom;
int iso = coviso;
sub = cov->sub[i];
// PMI(sub); assert(sub->nr != CONST || sub->maxdim > -2);
// printf("> %s entering %s\n", NAME(cov), ISONAMES[coviso]);
//PMI(cov->calling);
// assert(covdom == KERNEL);
//print("stop\n"); if ( iso != EARTH_SYMMETRIC || covdom != KERNEL) return ERRORFAILED;
if (sub == NULL)
SERR("+ or *: named submodels are not given in a sequence!");
if (!plus && type==VariogramType) type=PosDefType;
err = ERRORTYPECONSISTENCY;
// printf("unten zu 1 !!");
for (j=0; j<variants; j++) { // nur trend als abweichender typus erlaubt
// printf(">> %s : %d %d type =%s %s \n", NAME(cov), j, variants, TYPENAMES[type], ISONAMES[ iso]);
// int tt; printf("%d %s\n", tt=TypeConsistency(type, sub, 0), NAME(sub));
if (TypeConsistency(type, sub, 0) &&
(err = CHECK(sub, dim, xdim, type, dom, iso, SUBMODEL_DEP,
type == TrendType ? ROLE_BASE : role))
== NOERROR) break;
if ((!isNegDef(type) && !isProcess(type)) || isTrend(type)) break;
// printf("trying trend %s %d %s %s %s!\n", NAME(sub), j, TYPENAMES[ type], DOMAIN_NAMES[dom], ISONAMES[iso]);
type = TrendType;
dom = XONLY;
iso = trendiso;
//
// if (j == variants -1) {
//printf("trying trend %s %d of %d %s!\n", NAME(sub), j, variants, ISONAMES[iso]);
// PMI(sub);
//
// }
}
if (err != NOERROR) {
// printf("here dddd %s %s %d of %d; %d\n", NAME(cov), NAME(sub), j, variants, TypeConsistency(type, sub, 0)); MERR(err);
return err;
}
if (cov->typus == sub->typus) {
setbackward(cov, sub);
} else {
updatepref(cov, sub);
cov->tbm2num |= sub->tbm2num;
cov->vdim[0] = sub->vdim[0];
cov->vdim[1] = sub->vdim[1];
cov->deterministic &= sub->deterministic;
};
for(j=0; j<2; j++) {
if (vdim[j] == 1) {
if (cov->vdim[j] != 1) vdim[j] = cov->vdim[j];
} else {
if (cov->vdim[j] != 1 && cov->vdim[j] != vdim[j])
SERR4("multivariate dimensionality is different in the submodels (%s is %d-variate; %s is %d-variate)", NICK(cov->sub[0]), cov->vdim[j], NICK(sub), sub->vdim[j]);
}
}
// if (vdim[1] > vdim[0]) SERR("unclear construction"); // at least currently not allowed
cov->matrix_indep_of_x &= sub->matrix_indep_of_x;
} // i, nsub
cov->vdim[0] = vdim[0];
cov->vdim[1] = vdim[1];
// !! incorrect !!
// cov->semiseparatelast = false;
//cov->separatelast = false;
return NOERROR;
}
// see private/old.code/operator.cc for some code including variable locations
void select(double *x, cov_model *cov, double *v) {
int len,
*element = PINT(SELECT_SUBNR);
cov_model *sub = cov->sub[*element];
assert(cov->vdim[0] == cov->vdim[1]);
if (*element >= cov->nsub) ERR("select: element out of range");
COV(x, sub, v);
if ( (len = cov->nrow[SELECT_SUBNR]) > 1) {
int i, m,
vdim = cov->vdim[0],
vsq = vdim * vdim;
ALLOC_EXTRA(z, vsq);
for (i=1; i<len; i++) {
sub = cov->sub[element[i]];
COV(x, sub, z);
for (m=0; m<vsq; m++) v[m] += z[m];
}
}
}
void covmatrix_select(cov_model *cov, double *v) {
int len = cov->nrow[SELECT_SUBNR];
if (len == 1) {
int element = P0INT(SELECT_SUBNR);
cov_model *next = cov->sub[element];
if (element >= cov->nsub) ERR("select: element out of range");
CovList[next->nr].covmatrix(next, v);
} else StandardCovMatrix(cov, v);
}
char iscovmatrix_select(cov_model VARIABLE_IS_NOT_USED *cov) { return 2; }
int checkselect(cov_model *cov) {
int err;
assert(cov->Splus == NULL);
if (!isCartesian(cov->isoown)) return ERRORNOTCARTESIAN;
kdefault(cov, SELECT_SUBNR, 0);
if ((err = checkplus(cov)) != NOERROR) return err;
if ((err = checkkappas(cov)) != NOERROR) return err;
EXTRA_STORAGE;
return NOERROR;
}
void rangeselect(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){
range->min[SELECT_SUBNR] = 0;
range->max[SELECT_SUBNR] = MAXSUB-1;
range->pmin[SELECT_SUBNR] = 0;
range->pmax[SELECT_SUBNR] = MAXSUB-1;
range->openmin[SELECT_SUBNR] = false;
range->openmax[SELECT_SUBNR] = false;
}
void plusStat(double *x, cov_model *cov, double *v){
cov_model *sub;
int i, m,
nsub=cov->nsub,
vdim = cov->vdim[0],
vsq = vdim * vdim;
// assert(cov->vdim[0] == cov->vdim[1]);
// printf("%d %d\n", cov->vdim[0], cov->vdim[1]);
ALLOC_EXTRA(z, vsq);
for (m=0; m<vsq; m++) v[m] = 0.0;
for(i=0; i<nsub; i++) {
sub = cov->sub[i];
// printf("PLUS %s %s ---> %s %s %d\n", NAME(sub), NICK(sub), TYPENAMES[cov->typus], TYPENAMES[sub->typus], TypeConsistency(cov->typus, sub->typus));
if (TypeConsistency(cov->typus, sub->typus)) {
// pmi(sub, 0);
//
FCTN(x, sub, z);
//printf("%s %d %s %d\n", NAME(cov), cov->typus, NAME(sub), cov->typus);
// assert(cov->typus != 9);
if (sub->vdim[0] == 1) for (m=0; m<vsq; m++) v[m] += z[0];
else for (m=0; m<vsq; m++) v[m] += z[m];
}
// if (!R_FINITE(z[0])) { PMI(sub); printf("plus i=%d m=%d x=%f %f %f\n", i, m, x[0], v[0], z[0]); }
//
}
// APMI(cov);
// assert(false);
}
void plusNonStat(double *x, double *y, cov_model *cov, double *v){
cov_model *sub;
int i, m,
nsub=cov->nsub,
vsq = cov->vdim[0] * cov->vdim[1];
assert(cov->vdim[0] == cov->vdim[1]);
ALLOC_EXTRA(z, vsq);
for (m=0; m<vsq; m++) v[m] = 0.0;
for(i=0; i<nsub; i++) {
sub = cov->sub[i];
if (TypeConsistency(cov->typus, sub->typus)) {
NONSTATCOV(x, y, sub, z);
if (sub->vdim[0] == 1) for (m=0; m<vsq; m++) v[m] += z[0];
else for (m=0; m<vsq; m++) v[m] += z[m];
}
}
}
void Dplus(double *x, cov_model *cov, double *v){
cov_model *sub;
int n = cov->nsub, i,
vsq = cov->vdim[0] * cov->vdim[1];
ALLOC_EXTRA(z, vsq);
for (int m=0; m<vsq; m++) v[m] = 0.0;
for (i=0; i<n; i++) {
sub = cov->sub[i];
if (TypeConsistency(cov->typus, sub->typus)) {
Abl1(x, sub, z);
if (sub->vdim[0] == 1) for (int m=0; m<vsq; m++) v[m] += z[0];
else for (int m=0; m<vsq; m++) v[m] += z[m];
}
}
}
void DDplus(double *x, cov_model *cov, double *v){
cov_model *sub;
int n = cov->nsub, i,
vsq = cov->vdim[0] * cov->vdim[1];
ALLOC_EXTRA(z, vsq);
for (int m=0; m<vsq; m++) v[m] = 0.0;
for (i=0; i<n; i++) {
sub = cov->sub[i];
if (TypeConsistency(cov->typus, sub->typus)) {
Abl2(x, sub, z);
if (sub->vdim[0] == 1) for (int m=0; m<vsq; m++) v[m] += z[0];
else for (int m=0; m<vsq; m++) v[m] += z[m];
}
}
}
int checkplus(cov_model *cov) {
int err, i;
if ((err = checkplusmal(cov)) != NOERROR) {
return err;
}
if (cov->domown == DOMAIN_MISMATCH) return ERRORNOVARIOGRAM;
if (cov->nsub == 0) cov->pref[SpectralTBM] = PREF_NONE;
if (isPosDef(cov) && cov->domown == XONLY) cov->logspeed = 0.0;
else if (isVariogram(cov) && cov->domown == XONLY) {
cov->logspeed = 0.0;
for (i=0; i<cov->nsub; i++) {
cov_model *sub = cov->sub[i];
if (TypeConsistency(cov->typus, sub->typus)) {
if (ISNAN(sub->logspeed)) {
cov->logspeed = RF_NA;
break;
} else cov->logspeed += sub->logspeed;
}
}
} else cov->logspeed = RF_NA;
EXTRA_STORAGE;
return NOERROR;
// spectral mit "+" funktioniert, falls alle varianzen gleich sind,
// d.h. nachfolgend DOLLAR die Varianzen gleich sind oder DOLLAR nicht
// auftritt; dann zufaellige Auswahl unter "+"
}
bool Typeplus(Types required, cov_model *cov, int depth) {
// assert(false);
bool allowed = TypeConsistency(ShapeType, required) || isTrend(required);
//||required==ProcessType ||required==GaussMethodType; not yet allowed;to do
if (!allowed) return false;
int i;
for (i=0; i<cov->nsub; i++) {
if (TypeConsistency(required, cov->sub[i], depth-1)) return true;
}
return false;
}
void spectralplus(cov_model *cov, gen_storage *s, double *e){
assert(cov->vdim[0] == 1);
int nr;
double dummy;
spec_properties *cs = &(s->spec);
double *sd_cum = cs->sub_sd_cum;
nr = cov->nsub - 1;
dummy = UNIFORM_RANDOM * sd_cum[nr];
if (ISNAN(dummy)) BUG;
while (nr > 0 && dummy <= sd_cum[nr - 1]) nr--;
cov_model *sub = cov->sub[nr];
SPECTRAL(sub, s, e); // nicht gatternr
}
int structplus(cov_model *cov, cov_model VARIABLE_IS_NOT_USED **newmodel){
int m, err;
switch(cov->role) {
case ROLE_COV : return NOERROR;
case ROLE_GAUSS :
if (isProcess(cov->typus)) {
//assert(cov->calling != NULL && (isInterface(cov->calling->typus) ||
/// isProcess(cov->calling->typus)));
assert(cov->nr != PLUS_PROC);
assert(cov->nr == PLUS);
//cov->nr = PLUS_PROC;
BUG;
//return structplusproc(cov, newmodel); // kein S-TRUCT !!
}
if (cov->Splus != NULL) BUG;
for (m=0; m<cov->nsub; m++) {
cov_model *sub = cov->sub[m];
if ((err = STRUCT(sub, newmodel)) > NOERROR) {
return err;
}
}
break;
default :
SERR2("role '%s' not allowed for '%s'", ROLENAMES[cov->role],
NICK(cov));
}
return NOERROR;
}
int initplus(cov_model *cov, gen_storage *s){
int i, err,
vdim = cov->vdim[0];
if (cov->vdim[0] != cov->vdim[1]) BUG; // ??
for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = RF_NA;
if (cov->role == ROLE_GAUSS) {
spec_properties *cs = &(s->spec);
double *sd_cum = cs->sub_sd_cum;
if (cov->vdim[0] == 1) {
for (i=0; i<cov->nsub; i++) {
cov_model *sub = cov->Splus == NULL ? cov->sub[i] : cov->Splus->keys[i];
if (sub->pref[Nothing] > PREF_NONE) { // to do ??
// for spectral plus only
COV(ZERO, sub, sd_cum + i);
if (i>0) sd_cum[i] += sd_cum[i-1];
}
cov->sub[i]->Sgen = (gen_storage *) MALLOC(sizeof(gen_storage));
if ((err = INIT(sub, cov->mpp.moments, s)) != NOERROR) {
// AERR(err);
return err;
}
sub->simu.active = true;
}
}
cov->fieldreturn = cov->Splus != NULL;
cov->origrf = false;
if (cov->Splus != NULL) cov->rf = cov->Splus->keys[0]->rf;
return NOERROR;
}
else if (cov->role == ROLE_COV) {
return NOERROR;
}
return ERRORFAILED;
}
void doplus(cov_model *cov, gen_storage *s) {
int i;
if (cov->role == ROLE_GAUSS && cov->method==SpectralTBM) {
ERR("error in doplus with spectral");
}
for (i=0; i<cov->nsub; i++) {
cov_model *sub = cov->Splus==NULL ? cov->sub[i] : cov->Splus->keys[i];
DO(sub, s);
}
}
void covmatrix_plus(cov_model *cov, double *v) {
location_type *loc = Loc(cov);
// cov_fct *C = CovList + cov->nr; // nicht gatternr
long i,
totalpoints = loc->totalpoints,
vdimtot = totalpoints * cov->vdim[0],
vdimtotSq = vdimtot * vdimtot;
int
nsub = cov->nsub;
bool is = iscovmatrix_plus(cov) >= 2;
double *mem = NULL;
if (is && nsub > 1) {
ALLOC_EXTRA2(MEM, vdimtotSq);
mem = MEM;
is = mem != NULL;
}
if (is) {
int j;
if (PisNULL(SELECT_SUBNR)) PALLOC(SELECT_SUBNR, 1, 1);
P(SELECT_SUBNR)[0] = 0;
CovList[SELECTNR].covmatrix(cov, v);
for (i=1; i<nsub; i++) {
if (Loc(cov->sub[i])->totalpoints != totalpoints) BUG;
P(SELECT_SUBNR)[0] = i;
CovList[SELECTNR].covmatrix(cov, mem);
for (j=0; j<vdimtotSq; j++) v[j] += mem[j];
}
} else StandardCovMatrix(cov, v);
}
char iscovmatrix_plus(cov_model *cov) {
char max=0, is;
int i,
nsub = cov->nsub;
for (i=0; i<nsub; i++) {
cov_model *sub = cov->sub[i];
is = CovList[sub->nr].is_covmatrix(sub);
if (is > max) max = is;
}
return max;
}
void malStat(double *x, cov_model *cov, double *v){
cov_model *sub;
int i, m,
nsub=cov->nsub,
vdim = cov->vdim[0],
vsq = vdim * vdim;
ALLOC_EXTRA(z, vsq);
// assert(x[0] >= 0.0 || cov->xdimown > 1);
for (m=0; m<vsq; m++) v[m] = 1.0;
for(i=0; i<nsub; i++) {
sub = cov->sub[i];
COV(x, sub, z);
if (sub->vdim[0] == 1) for (m=0; m<vsq; m++) v[m] *= z[0];
else for (m=0; m<vsq; m++) v[m] *= z[m];
//printf("%f %f\n", *x, *v);
}
}
void logmalStat(double *x, cov_model *cov, double *v, double *Sign){
cov_model *sub;
int i, m,
nsub=cov->nsub,
vdim = cov->vdim[0],
vsq = vdim * vdim;
ALLOC_EXTRA(z, vsq);
ALLOC_EXTRA2(zSign, vsq);
assert(cov->vdim[0] == cov->vdim[1]);
assert(x[0] >= 0.0 || cov->xdimown > 1);
for (m=0; m<vsq; m++) {v[m] = 0.0; Sign[m]=1.0;}
for(i=0; i<nsub; i++) {
sub = cov->sub[i];
LOGCOV(x, sub, z, zSign);
if (sub->vdim[0] == 1) {
for (m=0; m<vsq; m++) {
v[m] += z[0];
Sign[m] *= zSign[0];
}
} else {
for (m=0; m<vsq; m++) {
v[m] += z[m];
Sign[m] *= zSign[m];
}
}
}
}
void malNonStat(double *x, double *y, cov_model *cov, double *v){
cov_model *sub;
int i, m, nsub=cov->nsub,
vdim = cov->vdim[0],
vsq = vdim * vdim;
ALLOC_EXTRA(z, vsq);
assert(cov->vdim[0] == cov->vdim[1]);
assert(cov->vdim[0] == cov->vdim[1]);
for (m=0; m<vsq; m++) v[m] = 1.0;
for(i=0; i<nsub; i++) {
sub = cov->sub[i];
NONSTATCOV(x, y, sub, z);
if (sub->vdim[0] == 1) for (m=0; m<vsq; m++) v[m] *= z[0];
else for (m=0; m<vsq; m++) v[m] *= z[m];
}
}
void logmalNonStat(double *x, double *y, cov_model *cov, double *v,
double *Sign){
cov_model *sub;
int i, m, nsub=cov->nsub,
vdim = cov->vdim[0],
vsq = vdim * vdim;
assert(cov->vdim[0] == cov->vdim[1]);
ALLOC_EXTRA(z, vsq);
ALLOC_EXTRA2(zSign, vsq);
for (m=0; m<vsq; m++) {v[m] = 0.0; Sign[m]=1.0;}
for(i=0; i<nsub; i++) {
sub = cov->sub[i];
LOGNONSTATCOV(x, y, sub, z, zSign);
if (sub->vdim[0] == 1) {
for (m=0; m<vsq; m++) {
v[m] += z[0];
Sign[m] *= zSign[0];
}
} else {
for (m=0; m<vsq; m++) {
v[m] += z[m];
Sign[m] *= zSign[m];
}
}
}
}
void Dmal(double *x, cov_model *cov, double *v){
cov_model *sub;
int n = cov->nsub, i,
vsq = cov->vdim[0] * cov->vdim[1];
ALLOC_EXTRA3(c, vsq * MAXSUB);
ALLOC_EXTRA4(d, vsq * MAXSUB);
for (i=0; i<n; i++) {
sub = cov->sub[i];
COV(x, sub, c + i * vsq);
Abl1(x, sub, d + i* vsq);
}
*v = 0.0;
for (i=0; i<n; i++) {
double *zw = d + i *vsq;
int j;
for (j=0; j<n; j++) {
double *cc = c + j * vsq;
if (j!=i) {
for (int k=0; k<vsq; k++) zw[j] *= cc[j];
}
}
for (int k=0; k<vsq; k++) v[k] += zw[k];
}
}
int checkmal(cov_model *cov) {
cov_model *next1 = cov->sub[0];
cov_model *next2 = cov->sub[1];
int i, err,
nsub = cov->nsub;
if (next2 == NULL) next2 = next1;
// printf("cov = %s\n", NAME(next2));
if ((err = checkplusmal(cov)) != NOERROR) return err;
// printf("cov = %s no error \n", NAME(next2));
bool ok = cov->domown != DOMAIN_MISMATCH &&
(isTrend(cov->typus) ||
(isShape(cov->typus) && (!isNegDef(cov->typus) || isPosDef(cov->typus) )));
if (!ok) return ERRORNOVARIOGRAM;
// to do istcftype und allgemeiner typ zulassen
if (cov->typus == TrendType) {
ok = false;
for (i=0; i<nsub; i++)
if ((ok = cov->sub[i]->nr == CONST || cov->sub[i]->nr == BIND)) break;
if (!ok) SERR2("misuse as a trend function. At least one factor must be a constant (including 'NA') or a vector built with '%s(...)' or '%s(...).",
CovList[BIND].name, CovList[BIND].nick);
}
cov->logspeed = cov->domown == XONLY ? 0 : RF_NA;
if (cov->xdimown >= 2) cov->pref[TBM] = PREF_NONE;
if (cov->xdimown==2 && cov->nsub == 2 &&
isAnyDollar(next1) && isAnyDollar(next2)) {
double *aniso1 = PARAM(next1, DANISO),
*aniso2 = PARAM(next2, DANISO);
if (aniso1 != NULL && aniso2 != NULL) {
if (aniso1[0] == 0.0 && next1->ncol[DANISO] == 1) {
cov->pref[TBM] = next2->pref[TBM];
} else if (aniso2[0] == 0.0 && next2->ncol[DANISO] == 1) {
cov->pref[TBM] = next1->pref[TBM];
}
}
}
if (cov->ptwise_definite <= last_pt_definite) {
cov->ptwise_definite = next1->ptwise_definite;
if (cov->ptwise_definite != pt_zero) {
for (i=1; i<cov->nsub; i++) {
cov_model *sub = cov->sub[i];
if (sub->ptwise_definite == pt_zero) {
cov->ptwise_definite = pt_zero;
break;
}
if (sub->ptwise_definite != pt_posdef) {
if (sub->ptwise_definite == pt_negdef) {
cov->ptwise_definite =
cov->ptwise_definite == pt_posdef ? pt_negdef
: cov->ptwise_definite == pt_negdef ? pt_posdef
: pt_indef;
} else { // sub = indef
cov->ptwise_definite = sub->ptwise_definite;
break;
}
}
}
}
}
EXTRA_STORAGE;
return NOERROR;
}
bool Typemal(Types required, cov_model *cov, int depth) {
// printf("%d %d\n", isShape(required), isTrend(required));
if (!isShape(required) && !isTrend(required)) return false;
int i;
for (i=0; i<cov->nsub; i++) {
// print("Typemal %s %s %d\n", TYPENAMES[required], NAME(cov->sub[i]), TypeConsistency(required, cov->sub[i], depth-1));
if (!TypeConsistency(required, cov->sub[i], depth-1)) return false;
}
// print("Typemal OK\n");
return true;
}
int initmal(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){
// int err;
// return err;
return ERRORFAILED;
int i,
vdim = cov->vdim[0];
if (cov->vdim[0] != cov->vdim[1]) BUG;
for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = RF_NA;
}
void domal(cov_model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s){
BUG;
}
//////////////////////////////////////////////////////////////////////
// mpp-plus
#define PLUS_P 0 // parameter
int CheckAndSetP(cov_model *cov){
int n,
nsub = cov->nsub;
double
cump = 0.0;
if (PisNULL(PLUS_P)) {
assert(cov->nrow[PLUS_P] == 0 && cov->ncol[PLUS_P] == 0);
PALLOC(PLUS_P, nsub, 1);
for (n=0; n<nsub; n++) P(PLUS_P)[n] = 1.0 / (double) nsub;
} else {
cump = 0.0;
for(n = 0; n<nsub; n++) {
cump += P(PLUS_P)[n];
if (cump > 1.0 && n+1<nsub) return ERRORATOMP;
}
if (cump != 1.0) {
if (nsub == 1) {
warning("the p-values do not sum up to 1.\nHere only one p-value is given which must be 1.0");
P(PLUS_P)[0] = 1.0;
} else {
if (cump < 1.0 && P(PLUS_P)[nsub-1] == 0) {
WARNING1("The value of the last component of '%s' is increased.",
KNAME(PLUS_P));
}
else SERR1("The components of '%s' do not sum up to 1.", KNAME(PLUS_P));
P(PLUS_P)[nsub-1] = 1.0 - (cump - P(PLUS_P)[nsub-1]);
}
}
}
return NOERROR;
}
void kappamppplus(int i, cov_model *cov, int *nr, int *nc){
*nr = cov->nsub;
*nc = i < CovList[cov->nr].kappas ? 1 : -1;
}
void mppplus(double *x, cov_model *cov, double *v) {
int i, n,
vdim = cov->vdim[0],
vdimSq = vdim * vdim;
cov_model *sub;
ALLOC_EXTRA(z, vdimSq);
if (cov->role == ROLE_COV) {
for (i=0; i<vdimSq; i++) v[i] = 0.0;
for (n=0; n<cov->nsub; n++, sub++) {
sub = cov->sub[n];
COV(x, sub, z); // urspruenglich : covlist[sub].cov(x, cov, z); ?!
for (i=0; i<vdimSq; i++) v[i] += P(PLUS_P)[n] * z[i];
}
} else {
assert(hasPoissonRole(cov));
BUG;
}
}
int checkmppplus(cov_model *cov) {
int err,
size = 1;
cov->maxdim = MAXMPPDIM;
SERR("the current version does not support RMmppplus\n");
if ((err = checkplusmal(cov)) != NOERROR) {
return err;
}
if ((err = CheckAndSetP(cov)) != NOERROR) return(err);
if (cov->q == NULL) QALLOC(size);
EXTRA_STORAGE;
return NOERROR;
}
void rangempplus(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){
range->min[PLUS_P] = 0.0;
range->max[PLUS_P] = 1.0;
range->pmin[PLUS_P] = 0.0;
range->pmax[PLUS_P] = 1.0;
range->openmin[PLUS_P] = false;
range->openmax[PLUS_P] = false;
}
int struct_mppplus(cov_model *cov, cov_model **newmodel){
int m,
//nsub = cov->nsub,
err = NOERROR;
//if (cov->calling == NULL || cov->calling->ownloc == NULL)
// SERR("mppplus does not seem to be used in the right context");
if (!hasMaxStableRole(cov) && !hasPoissonRole(cov)) {
SERR("method is not based on Poisson point process");
}
return ERRORNOTPROGRAMMEDYET;
// Ausnahme: mppplus wird separat behandelt:
// if (nr == MPPPLUS) return S TRUCT(shape, NULL);
ASSERT_NEWMODEL_NOT_NULL;
NEW_STORAGE(plus);
plus_storage *s = cov->Splus;
for (m=0; m<cov->nsub; m++) {
cov_model *sub = cov->sub[m];
if (s->keys[m] != NULL) COV_DELETE(s->keys + m);
if ((err = covCpy(s->keys + m, sub)) != NOERROR) return err;
if ((err = addShapeFct(s->keys + m)) != NOERROR) return err;
s->keys[m]->calling = cov;
}
return NOERROR;
}
int init_mppplus(cov_model *cov, gen_storage *S) {
cov_model *sub;
double M2[MAXMPPVDIM], M2plus[MAXMPPVDIM], Eplus[MAXMPPVDIM],
maxheight[MAXMPPVDIM];
int i,n, err,
vdim = cov->vdim[0];
if (cov->vdim[0] != cov->vdim[1]) BUG;
if (vdim > MAXMPPVDIM) BUG;
ext_bool
loggiven = SUBMODEL_DEP,
fieldreturn = SUBMODEL_DEP;
pgs_storage *pgs = NULL;
for (i=0; i<vdim; i++) {
maxheight[i] = RF_NEGINF;
M2[i] = M2plus[i] = Eplus[i] = 0.0;
}
NEW_STORAGE(pgs);
pgs = cov->Spgs;
pgs->totalmass = 0.0;
for (n=0; n<cov->nsub; n++) {
sub = cov->sub[n];
//if (!sub->mpp.loc_done)
//SERR1("submodel %d of '++': the location of the shapes is not defined",
// n);
if ((err = INIT(sub, cov->mpp.moments, S)) != NOERROR) return err;
//if (!sub->mpp.loc_done) SERR("locations are not initialised");
if (n==0) {
loggiven = sub->loggiven;
fieldreturn = sub->fieldreturn;
} else {
if (loggiven != sub->loggiven) loggiven = SUBMODEL_DEP;
if (fieldreturn != sub->loggiven) fieldreturn = SUBMODEL_DEP;
}
pgs->totalmass += sub->Spgs->totalmass * P(PLUS_P)[n];
for (i=0; i<vdim; i++)
if (cov->mpp.maxheights[i] > maxheight[i])
maxheight[i] = cov->mpp.maxheights[i];
loggiven &= cov->loggiven;
// Achtung cov->mpp.mM2 und Eplus koennten nicht direkt gesetzt
// werden, da sie vom CovList[sub->nr].Init ueberschrieben werden !!
if (cov->mpp.moments >= 1) {
int nmP1 = sub->mpp.moments + 1;
for (i=0; i<vdim; i++) {
int idx = i * nmP1;
Eplus[i] += PARAM0(sub, PLUS_P) * sub->mpp.mMplus[idx + 1];
}
if (cov->mpp.moments >= 2) {
for (i=0; i<vdim; i++) {
int idx = i * nmP1;
M2[i] += PARAM0(sub, PLUS_P) * sub->mpp.mM[idx + 2];
M2plus[i] += PARAM0(sub, PLUS_P) * sub->mpp.mM[idx + 2];
}
}
}
//assert(sub->mpp.loc_done);
}
for (i=0; i<vdim; i++) cov->mpp.maxheights[i] = maxheight[i];
//cov->mpp.refsd = RF_NA;
// cov->mpp.refradius = RF_NA;
if (cov->mpp.moments >= 1) {
int nmP1 = cov->mpp.moments + 1;
for (i=0; i<vdim; i++) {
int idx = i * nmP1;
cov->mpp.mMplus[idx + 1] = Eplus[i];
cov->mpp.mM[idx + 1] = RF_NA;
}
if (cov->mpp.moments >= 2) {
for (i=0; i<vdim; i++) {
int idx = i * nmP1;
cov->mpp.mM[idx + 2] = M2[i];
cov->mpp.mMplus[idx + 2] = M2plus[i];
}
}
}
cov->loggiven = loggiven;
cov->fieldreturn = fieldreturn;
//cov->mpp.loc_done = true;
cov->origrf = false;
cov->rf = NULL;
return NOERROR;
}
void do_mppplus(cov_model *cov, gen_storage *s) {
cov_model *sub;
double subselect = UNIFORM_RANDOM;
int i, subnr,
vdim = cov->vdim[0];
assert(cov->vdim[0] == cov->vdim[1]);
for (subnr=0; (subselect -= PARAM0(cov->sub[subnr], PLUS_P)) > 0; subnr++);
cov->q[0] = (double) subnr;
sub = cov->sub[subnr];
CovList[sub->nr].Do(sub, s); // nicht gatternr
for (i=0; i<vdim; i++)
cov->mpp.maxheights[i] = sub->mpp.maxheights[i];
cov->fieldreturn = sub->fieldreturn;
cov->loggiven = sub->loggiven;
}
//////////////////////////////////////////////////////////////////////
// PROCESSES
//////////////////////////////////////////////////////////////////////
int structSproc(cov_model *cov, cov_model **newmodel) {
// printf("entering\n");
cov_model
*next = cov->sub[DOLLAR_SUB],
*Scale = cov->kappasub[DSCALE],
*Aniso = cov->kappasub[DAUSER];
int
dim = Gettimespacedim(cov),
newdim = dim,
err;
// cov_model *sub;
assert(isDollarProc(cov));
if ((Aniso != NULL && !Aniso->deterministic) ||
(Scale != NULL && !Scale->deterministic)
) {
SERR1("complicated models including arbitrary functions for '%s' cannot be simulated yet", KNAME(DAUSER));
}
switch (cov->role) {
case ROLE_GAUSS :
ASSERT_NEWMODEL_NULL;
if (cov->key != NULL) COV_DELETE(&(cov->key));
if (PrevLoc(cov)->distances)
SERR("distances do not allow for more sophisticated simulation methods");
if (Aniso!= NULL) {
//A
// crash();
TransformLoc(cov, false, true, true);
newdim = Aniso->vdim[0];
if (newdim != dim)
ERR("change of dimension in struct S not programmed yet");
location_type *loc = Loc(cov); // do not move before TransformLoc!!
int
bytes = newdim * sizeof(double);
long i,
total = loc->totalpoints;
// PMI(Aniso);
// printf("%d\n", dim);
double *v = NULL,
*x = loc->x,
*xnew = x; // not correct for newdim != dim;
// instead partial_loc_set should be used to reset loc
assert(x != NULL);
assert(!loc->grid);
if ((v = (double*) MALLOC(bytes)) == NULL) return ERRORMEMORYALLOCATION;
for (i=0; i<total; i++, x+=dim, xnew += newdim) {
FCTN(x, Aniso, v);
MEMCOPY(xnew, v, bytes);
}
FREE(v);
} else if (Scale != NULL && !isRandom(Scale)) {
SERR1("Simulation algorithms for arbitrary scale functions do not exist yet -- try using arbitrary function for '%s'", KNAME(DANISO));
} else {
//pmi(cov, 0); printf("next==intern %d\n", next->nr == TBM_PROC_INTERN);
int gridexpand = next->nr == TBM_PROC_INTERN || next->nr == NUGGET ||
next->nr == NUGGET_USER || next->nr == NUGGET_INTERN
? false : GRIDEXPAND_AVOID;
//print("here %d %d\n", (int) gridexpand, (int) GRIDEXPAND_AVOID);
TransformLocReduce(cov, true /* timesep*/, gridexpand,
true /* involveddollar */);
newdim = Gettimespacedim(cov);
//APMI(cov);
}
// APMI(cov);
if ((err = covCpy(&(cov->key), next)) != NOERROR) return err;
if (!isGaussProcess(cov->key)) addModel(&(cov->key), GAUSSPROC);
SetLoc2NewLoc(cov->key, PLoc(cov));
cov_model *key;
key = cov->key;
assert(key->calling == cov);
//dim = PrevLoc(key)->timespacedim;
if ((err = CHECK_NO_TRAFO(key, newdim, newdim, ProcessType, XONLY,
CoordinateSystemOf(cov->Sdollar->isoown),
cov->vdim[0], cov->role)) != NOERROR) {
return err;
}
err = STRUCT(key, NULL);
// MERR(err);
// APMI(key);
return err;
default :
SERR2("%s: changes in scale/variance not programmed yet for '%s'",
NICK(cov), ROLENAMES[cov->role]);
}
return NOERROR;
}
int initSproc(cov_model *cov, gen_storage *s){
// am liebsten wuerde ich hier die Koordinaten transformieren;
// zu grosser Nachteil ist dass GetDiameter nach trafo
// grid def nicht mehr ausnutzen kann -- umgehbar?!
//cov_model *next = cov->sub[DOLLAR_SUB];
//mppinfotype *info = &(s->mppinfo);
// location_type *loc = cov->prevloc;
cov_model
*key = cov->key;
//*sub = key == NULL ? next : key;
location_type *prevloc = PrevLoc(cov);
int
prevdim = prevloc->timespacedim,
dim = Gettimespacedim(cov),
err = NOERROR;
assert(key != NULL);
if ((err = INIT(key, 0, s)) != NOERROR) {
return err;
}
key->simu.active = true;
assert(s != NULL);
cov->fieldreturn = true;
// APMI(cov->calling->calling->calling->calling->calling->calling);
// printf("%d, %d %d\n",cov->ownloc != NULL, Loc(cov)->totalpoints, prevloc->totalpoints);
if ((cov->origrf = cov->ownloc != NULL && prevdim > dim)) {
// printf("hier\n");
if (cov->vdim[0] != cov->vdim[1]) BUG;
cov->rf = (double*) MALLOC(sizeof(double) *
cov->vdim[0] *
prevloc->totalpoints);
DOLLAR_STORAGE;
int d,
*proj = PINT(DPROJ),
bytes = prevdim * sizeof(int),
*cumsum = cov->Sdollar->cumsum = (int*) MALLOC(bytes),
*total = cov->Sdollar->total = (int*) MALLOC(bytes),
*len = cov->Sdollar->len = (int*) MALLOC(bytes);
cov->Sdollar->nx = (int*) MALLOC(bytes);
if (prevloc->grid) {
for (d=0; d<prevdim; d++) {
cumsum[d] = 0;
len[d] = prevloc->xgr[d][XLENGTH];
}
if (proj != NULL) {
int
nproj = cov->nrow[DPROJ];
d = 0;
assert(proj[d] > 0);
cumsum[proj[d] - 1] = 1;
for (d = 1; d < nproj; d++) {
cumsum[proj[d] - 1] = cumsum[proj[d - 1] - 1] * len[d-1];
}
} else {
int i,
iold = 0,
nrow = cov->nrow[DANISO],
ncol = cov->ncol[DANISO];
double *A = P(DANISO);
for (d=0; d<ncol; d++, A += nrow) {
for (i = 0; i < nrow && A[i] == 0.0; i++);
if (i == nrow) i = nrow - 1;
if (d > 0) {
cumsum[i] = cumsum[iold] * len[d-1];
} else { // d ==0
cumsum[i] = 1;
}
iold = i;
for (i++; i < nrow; i++) if (A[i] != 0.0) BUG; // just a check
}
}
} else { // !prevloc->grid
if (!prevloc->Time) goto Standard;
len[0] = prevloc->spatialtotalpoints;
len[1] = prevloc->T[XLENGTH];
int nproj = cov->nrow[DPROJ];
if (proj[0] != prevdim) { // spatial
for (d=1; d<nproj; d++) if (proj[d] == prevdim) goto Standard;
cumsum[0] = 1;
cumsum[1] = 0;
// pmi(cov->calling->calling->calling->calling->calling->calling, 0); APMI(cov);
} else {
if (nproj != 1) goto Standard;
cumsum[0] = 0;
cumsum[1] = 1;
//pmi(cov->calling->calling->calling->calling->calling->calling, 0);
//PMI(cov);
}
prevdim = 2;
}
for (d=0; d<prevdim; d++) total[d] = cumsum[d] * len[d];
return NOERROR;
}
Standard:
cov->origrf = false;
cov->rf = cov->key->rf;
//printf("cumsum %d\n", cov->Sdollar->cumsum[1]);
//PMI(cov);
return NOERROR;
}
void doSproc(cov_model *cov, gen_storage *s){
int i,
vdim = cov->vdim[0];
if (hasMaxStableRole(cov) || hasPoissonRole(cov)) {
assert(vdim == 1);
cov_model *next = cov->sub[DOLLAR_SUB];
cov_model
*varM = cov->kappasub[DVAR],
*scaleM = cov->kappasub[DSCALE];
assert(cov->vdim[0] == cov->vdim[1]);
if (varM != NULL && !varM->deterministic) {
assert(!PisNULL(DVAR) && isRandom(varM));
VTLG_R(NULL, varM, P(DVAR));
}
if (scaleM != NULL && !scaleM->deterministic) { // remote could be deterministic although local ist not
assert(!PisNULL(DSCALE));
VTLG_R(NULL, scaleM, P(DSCALE));
}
DO(next, s);// nicht gatternr
for (i=0; i<vdim; i++)
cov->mpp.maxheights[i] = next->mpp.maxheights[i] * P0(DVAR);
}
else if (cov->role == ROLE_GAUSS) {
assert(cov->key != NULL);
double
*res = cov->key->rf,
sd = sqrt(P0(DVAR));
int
totalpoints = Gettotalpoints(cov),
totptsvdim = totalpoints * vdim;
DO(cov->key, s);
cov_model
*varM = cov->kappasub[DVAR];
if (varM != NULL && !isRandom(varM)) {
ALLOC_NEW(Sdollar, var, totptsvdim, var);
Fctn(NULL, cov, var);
for (i=0; i<totptsvdim; i++) {
res[i] *= sqrt(var[i]);
}
} else if (sd != 1.0) {
for (i=0; i<totptsvdim; i++) res[i] *= sd;
}
}
else BUG;
if (cov->origrf) { // es wird in cov->key->rf nur
// der projezierte Teil simuliert. Dieser Teil
// muss auf das gesamte cov->rf hochgezogen werden
//if (vdim != 1) BUG;
int zaehler, d,
dim = PrevLoc(cov)->grid ? PrevLoc(cov)->timespacedim : 2,
prevtotalpts = PrevLoc(cov)->totalpoints,
owntotalpts = Loc(cov)->totalpoints,
*cumsum = cov->Sdollar->cumsum,
*nx = cov->Sdollar->nx,
*len = cov->Sdollar->len,
*total = cov->Sdollar->total;
assert(cov->key != NULL);
assert(nx != NULL && total != NULL && cumsum != NULL);
for (d=0; d<dim; d++) {
nx[d] = 0;
}
zaehler = 0;
i = 0;
for (int v=0; v<vdim; v++) {
double *res = cov->rf + v * prevtotalpts,
*rf = cov->key->rf + v * owntotalpts;
while (true) {
res[i++] = rf[zaehler];
d = 0;
nx[d]++;
zaehler += cumsum[d];
while (nx[d] >= len[d]) {
nx[d] = 0;
zaehler -= total[d];
if (++d >= dim) break;
nx[d]++;
zaehler += cumsum[d];
}
if (d >= dim) break;
}
}
}
}
int checkplusmalproc(cov_model *cov) {
cov_model *sub;
int i, err,
dim = cov->tsdim,
xdim = cov->xdimown,
role = cov->role;
Types type = ProcessType;
domain_type dom = cov->domown;
isotropy_type iso = cov->isoown;
assert(cov->Splus != NULL);
for (i=0; i<cov->nsub; i++) {
sub = cov->Splus->keys[i];
if (sub == NULL)
SERR("named submodels are not given in a sequence.");
if (!TypeConsistency(type, sub, 0)) {
return ERRORTYPECONSISTENCY;
}
if ((err= CHECK(sub, dim, xdim, type, dom, iso, SUBMODEL_DEP, role))
!= NOERROR) {
return err;
}
if (i==0) {
cov->vdim[0]=sub->vdim[0]; // to do: inkonsistent mit vorigen Zeilen !!
cov->vdim[1]=sub->vdim[1]; // to do: inkonsistent mit vorigen Zeilen !!
} else {
if (cov->vdim[0] != sub->vdim[0] || cov->vdim[1] != sub->vdim[1]) {
SERR("multivariate dimensionality must be equal in the submodels");
}
}
}
return NOERROR;
}
int checkplusproc(cov_model *cov) {
int err;
if ((err = checkplusmalproc(cov)) != NOERROR) {
return err;
}
return NOERROR;
}
int structplusmalproc(cov_model *cov, cov_model VARIABLE_IS_NOT_USED**newmodel){
int m, err,
dim = Gettimespacedim(cov);
// bool plus = cov->nr == PLUS_PROC ;
switch(cov->role) {
case ROLE_GAUSS : {
location_type *loc = Loc(cov);
NEW_STORAGE(plus);
plus_storage *s =cov->Splus;
for (m=0; m<cov->nsub; m++) {
cov_model *sub = cov->sub[m];
if (s->keys[m] != NULL) COV_DELETE(s->keys + m);
if ((err = covCpy(s->keys + m, sub)) != NOERROR) {
return err;
}
assert(s->keys[m] != NULL);
assert(s->keys[m]->calling == cov);
if (PL >= PL_STRUCTURE) {
LPRINT("plus: trying initialisation of submodel #%d (%s).\n", m+1,
NICK(sub));
}
isotropy_type newiso = UpgradeToCoordinateSystem(cov->isoprev); //!prev
if (newiso == ISO_MISMATCH) {
SERR2("'%s' for '%s' cannot be upgraded to a coordinate system",
NAME(sub), ISONAMES[cov->isoown]);
}
addModel(s->keys + m, isTrend(sub->typus) ? TRENDEVAL : GAUSSPROC);
if (isTrend(sub->typus) && sub->Spgs == NULL) {
if ((err = alloc_cov(sub, dim, sub->vdim[0], sub->vdim[1])) != NOERROR)
return err;
}
s->keys[m]->calling = cov;
if ((err = CHECK(s->keys[m], loc->timespacedim, loc->timespacedim,
ProcessType, XONLY,
newiso,
cov->vdim,
ROLE_GAUSS)) != NOERROR) {
return err;
}
if ((s->struct_err[m] =
err = STRUCT(s->keys[m], NULL)) > NOERROR) {
//PMI(s->keys[m]);
// XERR(err);
return err;
}
}
return NOERROR;
}
default :
SERR2("role '%s' not allowed for '%s'", ROLENAMES[cov->role],
NICK(cov));
}
return NOERROR;
}
int structplusproc(cov_model *cov, cov_model **newmodel){
assert(cov->nr == PLUS_PROC);
return structplusmalproc(cov, newmodel);
}
int structmultproc(cov_model *cov, cov_model **newmodel){
assert(CovList[cov->nr].check == checkmultproc);
return structplusmalproc(cov, newmodel);
}
int initplusmalproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){
int i, err,
vdim = cov->vdim[0];
assert(cov->vdim[0] == cov->vdim[1]);
bool plus = cov->nr == PLUS_PROC ;
for (i=0; i<vdim; i++)
cov->mpp.maxheights[i] = RF_NA;
if (cov->Splus == NULL) BUG;
if (cov->role == ROLE_GAUSS) {
for (i=0; i<cov->nsub; i++) {
cov_model *sub = cov->Splus == NULL ? cov->sub[i] : cov->Splus->keys[i];
if (!plus &&
(sub->nr == CONST
//|| CovList[sub[0]->nr].check == checkconstant ||
// (isDollar(sub) && CovList[sub->sub[0]->nr].check == checkconstant)
))
continue;
assert(cov->sub[i]->Sgen==NULL);
cov->sub[i]->Sgen = (gen_storage *) MALLOC(sizeof(gen_storage));
if ((err = INIT(sub, 0, cov->sub[i]->Sgen)) != NOERROR) {
return err;
}
sub->simu.active = true;
}
cov->simu.active = true;
return NOERROR;
}
else {
BUG;
}
return ERRORFAILED;
}
int initplusproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){
int err;
if ((err = initplusmalproc(cov, s)) != NOERROR) return err;
if (cov->role == ROLE_GAUSS) {
cov->fieldreturn = cov->Splus != NULL;
cov->origrf = false;
if (cov->Splus != NULL) cov->rf = cov->Splus->keys[0]->rf;
return NOERROR;
}
else {
BUG;
}
return ERRORFAILED;
}
void doplusproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
int m, i,
total = Loc(cov)->totalpoints * cov->vdim[0];
double *res = cov->rf;
assert(cov->vdim[0] == cov->vdim[1]);
if (cov->role == ROLE_GAUSS && cov->method==SpectralTBM) {
ERR("error in doplus with spectral");
}
assert(cov->Splus != NULL);
for (m=0; m<cov->nsub; m++) {
cov_model *key = cov->Splus->keys[m],
*sub = cov->sub[m];
double *keyrf = key->rf;
// PMIL(key, 2);
DO(key, sub->Sgen);
if (m > 0)
for(i=0; i<total; i++) res[i] += keyrf[i];
}
return;
}
#define MULTPROC_COPIES 0
int checkmultproc(cov_model *cov) {
int err;
kdefault(cov, MULTPROC_COPIES, GLOBAL.special.multcopies);
if ((err = checkplusmalproc(cov)) != NOERROR) {
return err;
}
EXTRA_STORAGE;
return NOERROR;
}
int initmultproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){
int err;
if ((err = initplusmalproc(cov, s)) != NOERROR) {
BUG;
return err;
}
if (cov->role == ROLE_GAUSS) {
FieldReturn(cov);
return NOERROR;
}
else {
BUG;
}
return ERRORFAILED;
}
void domultproc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s) {
double *res = cov->rf;
assert(cov->vdim[0] == cov->vdim[1]);
int m, i, c, idx,
vdim = cov->vdim[0],
// vdimSq= vdim * vdim,
total = Loc(cov)->totalpoints,
totalvdim = total * vdim,
copies = GLOBAL.special.multcopies,
factors = 0;
if (cov->role == ROLE_GAUSS && cov->method==SpectralTBM) {
ERR("error in do_mult with spectral");
}
if (cov->nsub==2 &&
((cov->sub[0]->nr==PROD) xor (idx=cov->sub[1]->nr==PROD)) &&
cov->sub[0]->nr != CONST &&
cov->sub[1]->nr != CONST) {
// koennte noch allgemeiner gemacht werden in verbindung mit CONST; todo
copies = 1;
cov->sub[idx]->Sgen->prodproc_random = false;
}
/// APMI(cov);
assert(cov->Splus != NULL);
SAVE_GAUSS_TRAFO;
for (c=0; c<copies; c++) {
for(i=0; i<totalvdim; res[i++] = 1.0);
for (m=0; m<cov->nsub; m++) {
if (PL > PL_RECURSIVE) {
PRINTF("\rcopies=%d sub=%d\n", c, m);
}
cov_model *key = cov->Splus->keys[m],
*sub = cov->sub[m];
double *keyrf = key->rf;
if (sub->nr == CONST) {
double
cc = isTrend(sub->typus) ? PARAM0(sub, CONST_C)
: sqrt(PARAM0(sub, CONST_C));
for(i=0; i<totalvdim; i++) res[i] *= cc;
}
/* else {
bool dollar = isDollar(sub);
cov_model *Const = dollar ? Const->sub[0] : sub;
if (CovList[Check->nr].check == checkconstant) {
if (dollar) {
double var = PARAM0(sub, DVAR);
bool random = false;
if (sub->kappasub[DVAR] != NULL) {
if (random = isRandom(sub->kappasub[DVAR])) {
Do(sub->kappasub[DVAR], sub->Sgen);
} else {
ALLOC_EXTRA2(VarMem, totalvdim);
F ct n(NULL, sub->kappasub[DVAR], VarMem);
}
}
if (var != 1.0) {
double sd = sqrt(var);
for(i=0; i<totalvdim; i++) res[i] *= sd;
}
} else {
ALLOC_EXTRA3(ConstMem, vdimSq);
}
} */
else {
factors ++;
DO(key, sub->Sgen);
for(i=0; i<totalvdim; i++) {
res[i] *= keyrf[i];
}
}
}
if (factors == 1) return; // no error, just exit
if (c == 0) {
ALLOC_EXTRA(z, totalvdim);
res = z;
} else {
for(i=0; i<totalvdim; i++) cov->rf[i] += res[i];
}
}
double f;
f = 1 / sqrt((double) copies);
for(i=0; i<totalvdim; i++) cov->rf[i] *= f;
}
void rangemultproc(cov_model VARIABLE_IS_NOT_USED *cov, range_type* range){
range->min[MULTPROC_COPIES] = 1.0;
range->max[MULTPROC_COPIES] = RF_INF;
range->pmin[MULTPROC_COPIES] = 1.0;
range->pmax[MULTPROC_COPIES] = 1000;
range->openmin[MULTPROC_COPIES] = false;
range->openmax[MULTPROC_COPIES] = true;
}
// $power
void PowSstat(double *x, cov_model *cov, double *v){
logPowSstat(x, cov, v, NULL);
}
void logPowSstat(double *x, cov_model *cov, double *v, double *Sign){
cov_model *next = cov->sub[POW_SUB];
double
factor,
var = P0(POWVAR),
scale =P0(POWSCALE),
p = P0(POWPOWER),
invscale = 1.0 / scale;
int i,
vdim = cov->vdim[0],
vdimSq = vdim *vdim,
xdimown = cov->xdimown;
assert(cov->vdim[0] == cov->vdim[1]);
ALLOC_DOLLAR(z, xdimown);
for (i=0; i < xdimown; i++) z[i] = invscale * x[i];
if (Sign==NULL) {
COV(z, next, v);
factor = var * pow(scale, p);
for (i=0; i<vdimSq; i++) v[i] *= factor;
} else {
LOGCOV(z, next, v, Sign);
factor = log(var) + p * log(scale);
for (i=0; i<vdimSq; i++) v[i] += factor;
}
}
void PowSnonstat(double *x, double *y, cov_model *cov, double *v){
logPowSnonstat(x, y, cov, v, NULL);
}
void logPowSnonstat(double *x, double *y, cov_model *cov, double *v,
double *Sign){
cov_model *next = cov->sub[POW_SUB];
double
factor,
var = P0(POWVAR),
scale =P0(POWSCALE),
p = P0(POWPOWER),
invscale = 1.0 / scale;
int i,
vdim = cov->vdim[0],
vdimSq = vdim * vdim,
xdimown = cov->xdimown;
assert(cov->vdim[0] == cov->vdim[1]);
ALLOC_DOLLARY(z1, z2, xdimown);
for (i=0; i<xdimown; i++) {
z1[i] = invscale * x[i];
z2[i] = invscale * y[i];
}
if (Sign==NULL) {
NONSTATCOV(z1, z2, next, v);
factor = var * pow(scale, p);
for (i=0; i<vdimSq; i++) v[i] *= factor;
} else {
LOGNONSTATCOV(z1, z2, next, v, Sign);
factor = log(var) + p * log(scale);
for (i=0; i<vdimSq; i++) v[i] += factor;
}
}
void inversePowS(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[POW_SUB];
int i,
vdim = cov->vdim[0],
vdimSq = vdim * vdim;
double y,
scale =P0(POWSCALE),
p = P0(POWPOWER),
var = P0(POWVAR);
assert(cov->vdim[0] == cov->vdim[1]);
y= *x / (var * pow(scale, p)); // inversion, so variance becomes scale
if (CovList[next->nr].inverse == ErrCov) BUG;
INVERSE(&y, next, v);
for (i=0; i<vdimSq; i++) v[i] *= scale;
}
int TaylorPowS(cov_model *cov) {
if (cov->vdim[0] != 1) SERR("Taylor only known in the unvariate case");
cov_model
*next = cov->sub[POW_SUB];
int i;
double scale = PisNULL(POWSCALE) ? 1.0 : P0(POWSCALE);
cov->taylorN = next->taylorN;
for (i=0; i<cov->taylorN; i++) {
cov->taylor[i][TaylorPow] = next->taylor[i][TaylorPow];
cov->taylor[i][TaylorConst] = next->taylor[i][TaylorConst] *
P0(POWVAR) * pow(scale, P0(POWPOWER) - next->taylor[i][TaylorPow]);
}
cov->tailN = next->tailN;
for (i=0; i<cov->tailN; i++) {
cov->tail[i][TaylorPow] = next->tail[i][TaylorPow];
cov->tail[i][TaylorExpPow] = next->tail[i][TaylorExpPow];
cov->tail[i][TaylorConst] = next->tail[i][TaylorConst] *
P0(POWVAR) * pow(scale, P0(POWPOWER) - next->tail[i][TaylorPow]);
cov->tail[i][TaylorExpConst] = next->tail[i][TaylorExpConst] *
pow(scale, -next->tail[i][TaylorExpPow]);
}
return NOERROR;
}
int checkPowS(cov_model *cov) {
// hier kommt unerwartet ein scale == nan rein ?!!
cov_model
*next = cov->sub[POW_SUB];
int err,
tsdim = cov->tsdim,
xdimown = cov->xdimown,
xdimNeu = xdimown;
if (!isCartesian(cov->isoown)) return ERRORNOTCARTESIAN;
kdefault(cov, POWVAR, 1.0);
kdefault(cov, POWSCALE, 1.0);
kdefault(cov, POWPOWER, 0.0);
if ((err = checkkappas(cov)) != NOERROR) {
return err;
}
if ((err = CHECK(next, tsdim, xdimNeu, cov->typus, cov->domown,
cov->isoown, SUBMODEL_DEP, cov->role)) != NOERROR) {
return err;
}
setbackward(cov, next);
if ((err = TaylorPowS(cov)) != NOERROR) return err;
DOLLAR_STORAGE;
return NOERROR;
}
bool TypePowS(Types required, cov_model *cov, int depth) {
cov_model *next = cov->sub[0];
return (isShape(required) || isTrend(required) || isProcess(required))
&& TypeConsistency(required, next, depth-1);
}
void rangePowS(cov_model *cov, range_type* range){
range->min[POWVAR] = 0.0;
range->max[POWVAR] = RF_INF;
range->pmin[POWVAR] = 0.0;
range->pmax[POWVAR] = 100000;
range->openmin[POWVAR] = false;
range->openmax[POWVAR] = true;
range->min[POWSCALE] = 0.0;
range->max[POWSCALE] = RF_INF;
range->pmin[POWSCALE] = 0.0001;
range->pmax[POWSCALE] = 10000;
range->openmin[POWSCALE] = true;
range->openmax[POWSCALE] = true;
range->min[POWPOWER] = RF_NEGINF;
range->max[POWPOWER] = RF_INF;
range->pmin[POWPOWER] = -cov->tsdim;
range->pmax[POWPOWER] = +cov->tsdim;
range->openmin[POWPOWER] = true;
range->openmax[POWPOWER] = true;
}
void PowScaleToLoc(cov_model *to, cov_model *from, int VARIABLE_IS_NOT_USED depth) {
assert(!PARAMisNULL(to, LOC_SCALE) && !PARAMisNULL(from, POWSCALE));
PARAM(to, LOC_SCALE)[0] = PARAM0(from, POWSCALE);
}
int structPowS(cov_model *cov, cov_model **newmodel) {
cov_model
*next = cov->sub[POW_SUB],
*scale = cov->kappasub[POWSCALE];
int err;
if (!next->deterministic) SERR("random shapes not programmed yet");
switch (cov->role) {
case ROLE_SMITH : case ROLE_GAUSS :
ASSERT_NEWMODEL_NOT_NULL;
if ((err = STRUCT(next, newmodel)) > NOERROR) return err;
addModel(newmodel, POWER_DOLLAR);
if (!PisNULL(POWVAR)) kdefault(*newmodel, POWVAR, P0(POWVAR));
if (!PisNULL(POWSCALE)) kdefault(*newmodel, POWSCALE, P0(POWSCALE));
if (!PisNULL(POWPOWER)) kdefault(*newmodel, POWPOWER, P0(POWPOWER));
break;
case ROLE_MAXSTABLE : {
ASSERT_NEWMODEL_NOT_NULL;
if ((err = STRUCT(next, newmodel)) > NOERROR) return err;
if (!isRandom(scale)) SERR("unstationary scale not allowed yet");
addModel(newmodel, LOC);
addSetDistr(newmodel, scale, PowScaleToLoc, true, MAXINT);
}
break;
default :
SERR2("'%s': changes in scale/variance not programmed yet for '%s'",
NICK(cov), ROLENAMES[cov->role]);
}
return NOERROR;
}
int initPowS(cov_model *cov, gen_storage *s){
// am liebsten wuerde ich hier die Koordinaten transformieren;
// zu grosser Nachteil ist dass GetDiameter nach trafo
// grid def nicht mehr ausnutzen kann -- umgehbar?!
cov_model *next = cov->sub[POW_SUB],
*varM = cov->kappasub[POWVAR],
*scaleM = cov->kappasub[POWSCALE];
int
vdim = cov->vdim[0],
nm = cov->mpp.moments,
nmvdim = (nm + 1) * vdim,
err = NOERROR;
bool
maxstable = hasExactMaxStableRole(cov);// Realisationsweise
assert(cov->vdim[0] == cov->vdim[1]);
assert(cov->key == NULL || ({PMI(cov);false;}));//
if (hasAnyShapeRole(cov)) { // !! ohne maxstable selbst !!
double
var[MAXMPPVDIM],
p = P0(POWPOWER),
scale = P0(POWSCALE);
int i,
intp = (int) p,
dim = cov->tsdim;
// Achtung I-NIT_RANDOM ueberschreibt mpp.* !!
if (varM != NULL) {
int nm_neu = nm == 0 && !maxstable ? 1 : nm;
if ((err = INIT_RANDOM(varM, nm_neu, s, P(POWVAR))) != NOERROR)
return err;
int nmP1 = varM->mpp.moments + 1;
for (i=0; i<vdim; i++) {
int idx = i * nmP1;
var[i] = maxstable ? P0(DVAR) : varM->mpp.mM[idx + 1];
}
} else for (i=0; i<vdim; var[i++] = P0(POWVAR));
if (scaleM != NULL) {
if (p != intp)
SERR1("random scale can be initialised only for integer values of '%s'",
KNAME(POWPOWER));
if (dim + intp < 0) SERR("negative power cannot be calculated yet");
int idx_s = maxstable ? nm : dim + nm + intp < 1 ? 1 : dim + nm + intp;
if ((err = INIT_RANDOM(scaleM, idx_s, s, P(POWSCALE))) != NOERROR) return err;
assert(scaleM->mpp.moments == 1);
scale = maxstable ? P0(DSCALE) : scaleM->mpp.mM[1];
}
if ((err = INIT(next, nm, s)) != NOERROR) return err;
for (i=0; i < nmvdim; i++) {
cov->mpp.mM[i] = next->mpp.mM[i];
cov->mpp.mMplus[i] = next->mpp.mMplus[i];
}
if (varM != NULL && !maxstable) {
int j,
nmP1 = varM->mpp.moments + 1;
for (j=0; j<vdim; j++) {
int idx = j * nmP1;
for (i=0; i <= nm; i++) {
cov->mpp.mM[i] *= varM->mpp.mM[idx + i];
cov->mpp.mMplus[i] *= varM->mpp.mMplus[idx + i];
}
}
} else {
int j, k;
double pow_var;
for (k=j=0; j<vdim; j++) {
pow_var = 1.0;
for (i=0; i<=nm; i++, pow_var *= var[j], k++) {
cov->mpp.mM[k] *= pow_var;
cov->mpp.mMplus[k] *= pow_var;
}
}
}
if (scaleM != NULL && !maxstable) {
if (dim + nm * intp < 0 || dim + intp * nm > scaleM->mpp.moments)
SERR("moments cannot be calculated");
assert(scaleM->vdim[0] == 1 && scaleM->vdim[1] == 1 );
for (i=0; i <= nm; i++) {
int idx = dim + i * intp;
cov->mpp.mM[i] *= scaleM->mpp.mM[idx];
cov->mpp.mMplus[i] *= scaleM->mpp.mMplus[idx];
}
} else {
int j,k;
double
pow_scale,
pow_s = pow(scale, dim),
pow_p = pow(scale, p);
for (k=j=0; j<vdim; j++) {
pow_scale = pow_s;
for (i=0; i <= nm; i++, pow_scale *= pow_p, k++) {
cov->mpp.mM[k] *= pow_scale;
cov->mpp.mMplus[k] *= pow_scale;
}
}
}
if (R_FINITE(cov->mpp.unnormedmass)) {
if (vdim > 1) BUG;
cov->mpp.unnormedmass = next->mpp.unnormedmass * var[0] / pow(scale, p);
} else {
double pp = pow(scale, -p);
for (i=0; i<vdim; i++)
cov->mpp.maxheights[i] = next->mpp.maxheights[i] * var[i] * pp;
}
}
else if (cov->role == ROLE_GAUSS) {
if ((err=INIT(next, 0, s)) != NOERROR) return err;
}
else if (cov->role == ROLE_BASE) {
if ((err=INIT(next, 0, s)) != NOERROR) return err;
}
else SERR("Initiation of scale and/or variance failed");
if ((err = TaylorPowS(cov)) != NOERROR) return err;
return NOERROR;
}
void doPowS(cov_model *cov, gen_storage *s){
if (hasAnyShapeRole(cov)) {
cov_model *next = cov->sub[POW_SUB];
DO(next, s);// nicht gatternr
double factor = P0(POWVAR) / pow(P0(POWSCALE), P0(POWPOWER));
int i,
vdim = cov->vdim[0];
assert(cov->vdim[0] == cov->vdim[1]);
for (i=0; i<vdim; i++)
cov->mpp.maxheights[i] = next->mpp.maxheights[i] * factor;
return;
}
BUG;
}
