https://github.com/cran/RandomFields
Tip revision: 41d603eb8a5f4bfe82c56acee957c79e7500bfd4 authored by Martin Schlather on 18 January 2022, 18:12:52 UTC
version 3.3.14
version 3.3.14
Tip revision: 41d603e
operator.gaussmethod.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
Definition of correlation functions and derivatives of hypermodels
Copyright (C) 2005 -- 2017 Martin Schlather
2015-2017 Olga Moreva (cutoff, modified)
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 "def.h"
#include <Basic_utils.h>
#include <R_ext/Applic.h>
#include <General_utils.h>
#include <zzz_RandomFieldsUtils.h>
#include "questions.h"
#include "operator.h"
#include "Processes.h"
#include "cubicsolver.h"
#include "extern.h"
typedef int (*set_local_q_type) (model *next, double a, double d);
void tbm3(double *x, model *cov, double *v, double tbmdim){
model *next = cov->sub[TBMOP_COV];
int i,
vdim = VDIM0,
vdimsq = vdim * vdim;
assert(VDIM0 == VDIM1);
TALLOC_XX1(v1, vdimsq);
COV(x, next, v); // x has dim 2, if turning planes
if (x[0] != 0.0) {
Abl1(x, next, v1);
for (i=0; i<vdimsq; i++) v[i] += *x * v1[i] / tbmdim;
}
END_TALLOC_XX1;
}
typedef struct TBM2_integr{
model *cov;
double *x;
} TBM2_integr;
void TBM2NumIntegrFct(double *u, int n, void *ex) {
int i;
double z[2];
TBM2_integr *info;
info = (TBM2_integr*) ex;
model *cov=info->cov;
double *x = info->x;
for (i=0; i<n; i++) {
z[0] = x[0] * SQRT(1.0 - u[i] * u[i]);
tbm3(z, cov, u + i, 1.0);
}
}
void tbm2num(double *x, model *cov, double *v){
TBM2_integr info;
#define MAXSUBDIVISIONS 100
double a = 0,
b = 1,
eps = 1e-10; // 1e-15
int maxsubdivisions = MAXSUBDIVISIONS,
lenw = 4 * MAXSUBDIVISIONS;
double abserr, work[4 * MAXSUBDIVISIONS];
int subdivisions, integralevaluations, err, iwork[MAXSUBDIVISIONS];
info.cov = cov; // a bit strange: this is tbm2, but will call tbm3...
// in TBM2NumIntegrFct ...
info.x = x;
Rdqags(TBM2NumIntegrFct, (void *) &info, &a, &b, &eps, &eps,
v, &abserr, &integralevaluations, &err,
&maxsubdivisions, &lenw, &subdivisions, iwork, work);
}
#define DO_NUMERIC 0
void tbm(double *x, model *cov, double *v){
model *next = cov->sub[TBMOP_COV];
int
fulldim = P0INT(TBMOP_FULLDIM),
tbmdim = P0INT(TBMOP_TBMDIM);
//vdimsq = cov->vdim * cov->vdim;
// if (!hasGaussMethodFrame(cov) && !hasAnyEvaluationFrame(cov)) { BUG; } else
if (fulldim == tbmdim + 2) tbm3(x, cov, v, (double) tbmdim);
else if (fulldim == 2 && tbmdim == 1) {
assert(cov->q);
if (cov->q[DO_NUMERIC]) tbm2num(x, cov, v);
else {
// APMI(next);
TBM2CALL(x, next, v);
// ASSERT_GATTER(next);
// DefList[MODELNR(next)].tbm2(x, next, v);
}
}
else XERR(ERRORTBMCOMBI);
}
void Dtbm(double *x, model *cov, double *v){
model *next = cov->sub[TBMOP_COV];
int i,
fulldim = P0INT(TBMOP_FULLDIM),
vdim = VDIM0,
vdimsq = vdim * vdim;
assert(VDIM0 == VDIM1);
TALLOC_XX1(v1, vdimsq);
double
tbmdim = (double) P0INT(TBMOP_TBMDIM),
f = 1.0 + 1.0 / tbmdim;
BUG;// to do : Taylor-Entwicklung im Ursprung
if (fulldim == tbmdim + 2) {
Abl1(x, next, v); // x has dim 2, if turning planes
Abl2(x, next, v1);
for (i=0; i<vdimsq; i++) {
v[i] = v[i] * f + *x * v1[i] / tbmdim;
}
}
else XERR(ERRORTBMCOMBI);
END_TALLOC_XX1;
}
bool numeric_tbm(model *cov) {
int n = cov->nsub;
for (int i=0; i<n; i++) if (numeric_tbm(cov->sub[i])) return true;
return !DefList[COVNR].tbm2;
}
int checktbmop(model *cov) {
// printf("entering checktbmop\n");
model *next=cov->sub[TBMOP_COV];
tbm_param *gp = &(GLOBAL.tbm);
int err;
ASSERT_ONESYSTEM;
kdefault(cov, TBMOP_FULLDIM, PisNULL(TBMOP_TBMDIM) || gp->tbmdim >= 0
? gp->fulldim : P0INT(TBMOP_TBMDIM) - gp->tbmdim);
kdefault(cov, TBMOP_TBMDIM, gp->tbmdim > 0
? gp->tbmdim : P0INT(TBMOP_FULLDIM) + gp->tbmdim);
kdefault(cov, TBMOP_LAYERS, (int) gp->layers);
if ((err = checkkappas(cov)) != NOERROR) RETURN_ERR(err);
//PMI0(cov);
if (!isVariogram(OWNTYPE(0))) {
// APMI(cov);
SERR("must be a variogram");
}
int
tbmdim = P0INT(TBMOP_TBMDIM),
fulldim = P0INT(TBMOP_FULLDIM),
vdim = VDIM0,
storedlayer = P0INT(TBMOP_LAYERS);
bool layers = storedlayer != NA_LOGICAL ? storedlayer :
OWNXDIM(0) == tbmdim + 1 && equalsSpaceIsotropic(OWN);
if(VDIM0 != VDIM1) BUG;
if (tbmdim >= fulldim)
SERR4("'%.50s' (=%d) must be less than '%.50s' (=%d)",
KNAME(TBMOP_TBMDIM), tbmdim, KNAME(TBMOP_FULLDIM), fulldim);
if (OWNLOGDIM(0) > fulldim + layers) RETURN_ERR(ERRORWRONGDIM);
if (OWNXDIM(0) > tbmdim + layers) {
SERR("dimension of coordinates does not match reduced dimension of tbm");
}
if ((err = CHECK_PASSFRAME(next, EvaluationType)) != NOERROR) {
// APMI(cov);
RETURN_ERR(err);
}
if (next->pref[TBM] == PREF_NONE) RETURN_ERR(ERRORPREFNONE);
assert(isSpaceIsotropic(OWN));
assert(isnowVariogram(cov));
assert(isXonly(OWN));
set_maxdim(OWN, 0, 0);
setbackward(cov, next);
cov->monotone=NOT_MONOTONE;
set_maxdim(OWN, 0, fulldim + layers);
cov->rese_derivs = next->rese_derivs - 1;
cov->finiterange =
(ext_bool) (((fulldim - tbmdim) % 2 == 0) && next->finiterange == true);
if (vdim > MAXTBMVDIM)
SERR2("vdim (%d) exceeds max. value of vdim in tbm3 (%d)", vdim,MAXTBMVDIM);
// only after being sure that the subsequent model does not cause
// problems. So the time dimension should be fixed.
// This is not absolutely safe programming, but should be OK.
// But difficult to get around for MLE calls that do not allow
// for NAs values in integer variables.
PINT(TBMOP_LAYERS)[0] = layers;
if (fulldim == 2 && tbmdim == 1) QALLOC1(numeric_tbm(cov));
EXTRA_STORAGE;
RETURN_NOERROR;
}
Types Typetbm(Types required, model *cov, isotropy_type i) {
bool layers = P0INT(TBMOP_LAYERS);
assert(cov->sub[0] != NULL);
// printf("i=%d %.50s\n",i, ISO_NAMES[i]);
if (!isCartesian(i) ||
((OWNXDIM(0) == 1) xor equalsIsotropic(i)) ||
((OWNXDIM(0) == 2) xor equalsSpaceIsotropic(i)) ||
(OWNXDIM(0) > 2) ||
(layers != NA_INTEGER && layers && !equalsSpaceIsotropic(i)) ||
!equalsXonly(OWNDOM(0))) return BadType;
return TypeConsistency(required, cov->sub[0], i);
}
bool settbm(model *cov) {
isotropy_type iso = CONDPREVISO(0);
if (!isFixed(iso)) return false;
tbm_param *gp = &(GLOBAL.tbm);
// critical parameter TBMOP_LAYER is already treated by SUBMODEL_I
kdefault(cov, TBMOP_LAYERS, (int) gp->layers);
set_type(OWN, 0, PREVTYPE(0));
set_iso(OWN, 0, P0INT(TBMOP_LAYERS) ? DOUBLEISOTROPIC : ISOTROPIC);
return true; // da nur PREVISO kritisch ist
}
bool allowedItbm(model *cov) {
tbm_param *gp = &(GLOBAL.tbm);
bool *I = cov->allowedI;
kdefault(cov, TBMOP_LAYERS, (int) gp->layers);
for (int i=(int) FIRST_ISOUSER; i<=(int) LAST_ISOUSER; I[i++] = false);
I[P0INT(TBMOP_LAYERS) ? DOUBLEISOTROPIC : ISOTROPIC] = true;
return false;
}
void rangetbm_common(model VARIABLE_IS_NOT_USED *cov, range_type *range, bool tbmop ){
int
TBMDIM = tbmop ? TBMOP_TBMDIM : TBM_TBMDIM,
FULLDIM = tbmop ? TBMOP_FULLDIM : TBM_FULLDIM,
LAYERS = tbmop ? TBMOP_LAYERS : TBM_LAYERS;
range->min[FULLDIM] = 1.0;
range->max[FULLDIM] = RF_INF;
range->pmin[FULLDIM] = 1.0;
range->pmax[FULLDIM] = 100;
range->openmin[FULLDIM] = false;
range->openmax[FULLDIM] = true;
range->min[TBMDIM] = RF_NEGINF;
range->max[TBMDIM] = RF_INF;
range->pmin[TBMDIM] = RF_NEGINF;
range->pmax[TBMDIM] = 100;
range->openmin[TBMDIM] = false;
range->openmax[TBMDIM] = true;
booleanRange(LAYERS);
}
void rangetbmop(model *cov, range_type *range){
rangetbm_common(cov, range, true);
}
int set_stein_q(model *cov, double r, double d) {
localCE_storage *s = cov->calling->SlocalCE;
double C0, phi0, phi1, phi2,
zero = 0.0,
rP1 = r + 1.0,
rM1 = r - 1.0,
dsq = d * d;
localvariab *q = s->q2;
COV(&zero, cov, &C0);
COV(&d, cov, &phi0);
Abl1(&d, cov, &phi1); // derivative of
phi1 *= d; // scaled fctn at 1
Abl2(&d, cov, &phi2); // 2nd derivative of
phi2 *= dsq; // scaled fctn at 1
q->R = r * d;
double dummy = (phi2 - phi1) / (3.0 * r * rP1);
q->intrinsic.B = (r == 1.0) ? 0.0 : dummy / (rM1 * dsq);
q->intrinsic.A2 = (dummy - phi1 / 3.0 - phi2 / 6.0) / dsq;
q->intrinsic.A0 = 0.5 * rM1 / rP1 * phi2 + phi1 / rP1 - phi0;
if ((q->intrinsic.B < 0.0) || (q->intrinsic.A2 < 0.0) ||
(q->intrinsic.A0 + C0 < 0.0)) {
return MSGLOCAL_INITINTRINSIC;
}
RETURN_NOERROR;
}
void getlocalconstants(double phi0, double phi1, double phi2, double phi3,
double radius, double m, double n, double l,
double *aa, double *bb, double *cc, double *constant) {
double phi3r2 = phi3*radius*radius,
phi2r = phi2 * radius;
*aa = -(phi3r2 + phi2r * (m + l - 3) + phi1 *
(m - 1) * (l - 1)) / (n * (m - n) * (l - n) * POW(radius, n - 1));
*bb = -(phi3r2 + phi2r * (n + l - 3) + phi1 *
(n - 1) * (l - 1)) / (m * (n - m)*(l - m) * POW(radius, m - 1));
*cc = -(phi3r2 + phi2r*(m + n - 3) + phi1*
(m - 1)*(n - 1)) / (l*(m - l)*(n - l) * POW(radius, l - 1));
*constant = -phi0 + *aa * POW(radius, n) +
*bb * POW(radius, m) + *cc * POW(radius, l);
}
int set_cutoff_q2variate(model *cov, double a, double d) {
// auf modell ebene, d.h. co->sub[0] oder stein->sub[0]
double
phi0[4], phi1[4], phi2[4],
a2 = a * a;
// radius = MAXINT;
//double roots[3][2];
assert(cov->calling->SlocalCE != NULL);
localCE_storage *s = cov->calling->SlocalCE;
localvariab *q = s->q2;
COV(&d, cov, phi0);
// COV(&d, cov->sub[0], &phi01);
Abl1(&d, cov, phi1);
Abl2(&d, cov, phi2);
//bivariate Whittle Cauchy Stable
//assert(false);
if (VDIM0 > LOCALCE_MAXVDIM || VDIM1 > LOCALCE_MAXVDIM) BUG;
if (phi1[1] != phi1[2] )
return MSGLOCAL_NOTSYMMETRICMULTIVARIATE;
for (int j = 0; j < 4; j++) {
q[j].a = a;
if ( phi0[j] > 0 ) {
q[j].cutoff.b = phi2[j]*phi2[j]*phi2[j]/(108.0*phi1[j]*phi1[j]);
q[j].cutoff.theor = POW(1.0 - 3.0 * a2 * phi1[j] / (d* phi2[j]), 1.0 / a);
q[j].R = d * q[j].cutoff.theor;
q[j].cutoff.asqrtr = POW(q[j].R, a);
q[j].cutoff.constant = phi0[j] -q[j].cutoff.b*POW(-3.0*phi1[j]/phi2[j], 4);
} else {
q[j].cutoff.b = 0;
q[j].cutoff.theor = 0;
q[j].R = 0;
q[j].cutoff.asqrtr = 0;
q[j].cutoff.constant = 0;
}
}
if (q[1].R > q[0].R || q[1].R > q[3].R ) {
// printf("%10g > %10g | %10g > %10g\n", q[1].R, q[0].R, q[1].R, q[3].R);
// APMI(cov->calling);
return MSGLOCAL_WRONGRADII;
}
//PMI(cov);
RETURN_NOERROR;
}
int set_cutoff_q(model *cov, double a, double d) {
assert(VDIM0 > 0 && VDIM0 <= LOCALCE_MAXVDIM);
if (VDIM0 > 1) return set_cutoff_q2variate(cov, a, d);
localCE_storage *s = cov->calling->SlocalCE;
// auf modell ebene, d.h. co->sub[0] oder stein->sub[0]
double
phi0, phi1,
phi2 = RF_NA,
phi3 = RF_NA,
phi4 = RF_NA,
a2 = a * a,
// d2 = d*d,
radius = UNSET;
localvariab *q = s->q2;
q->a = a;
// if ((err = covcpy(&neu, cov)) != NOERROR) RETURN_ERR(err);
// DefList[cov->gatternr].cov(&d, neu, &phi0);
// DefList[cov->gatternr].D(&d, neu, &phi1);
//FREE(neu);
COV(&d, cov, &phi0);
Abl1(&d, cov, &phi1);
//If it is a variogram, then we must determine a constant such that Const - Variogram is a valied covariance
// The constant depends on a and d
// cov->calling->typus
// if (cov->typus == V ariogramType) {
if (equalsnowVariogram(cov)) {
if (a == 0.5) {
//COV(&d, cov, &q[CUTOFF_CONSTANT]);
//q[CUTOFF_B] = -2*phi1*SQRT(d);
//q[CUTOFF_THEOR] = POW(1.0 - 0.5* (q[CUTOFF_CONSTANT] + phi0) / phi1 /d, 1/a);
//q[LOCAL_R] = d * q[CUTOFF_THEOR];
q->cutoff.constant = -phi0;
q->R = d;
q->cutoff.asqrtr = POW(q->R, a);
} else if (a == 1 ) {
//second derivative
Abl2(&d, cov, &phi2);
if (phi2 <= 0 || phi1 >= 0 || phi0 >= 0 ) {
return MSGLOCAL_SIGNPHISND;
}
//if second derivative is positive at d, then
//q[CUTOFF_CONSTANT] = -phi0+phi1*phi1/(2*phi2) + EPSILON_C;
//q[CUTOFF_B] = 0.25*phi1*phi1/(q[CUTOFF_CONSTANT] + phi0);
//q[CUTOFF_THEOR] = POW(1.0 - 2 * ( q[CUTOFF_CONSTANT] + phi0) / phi1 /d, 1/a);
q->cutoff.constant = -phi0+phi1*phi1/(2*phi2);
q->cutoff.b = phi2/2;
q->R = d - phi1/phi2;
q->cutoff.asqrtr = POW(q->R, a);
} else if (a == CUTOFF_THIRD_CONDITION) {
Abl2(&d, cov, &phi2);
Abl3(&d, cov, &phi3);
Abl4(&d, cov, &phi4);
double roots[3][2];
//de Wijsian model
int
n = 4,
m = 6,
l = 7;
q->cube.N = n;
q->cube.M = m;
q->cube.L = l;
cubicsolver( phi4, (n+m+l-6)*phi3,
(n*m + n*l + m*l - 3*(n+m+l) + 7)*phi2,
(n-1)*(m-1)*(l-1)*phi1 ,
roots);
//find a positive root of the cubic polynomial
for (int i = 0; i < 3; i++) {
if (roots[i][1] == 0 && roots[i][0] > radius ) {
radius = roots[i][0];
}
}
//if no positive root, cannot do anything
if (radius <= 0.0) return MSGLOCAL_NOPOSITIVEROOT;
getlocalconstants(phi0, phi1, phi2, phi3, radius, m, n, l,
&(q->cube.A), &(q->cube.B), &(q->cube.C),
&(q->cube.constant));
if (q->cube.constant <= 0.0) return MSGLOCAL_SIGNPHI;
q->R = radius + d;
// q->cutoff.asqrtr = radius; Olga, line looks strange
} else BUG;
} else {
//now we are in a positive definite function case
if (phi0 <= 0.0) return MSGLOCAL_SIGNPHI;
if (phi1 >= 0.0) return MSGLOCAL_SIGNPHIFST;
//find parameters of the method
if (a != CUTOFF_THIRD_CONDITION) {
phi1 *= d;
if (phi1 >= 0.0) return MSGLOCAL_SIGNPHIFST;
q->cutoff.b =//sound qconstant even if variance of submodel is not 1
POW(-phi1 / (2.0 * a2 * phi0), 2.0 * a) * phi0 / POW(d, 2.0 * a2);
// assert(false);
q->cutoff.theor = POW(1.0 - 2.0 * a2 * phi0 / phi1, 1.0 / a);
q->R = d * q->cutoff.theor;
q->cutoff.asqrtr = POW(q->R, a);
} else if (a == CUTOFF_THIRD_CONDITION){
//a = CUTOFF_THIRD_CONDITION
Abl2(&d, cov, &phi2);
Abl3(&d, cov, &phi3);
Abl4(&d, cov, &phi4);
double roots[3][2];
//Whittle model
int
n = 5,
m = 6,
l = 7;
q->cube.N = n;
q->cube.M = m;
q->cube.L = l;
cubicsolver( phi4, (n+m+l-6)*phi3,
(n*m + n*l + m*l - 3*(n+m-l) + 7)*phi2,
(n-1)*(m-1)*(l-1)*phi1 ,
roots);
//find a positive root of the cubic polynomial
for (int i = 0; i < 3; i++) {
if (roots[i][1] == 0 && roots[i][0] > radius ) {
radius = roots[i][0];
}
}
//if no positive root, cannot do anything
if (radius <= 0.0) return MSGLOCAL_NOPOSITIVEROOT;
getlocalconstants(phi0, phi1, phi2, phi3, radius, m, n, l,
&(q[0].cube.A), &(q[0].cube.B), &(q[0].cube.C),
&(q[0].cube.constant));
if (q->cube.constant <= 0.0) return MSGLOCAL_SIGNPHI;
q->R = radius + d;
} else BUG;
}
RETURN_NOERROR;
}
int check_local(model *cov,
Methods method,
getlocalparam init, // next->coinit
set_local_q_type set_local) {
// PMI(cov);
localCE_storage *s = cov->SlocalCE;
assert(s != NULL);
location_type *loc = Loc(cov);
int i, msg,
dim = OWNLOGDIM(0),
err=NOERROR;
localvariab
*q = s->q,
*q2 = s->q2;
double d=RF_NA;
model *next = cov->sub[0];
localinfotype li;
//ce_param *gp = &(GLOBAL.localce); // ok
//printf("\n \n entering check local from %d:%.50s\n", cov->CALLINGNR,
//DefList[cov->CALLINGNR].name);
if ((err = CHECK(next, dim, 1,
method == CircEmbedCutoff ? PosDefType : VariogramType,
OWNDOM(0), OWNISO(0), SUBMODEL_DEP, EvaluationType)) != NOERROR) {
// printf("\n \n \n ::: checl_local :: if ::: dim = %d --- next->vdim[0] = %d::: \n \n \n", dim, next->vdim[0]);
if ( method != CircEmbedCutoff ||
(err = CHECK(next, dim, 1,
VariogramType,
OWNDOM(0), OWNISO(0), SUBMODEL_DEP, EvaluationType)) != NOERROR)
RETURN_ERR(err);
}
if (VDIM0 > LOCALCE_MAXVDIM) SERR("vdim of submodel must be less than 3");
// no setbackward ?!
setbackward(cov, next);
if (next->pref[method] == PREF_NONE) RETURN_ERR(ERRORPREFNONE);
if (init == NULL) RETURN_ERR(ERRORUNKNOWNMETHOD);
if (PisNULL(pLOC_DIAM)) {
double
diameter = GetDiameter(loc);
if (PL>=PL_DETAILS) { LPRINT("diameter %10g\n", diameter); }
kdefault(cov, pLOC_DIAM, diameter);
} else {
d = P0(pLOC_DIAM);
}
if (PisNULL(pLOC_A)) {
if (DefList[NEXTNR].implemented[method]) {
assert(init != NULL);
init(next, &li);
if (li.instances == 0) {
SERR("parameter values do not allow for finding second parameter");
}
q->R = RF_INF;
msg = XMSGLOCAL_FAILED;
for (i=0; i<li.instances; i++) {
// printf("li %d %d %d %10g %d\n", i, li.instances, li.msg[i], li.value[i], msg);
if (li.msg[i] <= msg) {
//
err = set_local(next, li.value[i], d);
//
if (err == NOERROR && (li.msg[i] < msg || q2->R < q->R)){
MEMCOPY(q, q2, sizeof(localvariabArray));
msg = li.msg[i];
if (PisNULL(pLOC_A)) {
PALLOC(pLOC_A, 1, 1);
}
P(pLOC_A)[0] = li.value[i]; // co: a; ie: rawr
}
}
}
q->msg = msg;
if (msg == MSGLOCAL_FAILED) {
// APMI(cov);
SERR2("'%.50s' did not work out correctly. Likely, invalid parameter \
values have been chosen for '%.50s'", NICK(cov), NICK(next));
}
} else {
SERR("2nd parameter is neither given nor can be found automatically");
}
} else {
if (cov->ncol[pLOC_A] != 1 || cov->nrow[pLOC_A] != 1)
SERR1("'%.50s' must be a scale", KNAME(pLOC_A));
err = set_local(next, P0(pLOC_A), d);
MEMCOPY(q, q2, sizeof(localvariabArray));
}
cov->pref[CircEmbed] = 5;
RETURN_ERR(err);
}
void co(double *x, model *cov, double *v) {
model *next = cov->sub[0];
localCE_storage *s = cov->SlocalCE;
localvariab *q = s->q;
double y=*x,
diameter = P0(pLOC_DIAM),
a = P0(pLOC_A);
// assert(hasGaussMethodFrame(cov) || hasAnyEvaluationFrame(cov));
//VDIM0 or next->vdim[0]?
//if ( cov->SlocalCE->is_multivariate_cutoff == true ) {
if ( VDIM0 > 1 ) {
if (y <= diameter) {
COV(x, next, v);
for (int i = 0; i < 4; i++) {
v[i] = v[i] - q[i].cutoff.constant;
}
} else {
for (int i = 0; i < 4; i++) {
if (y >= q[i].R) v[i] = 0.0;
else v[i] = q[i].cutoff.b*POW(q[i].cutoff.asqrtr - POW(y, a), 4.0 * a);
}
}
}
else {
if (y <= diameter) {
//If it is a variogram (actually -variogram), add a constant. If it is positive definite, go below
if (isnowVariogram(next)) {
COV(x, next, v);
*v = *v + q->cutoff.constant;
} else
{
COV(x, next, v)
}
}
else {
if (y >= q->R) *v = 0.0;
else {
if (a != CUTOFF_THIRD_CONDITION)
*v = q->cutoff.b * POW(q->cutoff.asqrtr - POW(y, a), 2.0 * a);
else
*v = q->cube.A * POW((q->R - y), q->cube.N) +
q->cube.B * POW((q->R - y), q->cube.M) +
q->cube.C * POW((q->R - y), q->cube.L);
//BUG;
}
}
}
}
int check_co(model *cov) {
model *next = cov->sub[0];
int err;
NEW_STORAGE(localCE);
assert(cov->SlocalCE != NULL);
// 1 ? 2? More than 2? or inside set_cutoff_q2variate
//PMI(cov);assert(next->vdim[0]>0);
err = check_local(cov, CircEmbedCutoff,
DefList[NEXTNR].coinit, set_cutoff_q);
RETURN_ERR(err);
}
bool alternativeparam_co(model VARIABLE_IS_NOT_USED *cov){
return false;
}
void range_co(model VARIABLE_IS_NOT_USED *cov, range_type *range){
range->min[pLOC_DIAM] = 0.0; // CUTOFF_DIAM
range->max[pLOC_DIAM] = RF_INF;
range->pmin[pLOC_DIAM] = 1e-10;
range->pmax[pLOC_DIAM] = 1e10;
range->openmin[pLOC_DIAM] = true;
range->openmax[pLOC_DIAM] = true;
range->min[pLOC_A] = 0.0; // cutoff_a
range->max[pLOC_A] = RF_INF;
range->pmin[pLOC_A] = 0.5;
range->pmax[pLOC_A] = 2.0;
range->openmin[pLOC_A] = true;
range->openmax[pLOC_A] = true;
}
void Stein(double *x, model *cov, double *v) {
model *next = cov->sub[0];
localCE_storage *s = cov->SlocalCE;
localvariab *q = s->q;
double y=*x, z,
diameter = P0(pLOC_DIAM);
// assert(hasAnyEvaluationFrame(cov) || hasGaussMethodFrame(cov));
if (y <= diameter) {
COV(x, next, v);
*v += q->intrinsic.A0 + q->intrinsic.A2 * y * y;
} else {
z = q->R - y;
if (z <= 0.0) *v = 0.0;
else *v = q->intrinsic.B * z * z * z / y;
}
}
int check_Stein(model *cov)
{
model *next = cov->sub[0]; // dito
NEW_STORAGE(localCE);
assert(cov->SlocalCE != NULL);
return check_local(cov, CircEmbedIntrinsic, DefList[NEXTNR].ieinit,
set_stein_q);
}
bool alternativeparam_Stein(model *cov) {
P(pLOC_A)[0] *= 2.0;
return true;
}
void range_Stein(model VARIABLE_IS_NOT_USED *cov, range_type *range) {
range->min[pLOC_DIAM] = 0.0;
range->max[pLOC_DIAM] = RF_INF;
range->pmin[pLOC_DIAM] = 0.01;
range->pmax[pLOC_DIAM] = 100;
range->openmin[pLOC_DIAM] = true;
range->openmax[pLOC_DIAM] = true;
range->min[pLOC_R] = 1.0; // stein_r
range->max[pLOC_R] = RF_INF;
range->pmin[pLOC_R] = 1.0;
range->pmax[pLOC_R] = 20.0;
range->openmin[pLOC_R] = false;
range->openmax[pLOC_R] = true;
}