https://github.com/cran/RandomFields
Tip revision: bd298816a60ec4ca975f1289dc6ad3475d6247bf authored by Martin Schlather on 09 January 2016, 13:56:44 UTC
version 3.1.8
version 3.1.8
Tip revision: bd29881
gauss.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
Definition of correlation functions and derivatives of hypermodels
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 "RF.h"
#include "Operator.h"
#include "shape_processes.h"
#include "Coordinate_systems.h"
#include "variogramAndCo.h"
void location_rules(cov_model *cov, pref_type pref) {
// 1. rules that depend on the the locations and the user's preferences only,
// but not on the covariance model#
//PMI(cov);
if (cov->nr != GAUSSPROC && cov->nr != BINARYPROC) BUG;
location_type *loc = Loc(cov);
double exactness = GLOBAL.general.exactness; //
unsigned int maxmem=500000000;
int i;
Methods Standard[Nothing] = {
CircEmbed, CircEmbedIntrinsic, CircEmbedCutoff, SpectralTBM, TBM,
Direct, Specific, Sequential, TrendEval, Average, Nugget, RandomCoin,
Hyperplane
};
for (i=0; i<Nothing; i++) {
pref[Standard[i]] = Nothing - i;
}
if (P0INT(GAUSSPROC_STATONLY) > 0) {
pref[CircEmbedIntrinsic] = LOC_PREF_NONE - 1;
}
if (!ISNA(exactness) && exactness) {
pref[TBM] = pref[SpectralTBM] = pref[Average] = pref[RandomCoin] =
pref[Sequential] = pref[Hyperplane]
= LOC_PREF_NONE - 2;
}
if (loc->timespacedim == 1) pref[TBM] -= 2 * PREF_PENALTY;
if (loc->distances) {
if (loc->grid) BUG;
for (i=CircEmbed; i<Nothing; i++)
if (i!=Direct // && i!=Nugget
) pref[i] = LOC_PREF_NONE;
} else if (loc->grid) {
if ((!ISNA(exactness) || !exactness) &&
loc->totalpoints * (1 << loc->timespacedim) * sizeof(double) > maxmem){
pref[CircEmbed] -= PREF_PENALTY;
pref[CircEmbedIntrinsic] -= PREF_PENALTY;
pref[CircEmbedCutoff] -= PREF_PENALTY;
}
} else {
if (ISNA(exactness) || !exactness){
pref[CircEmbed] -= PREF_PENALTY;
pref[CircEmbedIntrinsic] = -3;
pref[CircEmbedCutoff] -= PREF_PENALTY;
} else {
pref[CircEmbed] = pref[CircEmbedIntrinsic] = pref[CircEmbedCutoff] = -3;
}
if (!loc->Time) pref[Sequential] = LOC_PREF_NONE;
}
}
void mixed_rules(cov_model *cov, pref_type locpref,
pref_type pref, int *order) {
assert(cov->nr == GAUSSPROC);
cov_model *sub = cov->sub[0];
location_type *loc = Loc(cov);
int i,
*covpref = sub->pref,
totalpref[Nothing],
best_dirct=GLOBAL.gauss.direct_bestvariables,
max_variab= GLOBAL.direct.maxvariables,
vdim = cov->vdim[0];
for (i=0; i<Nothing; i++) {
totalpref[i] = i == Nugget ? PREF_NUGGET * (PREF_FACTOR + 1): PREF_BEST;
if (totalpref[i] > covpref[i]) {
totalpref[i] = covpref[i];
}
pref[i] = (totalpref[i] > PREF_NONE && locpref[i] > LOC_PREF_NONE)
? totalpref[i] * PREF_FACTOR + locpref[i]
: locpref[i] <= LOC_PREF_NONE ? locpref[i] : LOC_PREF_NONE - 4;
}
if (loc->totalpoints * vdim > max_variab &&
(sub->finiterange == false || GLOBAL.solve.sparse == false))
pref[Direct] = LOC_PREF_NONE - 0;
// printf("spam %f %d %f %d %d\n", GLOBAL.solve.sparse, GLOBAL.solve.sparse == false, sub->finiterange,sub->finiterange==false, pref[Direct]);
// printf("spam %f %d %d %d %d\n", GLOBAL.solve.sparse, GLOBAL.solve.sparse == false, sub->finiterange == true, pref[Direct], PREF_FACTOR);
//PMI(sub);
int vdimtot = loc->totalpoints * vdim;
if (vdimtot <= best_dirct && totalpref[Direct] == PREF_BEST)
pref[Direct] = (PREF_BEST + 1) * PREF_FACTOR;
else if (pref[Direct] >= PREF_NONE && GLOBAL.solve.sparse != true) {
#define MONTRAFO(x) ((x) * (x))
//#define MONTRAFO log
// bool orig_max = max_variab == DIRECT_ORIG_MAXVAR;
double ratio = max_variab <= DIRECT_ORIG_MAXVAR
? (double) (vdimtot - best_dirct) / (double) max_variab
: -0.1;
ratio *= fabs(ratio);
pref[Direct] = (int) pow((double) pref[Direct], 1.0 - ratio);
if (pref[Direct] < PREF_BEST)
pref[Direct] = sub->finiterange ? PREF_BEST : PREF_BEST / 2;
}
if (P0INT(GAUSSPROC_STATONLY) < 0 && isPosDef(cov)) {
pref[CircEmbedIntrinsic] = LOC_PREF_NONE - 6;
}
if (!isCartesian(cov->isoown)) {
pref[CircEmbedIntrinsic] = pref[CircEmbed] = pref[CircEmbedCutoff] =
LOC_PREF_NONE - 7;
if (isAnySpherical(cov->isoown) && cov->xdimown < 3)
pref[Sequential] = LOC_PREF_NONE - 8;
}
orderingInt(pref, Nothing, 1, order);
}
bool NAequal(double X, double Y) {
return ((ISNA(X) || ISNAN(X)) && (ISNA(Y) || ISNAN(Y))) || X == Y;
}
int kappaBoxCoxParam(cov_model *cov, int BC) {
int vdim_2 = cov->vdim[0] * 2;
if (PisNULL(BC)) {
PALLOC(BC, 2, cov->vdim[0]);
if (GLOBAL.gauss.loggauss) {
for (int i=0; i<vdim_2; i++) P(BC)[i] = 0.0;
GLOBAL.gauss.loggauss = false;
} else {
for (int i=0; i<vdim_2; i++)
P(BC)[i] = GLOBAL.gauss.boxcox[i];
}
} else {
int total = cov->nrow[BC] * cov->ncol[BC];
if (total == 1) {
double _boxcox = P0(BC);
PFREE(BC);
PALLOC(BC, 2, 1);
P(BC)[0] = _boxcox;
P(BC)[1] = GLOBAL.gauss.boxcox[1];
total = 2;
}
if (total < vdim_2) SERR("too few parameters for boxcox given");
if (total % 2 != 0) SERR("number of parameters of boxcox must even");
if (cov->nrow[BC] > 2 && cov->ncol[BC] > 1) SERR("parameters of boxcox must be given as a matrix with 2 rows");
cov->ncol[BC] = total / 2;
cov->nrow[BC] = 2;
bool notok = false;
if (GLOBAL.gauss.loggauss) {
for (int i=0; i<cov->vdim[0]; i++) {
if ((notok = P(BC)[2*i] != RF_INF &&
(P(BC)[2*i] != 0.0 || P(BC)[2*i+1] != 0.0)))
break;
}
} else {
for (int i=0; i<cov->vdim[0]; i++) {
//printf("%d %f %f %f %f\n", i, GLOBAL.gauss.boxcox[2 * i], GLOBAL.gauss.boxcox[2 * i + 1], P(BC)[2 * i], P(BC)[2 * i +1]);
if ((notok = (GLOBAL.gauss.boxcox[2 * i] != RF_INF &&
!NAequal(P(BC)[2 * i], GLOBAL.gauss.boxcox[2 * i])) ||
!NAequal(P(BC)[2 * i +1], GLOBAL.gauss.boxcox[2 * i + 1])))
break;
}
}
if (notok)
SERR("Box Cox transformation is given twice, locally and through 'RFoptions'");
}
for (int i=0; i<cov->vdim[0]; i++) {
GLOBAL.gauss.boxcox[2 * i] = RF_INF;
GLOBAL.gauss.boxcox[2 * i + 1] = 0.0;
}
return NOERROR;
}
void kappaGProc(int i, cov_model *cov, int *nr, int *nc){
*nc = i == GAUSS_BOXCOX ? SIZE_NOT_DETERMINED : 1;
*nr = i == GAUSS_BOXCOX ? SIZE_NOT_DETERMINED
: i < CovList[cov->nr].kappas ? 1 : -1;
}
int checkgaussprocess(cov_model *cov) {
cov_model
*next=cov->sub[cov->sub[0] == NULL],
*key = cov->key;
//location_type *loc=Loc(cov);
int err, role,
xdim = cov->xdimown, // could differ from tsdim in case of distances!
dim = cov->tsdim;
gauss_param *gp = &(GLOBAL.gauss);
assert((Loc(cov)->distances && xdim==1) || xdim == dim);
ROLE_ASSERT(ROLE_GAUSS || cov->role == ROLE_BERNOULLI ||
cov->role == ROLE_MAXSTABLE || cov->role == ROLE_LIKELIHOOD);
kdefault(cov, GAUSSPROC_STATONLY,
gp->stationary_only >= 0 ? gp->stationary_only : -1);
if (GLOBAL.direct.maxvariables < gp->direct_bestvariables)
SERR("maximum variables less than bestvariables for direct method");
if ((err = checkkappas(cov, false)) != NOERROR) return err;
cov->maxdim=INFDIM;
role = isVariogram(next) ? ROLE_COV
: isTrend(next) ? ROLE_GAUSS
: isGaussMethod(next) ? ROLE_GAUSS
: ROLE_UNDEFINED;
ASSERT_ROLE_DEFINED(next);
if (key == NULL) {
if (isGaussMethod(next)) {
// wird verboten
SERR1("%s may not call a method", NICK(cov));
} else {
if ((err = CHECKPD2ND(next, dim, xdim,
SymmetricOf(cov->isoown),// Jan 2015 SYMMETRIC,
SUBMODEL_DEP, role))
!= NOERROR) {
if (CHECK(next, dim, dim, TrendType, XONLY, cov->isoown, // err mocjt verwemdet
SUBMODEL_DEP, role)) return err; // previous error
}
}
// next->delflag = DEL_COV - 21;
} else {
if (PL >= PL_COV_STRUCTURE) {
PRINTF("checking key in gauss process ...\n");
}
if ((err = CHECK(key, dim, xdim, ProcessType, XONLY, cov->isoown,
SUBMODEL_DEP,
cov->role == ROLE_BASE ? ROLE_BASE : ROLE_GAUSS))
!= NOERROR) return err;
}
cov_model *sub = cov->key==NULL ? next : key;
setbackward(cov, sub);
if ((err = kappaBoxCoxParam(cov, GAUSS_BOXCOX)) != NOERROR) return err;
if ((err = checkkappas(cov, true)) != NOERROR) return err;
return NOERROR;
}
void rangegaussprocess(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range) {
GAUSS_COMMON_RANGE;
range->min[GAUSSPROC_STATONLY] = -1;
range->max[GAUSSPROC_STATONLY] = 1;
range->pmin[GAUSSPROC_STATONLY] = -1;
range->pmax[GAUSSPROC_STATONLY] = 1;
range->openmin[GAUSSPROC_STATONLY] = false;
range->openmax[GAUSSPROC_STATONLY] = false;
}
int gauss_init_settings(cov_model *cov) {
cov_model
*next = cov->sub[cov->sub[0] == NULL],
*sub = cov->key == NULL ? next : cov->key;
double sigma, var, meanDsigma,
*mean = NULL,
*variance = NULL;
int w, v, nmP1,
err = NOERROR,
vdim = next->vdim[0],
vdimSq = vdim * vdim,
vdimP1 = vdim + 1;
if ((variance = (double*) MALLOC(sizeof(double) * vdimSq)) == NULL ||
(mean = (double*) CALLOC(vdim, sizeof(double))) == NULL
) goto ErrorHandling;
GetInternalMean(next, vdim, mean);
if (ISNAN(mean[0])) // GetInternalMean currently only allows ...
GERR("Mean equals NA. Did you try a multivariate model or an incomplete (mixed) model?");
if (false && cov->key != NULL) { // sind diese Zeilen notwendig?
err = CHECK(next, cov->tsdim, cov->xdimown, PosDefType, KERNEL, SYMMETRIC,
cov->vdim, ROLE_COV);
if (err != NOERROR) goto ErrorHandling;
}
if (next->domprev == XONLY) {
COV(ZERO, next, variance);
} else {
assert(next->domprev == KERNEL);
NONSTATCOV(ZERO, ZERO, next, variance);
}
if (cov->q == NULL) QALLOC(vdim);
if ((err = alloc_mpp_M(cov, 2)) != NOERROR) goto ErrorHandling;
nmP1 = cov->mpp.moments + 1;
for (v=w=0; w<vdimSq; w+=vdimP1, v++) {
sigma = sqrt(variance[w]);
meanDsigma = sigma == 0 ? RF_INF : mean[v] / sigma;
//the following line is sqrt(2 pi s^2) int_0^infty x exp(-0.5 (x-mu)^2/s^2)
var = sigma * INVSQRTTWOPI * exp(-0.5 * meanDsigma * meanDsigma) +
mean[v] * /* correct next lines the 2nd factor is given */
pnorm(0.0, mean[v], sigma, false, false);
var = 1.0 / (var * var);
cov->q[v] = var;
// trend um factor var korrigieren in rfmodel
// eigentliches model um factor var korrigieren
// so dass ein Gausssches Zufallsfeld rauskommt, dessen Positivteil
// zumindest im Ursprung den Erwartungswert 1 hat.
cov->mpp.maxheights[v] =
GLOBAL.extreme.standardmax * sigma + ((mean[v]>0) ? mean[v] : 0);//approx!
int idx = v * nmP1;
cov->mpp.mM[idx + 0] = cov->mpp.mMplus[idx + 0] = 1.0;
cov->mpp.mMplus[idx + 1] =
sigma * INVSQRTTWOPI * exp(-0.5 * mean[v] * mean[v])
+ mean[v] * pnorm(-mean[v], 0.0, 1.0, false, false);
//(2pi)^-0.5 int x exp(-(x-m)/2s^2)
cov->mpp.mM[idx + 1] = 0.0;
cov->mpp.mM[idx + 2] = variance[w];
// todo: cov->mpp.mMplus[2] berechnen
}
cov->fieldreturn = true;
cov->origrf = false;
cov->rf = sub->rf;
ErrorHandling:
FREE(variance);
FREE(mean);
return err;
}
int struct_extractdollar(cov_model *cov, cov_model **newmodel) {
// uebernimmt struct- und init-Aufgaben !!
cov_model *next = cov->sub[0]; // nicht cov->sub[cov->sub != NULL]
location_type *loc=Loc(cov);
int role,
nr = cov->nr,
err = ERROROUTOFMETHODLIST,
xdim = cov->xdimprev, // could differ from tsdim in case of distances!
dim = cov->tsdim;
assert((loc->distances && xdim==1) || xdim == dim);
cov->fieldreturn = true;
ASSERT_NEWMODEL_NULL;
assert(cov!=NULL);
ROLE_ASSERT_GAUSS;
if ((cov->tsdim != cov->xdimprev || cov->tsdim != cov->xdimown) &&
(!loc->distances || cov->xdimprev!=1)) {
return ERRORDIM;
}
if (next!= NULL && !isVariogram(next)) {
SERR("submodel not a covariance function");
}
if (cov->key != NULL) COV_DELETE(&(cov->key));
if ((err = covCpy(&(cov->key), cov)) != NOERROR) return err; //!!cov
cov->key->nr =
nr == AVERAGE_USER ? AVERAGE_INTERN :
nr == CE_CUTOFFPROC_USER ? CE_CUTOFFPROC_INTERN :
nr == CE_INTRINPROC_USER ? CE_INTRINPROC_INTERN :
nr == HYPERPLANE_USER ? HYPERPLANE_INTERN :
nr == NUGGET_USER ? NUGGET_INTERN :
nr == RANDOMCOIN_USER ? AVERAGE_INTERN :
nr == SPECTRAL_PROC_USER ? SPECTRAL_PROC_INTERN :
nr == TBM_PROC_USER ? TBM_PROC_INTERN : MISSING_COV;
role = nr == AVERAGE_USER || nr == RANDOMCOIN_USER ? ROLE_POISSON_GAUSS
: nr == HYPERPLANE_USER ? ROLE_GAUSS : ROLE_GAUSS;
// cov!! nicht cov->key!! ??? Nein, sonst ist Coin u.ae. nicht gesetzt.
if ((err = CHECK(cov, dim, xdim, GaussMethodType, cov->domown,
cov->isoown, cov->vdim,
// role // 19.5.2013 auskommentiert und ersetzt durch
ROLE_BASE // wegen spectral.cc
)) != NOERROR) {
return err;
}
err = STRUCT(cov->key, NULL);
cov->role = ROLE_GAUSS; // wichtig !! -> check_nugget_proc
cov_model *dollar,
*firstmethod = cov->key; // z.B. coins_intern
int idx = firstmethod->sub[0] == NULL;
cov_model
*sub = firstmethod->sub[idx]; // z.B. oberste fkt shape von coins
cov_model *lastmethod= isGaussMethod(sub->typus) ? sub : firstmethod;
if (err != NOERROR && (err != ERRORPREFNONE || !isAnyDollar(sub))){
return err;
}
if (err != NOERROR) {
dollar = cov->key = lastmethod->sub[idx]; // $ jetzt im cov->key
lastmethod->sub[idx] = dollar->sub[0]; // eigentl. Modell jetzt anhaengen
lastmethod->sub[idx]->calling = lastmethod;
dollar->sub[0] = firstmethod; // auf $ setzen
firstmethod->calling = dollar;
dollar->calling = cov; // $ mit cov verbinden
dollar->prevloc = cov->prevloc;
/// cov nicht cov->key !!! ???? OK. Da sonst cov u.U. nicht gesetzt
if ((err = CHECK(cov, dim, xdim, ProcessType, cov->domown, cov->isoown,
cov->vdim, role)) != NOERROR) {
return err;
}
if ((err = STRUCT(cov->key, NULL)) != NOERROR) { return err; }
}
Methods m;
for (m=CircEmbed; m <= Hyperplane; m = (Methods) ((int) (m + 1))) {
if (gaussmethod[(Methods) m] == cov->nr) break;
}
cov->key->method = m;
return NOERROR;
}
int struct_gaussprocess(cov_model *cov, cov_model **newmodel) {
// uebernimmt struct- und init-Aufgaben !!
ASSERT_NEWMODEL_NULL;
assert(cov!=NULL);
if (cov->role == ROLE_LIKELIHOOD) {
if (cov->vdim[0] > MAXGAUSSVDIM)
SERR2("The number of variables, %d, exceeds the maximum number %d",
cov->vdim[0], MAXGAUSSVDIM);
return struct_gauss_logli(cov);
} else if (cov->role!= ROLE_GAUSS) return ERRORFAILED;
cov->fieldreturn = true;
location_type *loc = Loc(cov);
cov_model *next = cov->sub[cov->sub[0] == NULL];
pref_type locpref, pref;
int order[Nothing], i,
err = ERROROUTOFMETHODLIST;
char dummy[100 * Nothing];
int
xdim = cov->xdimprev, // could differ from tsdgetim in case of distances!
dim = cov->tsdim;
assert((loc->distances && xdim==1) || xdim == dim);
static char FailureMsg[][80] =
{"total number of points > maximum value",
"non-domain model not allowed",
"not an exact method as required by user",
"locations not on a grid",
"denied by cov. model or by user's pref.",
"no analytic solutn of Abelian integral",
"intrinsic fields denied by user"
};
bool all_PREF_NONE;
Methods unimeth = Forbidden;
int
meth,
zaehler = 0,
nr = next->nr;
#define nsel 4
if( (cov->tsdim != cov->xdimprev || cov->tsdim != cov->xdimown) &&
(!loc->distances || cov->xdimprev!=1)) {
return ERRORDIM;
}
if (cov->key != NULL) COV_DELETE(&(cov->key));
if (!isVariogram(next) && !isTrend(next)) {
SERR("submodel must be a covariance function");
}
//////////////////////////////////////////////////////////////////////
// ROLE_GAUSS; cov->key not given
{ // only for error reporting !
cov_model *sub = cov;
while(isAnyDollar(sub)) sub = sub->sub[0];
sprintf(ERROR_LOC, "Searching a simulation method for '%s': ",
NICK(sub));
}
location_rules(cov, locpref); // hier cov
mixed_rules(cov, locpref, pref, order);
if (PL >= PL_RECURSIVE) {
cov_fct *N = CovList + next->nr;
LPRINT("\n%s:%s\n", NAME(cov), NAME(next));
for (i=0; i < Nothing; i++) {
LPRINT("%-15s: base=%1d covprev=%1d locpref=%4d totpref=%6d\n",
METHODNAMES[i], N->pref[i], next->pref[i],
(locpref[i] < -9000) ? -999 : locpref[i], pref[i]);
}
LPRINT("initstandard %s (%d) [pref value %d]\n",
CovList[nr].name, nr, pref[order[Nothing -1]]);
}
all_PREF_NONE = true;
for (i=0; i<Nothing; i++) all_PREF_NONE &= pref[i] == PREF_NONE;
if (cov->key != NULL) COV_DELETE(&(cov->key));
if ((err = covCpy(&(cov->key), next)) != NOERROR) return err;
if (all_PREF_NONE) {
err = pref[Nothing] == PREF_NONE ? ERRORINVALIDMODEL : ERRORODDMODEL;
goto ErrorHandling;
}
err = ERROROUTOFMETHODLIST; // in case none of the methods are available
for (i = Nothing - 1; i>=0 && pref[order[i]] > PREF_NONE; i--) {
if (PL >= PL_BRANCHING && i < Nothing - 1) {
LPRINT("'%s%s' failed for '%s'\n",
CAT_TYPENAMES[GaussMethodType], METHODNAMES[unimeth], NICK(next));
}
zaehler++;
unimeth = (Methods) order[i];
if (PL >= PL_BRANCHING) {
LPRINT("trying '%s%s' for '%s'\n",
CAT_TYPENAMES[GaussMethodType], METHODNAMES[unimeth], NICK(next));
}
meth = gaussmethod[unimeth];
addModel(&(cov->key), meth);
cov->key->calling = cov;
cov_model *key = cov->key;
key->prevloc = PLoc(cov);
int role = cov->role != ROLE_GAUSS || key->nr != AVERAGE_INTERN
? cov->role
: ROLE_POISSON_GAUSS;
if ((err = CHECK(key, dim, xdim, GaussMethodType, cov->domown,
cov->isoown, cov->vdim, role)) == NOERROR) {
if ((err = STRUCT(key, NULL)) <= NOERROR) {
if (!key->initialised) {
if ((err = CHECK(key, dim, xdim, GaussMethodType, cov->domown,
cov->isoown, cov->vdim, role)) == NOERROR){
key->method = unimeth;
NEW_STORAGE(gen);
err = INIT(key, role == ROLE_POISSON_GAUSS ? 2 : 0, cov->Sgen);
if (err == NOERROR) {
if (PL >= PL_RECURSIVE) {
LPRINT("returning to higher level...\n");
}
goto Initialized;
}
}
}
}
if (PL >= PL_RECURSIVE) {
LPRINT("back from %s: err =%d [%d, %s]\n",
METHODNAMES[unimeth], err, unimeth, METHODNAMES[unimeth]);
if (PL > PL_ERRORS) {
char msg[LENERRMSG]; FinalErrorMSG(err, msg); PRINTF("error in init: %s\n",msg);
}
} else if (PL >= PL_BRANCHING) MERR(err) //
} else {
if (PL > PL_ERRORS) {
char msg[LENERRMSG]; FinalErrorMSG(err,msg); PRINTF("error in check: %s\n",msg);
}
}
key = key->sub[0];
COV_DELETE_WITHOUTSUB(&(cov->key));
cov->key=key;
key->calling = cov;
}
if (PL >= PL_RECURSIVE) {
strcpy(PREF_FAILURE, "");
int p, lp;
#define NMAX 14
char lpd[255], pd[255], names[NMAX];
names[NMAX-1] = '\0';
for (i=0; i<Nothing; i++) {
lp = locpref[i] - LOC_PREF_NONE;
p = pref[i] - LOC_PREF_NONE;
if (lp>0) {
strcpy(lpd, "");
} else {
sprintf(lpd, "%s (locp) ", FailureMsg[-lp]);
}
if (p>0 || p == lp) {
if (lp>0) strcpy(pd, "(method specific failure)");
else strcpy(pd, ""); // strcpy(pd, "(method specific failure)");
} else {
sprintf(pd, "%s (pref)", FailureMsg[-p]);
}
strcopyN(names, METHODNAMES[i], NMAX-1);
sprintf(dummy, "%s %-13s: %s%s\n", PREF_FAILURE, names, lpd, pd);
strcpy(PREF_FAILURE, dummy);
}
}
if (err != NOERROR) goto ErrorHandling;
Initialized:
// nachfolgende Zeile muss stehen, da Init in struct aufgerufen!
err = gauss_init_settings(cov);
cov->initialised = err == NOERROR;
ErrorHandling:
/*
if (false)
if (err > 0 && ISNA(GLOBAL.general.exactness)) {
if (*newmodel != NULL) COV_DELETE(newmodel);
GLOBAL.general.exactness = FALSE;
int newerr = struct_gaussprocess(cov, newmodel);
GLOBAL.general.exactness = RF_NA;
if (newerr == NOERROR) {
char msg[LENERRMSG];
sprintf(msg, "A somehow rougher approximation has been used. Set '%s=FALSE' to avoid this warning. See 'RFoptions' for more details.", general[GENERAL_EXACTNESS]);
warning(msg);
return NOERROR;
} else return err;
}
*/
if (err <= NOERROR || zaehler==0){
if (err <= NOERROR) {
if (cov->key == NULL) next->simu.active = true;
else cov->key->simu.active = true;
}
return err;
}
if (zaehler==1) {
sprintf(ERROR_LOC, "Only 1 method found for '%s', namely '%s', that comes into question.\nHowever", NICK(next), METHODNAMES[unimeth]);
return err;
}
{ // only for error reporting !
cov_model *sub = next;
while(isAnyDollar(sub)) sub = sub->sub[0];
sprintf(ERROR_LOC, "searching a simulation method for '%s': ",
NICK(sub));
}
return ERROROUTOFMETHODLIST;
}
int init_gaussprocess(cov_model *cov, gen_storage *s) {
/// ACHTUNG struct_gaussprocess hat bereits init aufgerufen!
// direkter Aufruf, u.a. durch RPtbm, average, cutoff, intrinsic, hyperplane,
// nugget, spectral !!!
cov_model
*key = cov->key;
int err;
ROLE_ASSERT_GAUSS;
assert(key != NULL);
if ((err = INIT(key, 0, s)) != NOERROR) return err;
if ((err = gauss_init_settings(cov)) != NOERROR) {
return err;
}
key->simu.active = true;
return NOERROR;
}
void do_gaussprocess(cov_model *cov, gen_storage *s) {
assert(s != NULL);
// reopened by internal_dogauss
char errorloc_save[nErrorLoc];
double *res = cov->rf;
int i,
vdimtot = Gettotalpoints(cov) * cov->vdim[0] ;
cov_model *key = cov->key ;
SAVE_GAUSS_TRAFO;
strcpy( errorloc_save, ERROR_LOC);
if (cov->simu.pair) {
for (i=0; i<vdimtot; i++) res[i] = -res[i];
cov->simu.pair = false;
return;
} else {
cov->simu.pair = GLOBAL.gauss.paired;
}
assert(key != NULL);
// falls gaussprocess als solcher aufgerufen, dann cov->Sgen != NULL
// falls ueber *_USER, dann i.a.(?!) nicht gesetzt, sondern s:
DO(key, cov->Sgen == NULL ? s : cov->Sgen);
// (x^\lambda_1-1)/\lambda_1+\lambda_2
BOXCOX_INVERSE;
strcpy( ERROR_LOC, errorloc_save);
}
//////////////////////////////////////////////////////////////////////
// Binary Gauss
#define BINARY_THRESHOLD (GAUSSPROC_LAST + 1)
void kappa_binaryprocess(int i, cov_model VARIABLE_IS_NOT_USED *cov,
int *nr, int *nc){
kappaGProc(i, cov, nr, nc);
if (i == BINARY_THRESHOLD) *nr = SIZE_NOT_DETERMINED;
}
int checkbinaryprocess(cov_model *cov) {
// Achtung! binary hat Doppelbedeutung: binary-Gaussian als default
// und binary angewandt auf x-bel. process.
cov_model
*key = cov->key,
*next = cov->sub[0],
*sub = key != NULL ? key : next;
double v;
int role,
err = NOERROR;
if (PisNULL(BINARY_THRESHOLD)) kdefault(cov, BINARY_THRESHOLD, 0);
// PMI(cov); printf("%d %d\n", key==NULL, isVariogram(next));
if (key == NULL && isVariogram(next)) {
if ((err = checkgaussprocess(cov)) != NOERROR) {
return err;
}
COV(ZERO, sub, &v);
if (v != 1.0)
SERR("binaryian requires a correlation function as submodel.");
} else if (isProcess(sub)) {
role = cov->role == ROLE_BASE ? cov->role : role_of_process(sub->nr);
if ((err = CHECK(sub, cov->tsdim, cov->xdimprev, ProcessType,
cov->domown, cov->isoown,
SUBMODEL_DEP,
role)) != NOERROR) return err;
setbackward(cov, sub);
} else SERR1("process type model required, but '%s' obtained", NICK(sub));
cov->vdim[0] = sub->vdim[0];
cov->vdim[1] = sub->vdim[1];
return NOERROR;
}
int struct_binaryprocess(cov_model *cov, cov_model VARIABLE_IS_NOT_USED **newmodel) {
cov_model
*next = cov->sub[0];
int err;
ROLE_ASSERT(ROLE_BERNOULLI);
if (isVariogram(next)) {
assert(cov->key == NULL);
err = covCpy(&(cov->key), cov);
// darauffolgende Zeile absichern:
if (CovList[cov->nr].kappas != 3 || CovList[GAUSSPROC].kappas != 2) BUG;
if (cov->key != NULL && !PARAMisNULL(cov->key, BINARY_THRESHOLD)) {
PARAMFREE(cov->key, BINARY_THRESHOLD);
PARAMtoNULL(cov->key, BINARY_THRESHOLD);
}
if (err != NOERROR) return err; //!!cov
cov->key->nr = GAUSSPROC;
err = CHECK(cov->key, cov->tsdim, cov->xdimprev, ProcessType,
cov->domown, cov->isoown, SUBMODEL_DEP, ROLE_GAUSS);
if (err != NOERROR) return err;
err = STRUCT(cov->key, NULL);
return err;
}
else return STRUCT(next, NULL);
}
int init_binaryprocess( cov_model *cov, gen_storage *s) {
double sigma,
*mean = NULL,
*variance = NULL,
*p = P(BINARY_THRESHOLD);
cov_model
*key = cov->key,
*next=cov->sub[0],
*sub = key != NULL ? key : next;
int v, w, nmP1, pi,
npi = cov->nrow[BINARY_THRESHOLD],
err = NOERROR,
vdim = next->vdim[0],
vdimSq = vdim * vdim,
vdimP1 = vdim + 1;
assert(next->vdim[0] == next->vdim[1]);
if ((variance = (double*) MALLOC(sizeof(double) * vdimSq)) == NULL ||
(mean = (double*) CALLOC(vdim, sizeof(double))) == NULL
) goto ErrorHandling;
if ((err = INIT(sub, 0, s)) != NOERROR) goto ErrorHandling;
cov->rf = sub->rf;
cov->origrf = false;
if (isVariogram(next) || next->nr == GAUSSPROC) {
GetInternalMean(next, vdim, mean);
if (ISNAN(mean[0])) // GetInternalMean currently only allows ...
GERR1("'%s' currently only allows scalar fields - NA returned",
NICK(cov));
if (cov->mpp.moments >= 1)
COV(ZERO, next->nr==GAUSSPROC ? next->sub[0]: next, variance);
nmP1 = cov->mpp.moments + 1;
for (pi=v=w=0; w<vdimSq; w+=vdimP1, v++, pi = (pi + 1) % npi ) {
int idx = v * nmP1;
cov->mpp.maxheights[v] = 1.0;
if (cov->mpp.moments >= 0) {
assert(cov->mpp.mM != NULL);
cov->mpp.mM[idx + 0] = cov->mpp.mMplus[idx + 0] = 1.0;
if (cov->mpp.moments >= 1) {
if (variance[w]==0.0)
GERR1("Vanishing sill not allowed in '%s'", NICK(next));
sigma = sqrt(variance[w]);
cov->mpp.mM[idx + 1] = cov->mpp.mMplus[idx + 1] =
pnorm(p[pi], mean[v], sigma, false, false);
int i;
for (i=2; i<= cov->mpp.moments; i++)
cov->mpp.mM[idx + i] = cov->mpp.mMplus[idx + i] = cov->mpp.mM[idx + 1];
}
}
}
}
cov->fieldreturn = true;
cov->simu.active = err == NOERROR;
ErrorHandling:
FREE(variance);
FREE(mean);
return err;
}
void do_binaryprocess(cov_model *cov, gen_storage *s){
// reopened by internal_dogauss
long j,
tot = Loc(cov)->totalpoints,
endfor = tot;
assert(cov->ownloc == NULL);
int i,
pi,
npi = cov->nrow[BINARY_THRESHOLD],
vdim = cov->vdim[0];
double
threshold,
*p = P(BINARY_THRESHOLD),
*rf = cov->rf;
cov_model
*next=cov->sub[0];
assert(cov->vdim[0] == cov->vdim[1]);
if (isVariogram(next)) {
do_gaussprocess(cov, s);
} else DO(next, s);
for (j=pi=i=0; i<vdim; i++, pi = (pi + 1) % npi, endfor += tot) {
threshold = p[pi];
// printf("endfor %d %f\n", endfor, threshold);
if (threshold > RF_NEGINF && threshold < RF_INF)
for ( ; j<endfor; j++) rf[j] = (double) (rf[j] >= threshold);
}
}
void rangebinaryprocess(cov_model *cov, range_type *range) {
rangegaussprocess(cov, range);
range->min[BINARY_THRESHOLD] = RF_NEGINF;
range->max[BINARY_THRESHOLD] = RF_INF;
range->pmin[BINARY_THRESHOLD] = -3;
range->pmax[BINARY_THRESHOLD] = 3;
range->openmin[BINARY_THRESHOLD] = false;
range->openmax[BINARY_THRESHOLD] = false;
}
//////////////////////////////////////////////////////////////////////
// Chisq Gauss
#define CHISQ_DEGREE 1
#define T_NU CHISQ_DEGREE
int checkchisqprocess(cov_model *cov) {
cov_model
*key = cov->key,
*next = cov->sub[0];
double *v = NULL;
int err = NOERROR,
xdim = cov->xdimprev, // could differ from tsdim in case of distances!
dim = cov->tsdim;
assert((Loc(cov)->distances && xdim==1) || xdim == dim);
if (PisNULL(CHISQ_DEGREE)) SERR("degree of freedom must be given");
if (key == NULL) {
if (!isGaussProcess(next) && !isVariogram(next))
SERR1("Gaussian process required, but '%s' obtained", NICK(cov));
if ((err = CHECK(next, dim, xdim, ProcessType, XONLY, cov->isoown,
SUBMODEL_DEP, cov->role)) != NOERROR) {
if ((err = CHECK(next, dim, xdim,
isCartesian(cov->isoown) ? VariogramType : PosDefType,
KERNEL,
SymmetricOf(cov->isoown),
SUBMODEL_DEP, ROLE_COV)) != NOERROR) {
return err;
}
}
int
vdim = next->vdim[0],
vdimSq = vdim * vdim;
assert(vdim > 0);
if ((v = (double*) MALLOC(sizeof(double) * vdimSq)) == NULL)
return ERRORMEMORYALLOCATION;
cov_model *sub = next;
while (sub != NULL && isProcess(sub)) sub=sub->sub[0];
if (sub == NULL) BUG;
if (sub->domprev == XONLY) {
COV(ZERO, next, v);
} else {
assert(sub->domprev == KERNEL);
NONSTATCOV(ZERO, ZERO, next, v);
}
int w,
vdimP1 = vdim + 1;
for (w=0; w<vdimSq; w+=vdimP1) {
if (v[w] != 1.0) {
FREE(v);
SERR("chisq requires a correlation function as submodel.");
}
}
FREE(v);
cov->vdim[0] = sub->vdim[0];
cov->vdim[1] = sub->vdim[1];
//no need to setbackward
} else {
if ((err = CHECK(key, cov->tsdim, cov->xdimprev, ProcessType,
cov->domown, cov->isoown,
SUBMODEL_DEP, cov->role)) != NOERROR) return err;
setbackward(cov, key);
}
if ((err = kappaBoxCoxParam(cov, GAUSS_BOXCOX)) != NOERROR) return err;
return NOERROR;
}
int struct_chisqprocess(cov_model *cov,
cov_model VARIABLE_IS_NOT_USED **newmodel) {
cov_model
*next = cov->sub[0];
int err;
ROLE_ASSERT_GAUSS;
if (isVariogram(next)) {
assert(cov->key == NULL);
if ((err = covCpy(&(cov->key), next)) > NOERROR) return err;
addModel(&(cov->key), GAUSSPROC);
cov->key->calling = cov;
if ((err = CHECK(cov->key, cov->tsdim, cov->tsdim, ProcessType,
cov->domprev, cov->isoprev,
SUBMODEL_DEP, ROLE_GAUSS)) != NOERROR) return err;
return STRUCT(cov->key, NULL);
} else {
return STRUCT(next, NULL);
}
}
int init_chisqprocess(cov_model *cov, gen_storage *s) {
double // sigma,
mean, m2,
variance;
cov_model
*key = cov->key,
*next = cov->sub[0],
*sub = key == NULL ? next : key;
int i, err, nmP1,
nmP1sub = sub->mpp.moments + 1,
vdim = cov->vdim[0];
cov->simu.active = false;
if ((err = INIT(sub,
CovList[cov->nr].range == rangechisqprocess ? 2 :
CovList[cov->nr].range == rangetprocess ? 1 : 9999
, s)) != NOERROR) return err;
nmP1 = cov->mpp.moments + 1;
for (i=0; i<vdim; i++) {
int
idx = i * nmP1,
idxsub = i * nmP1sub;
mean = sub->mpp.mM[idxsub + 1];
m2 = sub->mpp.mM[idxsub + 2];
variance = m2 - mean * mean;
if (variance==0.0)
SERR1("Vanishing sill not allowed in '%s'", NICK(sub));
if (ISNAN(mean)) // GetInternalMean currently only allows ...
SERR1("'%s' currently only allows scalar fields -- NA returned",
NICK(cov));
cov->mpp.maxheights[i] = GLOBAL.extreme.standardmax *
GLOBAL.extreme.standardmax * m2;
if (cov->mpp.moments >= 0) {
assert(cov->mpp.mM != NULL);
cov->mpp.mM[idx + 0] = cov->mpp.mMplus[idx + 0] = 1.0;
if (cov->mpp.moments >= 1) {
cov->mpp.mMplus[idx + 1] =
CovList[cov->nr].range == rangechisqprocess ? m2
: CovList[cov->nr].range == rangetprocess ? RF_NAN
: RF_NAN;
cov->mpp.mM[idx + 1] = RF_NA; // to do: richtiger Wert
if (cov->mpp.moments >= 2)
cov->mpp.mM[idx + 2] = 3 * variance * RF_NA;//check the 4th moment of a Gaussian, to do
}
}
}
if (CovList[cov->nr].range == rangechisqprocess) FieldReturn(cov);
else if (CovList[cov->nr].range == rangetprocess){
cov->fieldreturn = true;
cov->origrf = false;
cov->rf = sub->rf;
} else BUG;
cov->simu.active = true;
return NOERROR;
}
void do_chisqprocess(cov_model *cov, gen_storage *s){
// reopened by internal_dogauss
int i, f,
degree = P0INT(CHISQ_DEGREE),
totvdim = Loc(cov)->totalpoints * cov->vdim[0];
assert(cov->ownloc == NULL);
cov_model
*next=cov->sub[0],
*key = cov->key,
*sub = key == NULL ? next : key;
double
*subrf = sub->rf,
*res = cov->rf;
for (i=0; i<totvdim; i++) res[i] = 0.0;
// APMI(cov);
for (f=0; f<degree; f++) {
// printf("f=%d %d\n", f, degree);
DO(sub, s);
for (i=0; i<totvdim; i++) {
double dummy = subrf[i];
res[i] += dummy * dummy;
}
}
BOXCOX_INVERSE;
}
void rangechisqprocess(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range) {
GAUSS_COMMON_RANGE;
range->min[CHISQ_DEGREE] = 1;
range->max[CHISQ_DEGREE] = RF_INF;
range->pmin[CHISQ_DEGREE] = 1;
range->pmax[CHISQ_DEGREE] = 1000;
range->openmin[CHISQ_DEGREE] = false;
range->openmax[CHISQ_DEGREE] = true;
}
//////////////////////////////////////////////////////////////////////
// multivariate t distribution
void do_tprocess(cov_model *cov, gen_storage *s){
// reopened by internal_dogauss
int i,
n = Loc(cov)->totalpoints * cov->vdim[0];
cov_model
*next=cov->sub[0],
*key = cov->key,
*sub = key == NULL ? next : key;
double
nu = P0(T_NU),
factor = sqrt(nu / rgamma(0.5 * nu, 0.5)), // note inverse!! see def rgamma
*res = cov->rf;
DO(sub, s);
for (i=0; i<n; res[i++] *= factor);
BOXCOX_INVERSE;
}
void rangetprocess(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range) {
GAUSS_COMMON_RANGE;
range->min[T_NU] = 0;
range->max[T_NU] = RF_INF;
range->pmin[T_NU] = 1;
range->pmax[T_NU] = 1000;
range->openmin[T_NU] = true;
range->openmax[T_NU] = true;
}