https://github.com/cran/RandomFields
Tip revision: af7aa2881d9fd941be5411589f41ac2ed9d24e57 authored by Martin Schlather on 11 March 2011, 00:00:00 UTC
version 2.0.45
version 2.0.45
Tip revision: af7aa28
operator.cc
/*
Authors
Martin Schlather, martin.schlather@math.uni-goettingen.de
Definition of correlation functions and derivatives of hypermodels
Copyright (C) 2005 -- 2011 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 2
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 <assert.h>
#include "RF.h"
#include "Covariance.h"
#include <R_ext/Lapack.h>
#include <R_ext/Applic.h>
// ------- end internal functions
void co(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[0];
double y=*x, *q=cov->q,
diameter = cov->p[pLOC_DIAM][0],
a = cov->p[pLOC_A][0];
if (y <= diameter) CovList[next->nr].cov(x, next, v);
else {
*v = (y >= q[LOCAL_R]) ? 0.0
: q[CUTOFF_B] * pow(q[CUTOFF_ASQRTR] - pow(y, a), 2.0 * a);
}
}
bool alternativeparam_co(cov_model *cov){
return false;
}
void range_co(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = 0.0; // CUTOFF_DIAM
range->max[0] = RF_INF;
range->pmin[0] = 1e-10;
range->pmax[0] = 1e10;
range->openmin[0] = true;
range->openmax[0] = true;
range->min[1] = 0.0; // cutoff_a
range->max[1] = RF_INF;
range->pmin[1] = 0.5;
range->pmax[1] = 2.0;
range->openmin[1] = true;
range->openmax[1] = true;
}
int check_local(cov_model *cov,
SimulationType type, int maxq,
getlocalparam init, // next->coinit
set_local_q_type set_local) {
int i, err=NOERROR, msg; // timespacedim -- needed ?;
double *q, q2[LOCAL_MAX],
d=RF_NAN, **p = cov->p;
cov_model *next = cov->sub[0],
*truenext = next->sub[0];
localinfotype li;
cov_fct *C;
cov->qlen = maxq;
// checkkappas(cov) is not called, but checked directly !!
if ((err = check2X(next, cov->tsdim, 1,
type == CircEmbedCutoff ? STATIONARY : VARIOGRAM,
ISOTROPIC, UNIVARIATE)) != NOERROR)
return err;
next->tsdim = DEL_COV;
// setbackward(cov, next);
if (next->maxdim < cov->maxdim) cov->maxdim = next->maxdim;
cov->manipulating_x=true;
cov->normalmix = false;
// cov->finiterange = true by default
cov->diag &= next->diag;
// cov->quasidiag &= next->quasidiag;
cov->separatelast &= next->separatelast;
cov->semiseparatelast &= next->semiseparatelast;
// cov->variogram = false;
// no changing in pref by submodel !!
if (cov->q != NULL) {
free(cov->q);
// PrintModelInfo(cov);
// ERR("q not NULL in check_local -- ask author");
}
q = cov->q = (double*) malloc(maxq * sizeof(double));
for (i = 0; i < maxq; i++) q2[i] = RF_NAN; // q2 will be copied to cov->q;
// {int i; printf("%d %f\n", cov->qlen); for(i=0;i<9;i++)printf("%f ", cov->q[i]); printf("\n"); assert(false);}
C = CovList + truenext->nr; // next is GATTER only
if (p[pLOC_DIAM] == NULL) {
ERR("Diameter must always be given");
// cov->p[pLOC_DIAM] = (double*) malloc(sizeof(double));
// cov->ncol[0] = cov->nrow[0] = 1;
// cov->p[pLOC_DIAM][0] = d = GetScaledDiameter(cov);
} else {
if (cov->ncol[pLOC_DIAM] != 1 || cov->nrow[pLOC_DIAM] != 1)
ERR("diameter must be a scale");
d = cov->p[pLOC_DIAM][0];
}
if (p[pLOC_A] == NULL) {
if (C->implemented[type]) {
int i;
assert(init != NULL);
// PrintModelInfo(truenext);
init(truenext, &li);
if (li.instances == 0) {
ERR("parameter values do not allow for finding second parameter");
}
cov->p[pLOC_A] = (double*) malloc(sizeof(double));
cov->ncol[pLOC_A] = cov->nrow[pLOC_A] = 1;
q[LOCAL_R] = R_PosInf;
// printf("%f\n", d); assert(false);
msg = MSGLOCAL_FAILED;
for (i=0; i<li.instances; i++) {
if (li.msg[i] <= msg) {
err = set_local(truenext, li.value[i], d, q2);
if (err == NOERROR && (li.msg[i] < msg || q2[LOCAL_R] < q[LOCAL_R])){
memcpy(q, q2, sizeof(double) * maxq);
msg = li.msg[i];
cov->p[pLOC_A][0] = li.value[i]; // co: a; ie: rawr
}
}
}
cov->q[LOCAL_MSG] = msg;
if (msg == MSGLOCAL_FAILED) err = ERRORFAILED;
} else {
// PRINTF("%s\n", C->name);
ERR("2nd parameter is neither given nor can be found automatically");
}
} else {
if (cov->ncol[pLOC_A] != 1 || cov->nrow[pLOC_A] != 1)
ERR("`a' must be a scale");
err = set_local(truenext, cov->p[pLOC_A][0], d, q2);
memcpy(q, q2, sizeof(double) * maxq);
}
// printf("pdq %f %f %f\n", cov->p[pLOC_A][0], d, q2);
if (err != NOERROR) { XERR(err)
// char Msg[255];
// errorMessage(err);
// error(Msg);
}
return NOERROR;
}
int set_cutoff_q(cov_model *next, double a, double d, double *q) {
double phi0, phi1, a2; //, one=1.0;
a2 = a * a;
cov_fct *C = CovList + next->nr;
// PrintModelInfo(next);
C->cov(&d, next, &phi0);
C->D(&d, next, &phi1);
phi1 *= d;
if (phi0 <= 0.0) return MSGLOCAL_SIGNPHI;
if (phi1 >= 0.0) return MSGLOCAL_SIGNPHIFST;
// printf("check_co phi0=%f phi1=%f a=%f d=%f a2=%f\n", phi0, phi1, a, d, a2);
// assert(false);
// printf("p %f\n", phi1);
// printf("a2 %f\n", a2);
// printf("p0 %f\n", phi0);
// printf("a %f\n", a);
// printf("d %f\n", d);
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[LOCAL_R] = d * q[CUTOFF_THEOR];
q[CUTOFF_ASQRTR] = pow(q[LOCAL_R], a);
// printf("phi0=%f %f %f %f\n", phi0, phi1, a, d);
// printf("%f\n", d); assert(false);
// printf("%f %f %f %f\n", q[CUTOFF_B], q[CUTOFF_THEOR], q[LOCAL_R], q[CUTOFF_ASQRTR]); //assert(false);
return NOERROR;
}
int check_co(cov_model *cov) {
cov_model *next = cov->sub[0]->sub[0];// sub[0] is GATTER ONLY
cov_fct *C = CovList + next->nr;
return check_local(cov, CircEmbedCutoff, CUTOFF_MAX, C->coinit,
set_cutoff_q);
}
void Stein(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
double y=*x, z, *q=cov->q,
diameter = cov->p[pLOC_DIAM][0];
if (y <= diameter) {
C->cov(x, next, v);
*v += q[INTRINSIC_A0] + q[INTRINSIC_A2] * y * y;
// printf("%f %f %f %f %d %d \n ",
// y, v[0], q[INTRINSIC_A0], q[INTRINSIC_A2], INTRINSIC_A0, INTRINSIC_A2);
// PrintModelInfo(cov);
// assert(v[0] < RF_INF);
} else {
z = q[LOCAL_R] - y;
*v = (z <= 0.0) ? 0.0 : q[INTRINSIC_B] * z * z * z / y;
}
}
bool alternativeparam_Stein(cov_model *cov) {
cov->p[pLOC_A][0] *= 2.0;
return true;
}
void range_Stein(cov_model *cov, range_arraytype* ra)
{
range_type *range = ra->ranges;
range->min[0] = 0.0; // CUTOFF_DIAM
range->max[0] = RF_INF;
range->pmin[0] = 0.01;
range->pmax[0] = 100;
range->openmin[0] = true;
range->openmax[0] = true;
range->min[1] = 1.0; // stein_r
range->max[1] = RF_INF;
range->pmin[1] = 1.0;
range->pmax[1] = 20.0;
range->openmin[1] = false;
range->openmax[1] = true;
}
int set_stein_q(cov_model *next, double r, double d, double *q) {
double C0, phi0, phi1, phi2,
zero = 0.0,
rP1 = r + 1.0,
rM1 = r - 1.0,
dsq = d * d;
cov_fct *C = CovList + next->nr;
// printf("%ld %s\n", &zero, CovList[next->nr].name);
// PrintModelInfo(next);
C->cov(&zero, next, &C0);
C->cov(&d, next, &phi0);
C->D(&d, next, &phi1); // derivative of
phi1 *= d; // scaled fctn at 1
C->D2(&d, next, &phi2); // 2nd derivative of
phi2 *= dsq; // scaled fctn at 1
q[LOCAL_R] = r * d;
q[INTRINSIC_A2] = (phi2 - phi1) / (3.0 * r * rP1) ;
q[INTRINSIC_B] = (r == 1.0) ? 0.0 : q[INTRINSIC_A2] / (rM1 * dsq);
q[INTRINSIC_A2] = (q[INTRINSIC_A2] - phi1 / 3.0 - phi2 / 6.0) / dsq;
q[INTRINSIC_A0] = 0.5 * rM1 / rP1 * phi2 + phi1 / rP1 - phi0;
// printf("check: %f %e r=%f intr:%f %f %f\n",
// phi2, phi1, r, q[INTRINSIC_B],
// q[INTRINSIC_A0], q[INTRINSIC_A2]);
// assert(false);
if ((q[INTRINSIC_B] < 0.0) || (q[INTRINSIC_A2] < 0.0) ||
(q[INTRINSIC_A0] + C0 < 0.0))
return MSGLOCAL_INITINTRINSIC;
return NOERROR;
}
int check_Stein(cov_model *cov)
{
cov_model *next = cov->sub[0]->sub[0]; // dito
cov_fct *C = CovList + next->nr;
// printf("******************* check stein\n");
// PrintModelInfo(cov);
// assert(false);
return check_local(cov, CircEmbedIntrinsic, INTRINSIC_MAX, C->ieinit,
set_stein_q);
}
void MaStein(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[0];
double nuG, gammas, **p = cov->p, v1, v2,
nu = p[0][0], delta = p[1][0];
cov_fct *C = CovList + next->nr;
C->cov(ZERO, next, &v1);
C->cov(x + 1, next, &v2);
nuG = nu + v1 + v2;
if (nuG >= 80.0) {
error("Whittle Matern function cannot be evaluated with parameter value b+g(t) greater than 80.");
}
gammas = lgammafn(nu + delta) - lgammafn(nu) -lgammafn(nuG + delta);
double s = *x;
*v = (s==0.0) ? exp(lgammafn(nuG) + gammas)
: 2.0 * exp(nuG * log(0.5 * s) + gammas + log(bessel_k(s, nuG, 2.0)) - s);
}
int check_MaStein(cov_model *cov) {
cov_model *next = cov->sub[0];
int err;
if ((err = checkkappas(cov)) != NOERROR) return err;
cov->normalmix = cov->finiterange = false;
cov->manipulating_x=true;
if ((err = check2X(next, cov->tsdim, cov->xdim,
STATIONARY, ANISOTROPIC,
UNIVARIATE)) != NOERROR) return err;
if (cov->ncol[0] != 1 || cov->nrow[0] != 1) ERR("nu not scalar");
if (cov->ncol[1] != 1 || cov->nrow[1] != 1) ERR("d not scalar");
// if (cov->ncol[2] != 1 || cov->nrow[2] != cov->xdim - 1) {
// sprintf(ERRORSTRING, "%s",
// "dimension of vector V must be 1 less than the dim of the coordinates");
// return ERRORM;
// }
return NOERROR;
// no setbackward !
}
void range_MaStein(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = 0.0; //
range->max[0] = RF_INF;
range->pmin[0] = 1e-2;
range->pmax[0] = 10.0;
range->openmin[0] = true;
range->openmax[0] = true;
range->min[1] = 0.5 * (double) (cov->tsdim - 1); // d
range->max[1] = RF_INF;
range->pmin[1] = range->min[1];
range->pmax[1] = 10;
range->openmin[1] = false;
range->openmax[1] = true;
}
/* bivariate shift model */
void kappashift(int i, cov_model *cov, int *nr, int *nc){
*nc = 1;
*nr = (i<=0) ? cov->tsdim : -1;
}
void shift(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
double y[ShiftMaxDim],
*h = cov->p[0];
int i,
dim = cov->tsdim;
C->cov(x, next, v + 0);
for (i=0; i<dim; i++) y[i] = x[i] + h[i];
C->cov(y, next, v + 1);
for (i=0; i<dim; i++) y[i] = x[i] - h[i];
C->cov(y, next, v + 2);
v[3] = v[0];
}
int checkshift(cov_model *cov) {
cov_model *next = cov->sub[0];
// cov_fct *C = CovList + next->nr;
int err,
dim = cov->tsdim;
if (cov->xdim > ShiftMaxDim) {
sprintf(ERRORSTRING,
"For technical reasons max. dimension for ave is %d. Got %d.",
StpMaxDim, cov->xdim);
return ERRORMSG;
}
if ((err = checkkappas(cov)) != NOERROR) return err;
cov->manipulating_x=true;
cov->normalmix = cov->finiterange = false;
if ((err = check2X(next, dim, dim, STATIONARY,
(dim > 1) ? ANISOTROPIC : ISOTROPIC,
UNIVARIATE)) != NOERROR)
return err;
setbackward(cov, next);
return NOERROR;
}
void rangeshift(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = -RF_INF;
range->max[0] = RF_INF;
range->pmin[0] = -1000;
range->pmax[0] = 1000;
range->openmin[0] = true;
range->openmax[0] = true;
}
/* Vector - Delta Delta^T */
void vector(double *x, cov_model *cov, double *v) {
/*
r = \| x\|
d/d x_i f(r) = f'(r) x_i / r
d^2/d x_i^2 f(r) = f''(r) x_i^2/r^2 + f'(r) (r^2 - x_i^2) / r^3
d^2/d x_i x_j f(r) = f''(r) x_i x_j/r^2 + f'(r) (- x_i x_j) / r^3
r = 0:
d/d x_i f(r) |_0 = 0
d^2/d x_i^2 f(r) |_0 = f''(0)
d^2/d x_i x_j f(r) = 0
Vector (n = dimension of spatial field)
r=0
L = n * f''(0)
operator: -0.5 * (a + 1) laplace + a * hessian, a in [-1,1]; a=1
*/
cov_model *next = cov->sub[0];
cov_fct *N = CovList + next->nr;
double norm[2], normSq0, normL2, normT2, D, D2, diag,
a = cov->p[0][0],
b = - 0.5 * (1.0 + a);
int i,
td = cov->tsdim,
dim = ((int *)cov->p[1])[0],
dimP1 = dim + 1,
dimsq = dim * dim;
normL2 = normT2 = 0.0;
for (i=0; i<dim; i++) normL2 += x[i] * x[i];
for (; i<td; i++) normT2 += x[i] * x[i];
if (next->isoIn == ISOTROPIC) {
normSq0 = normL2 + normT2;
} else {
normSq0 = normL2;
norm[1] = sqrt(normT2);
}
norm[0] = sqrt(normSq0);
N->D(norm, next, &D);
N->D2(norm, next, &D2);
if (normSq0 == 0.0) {
diag = (b * dim + a) * D2;
for (i=0; i<dimsq; i++) {
// printf("%f D2 %f\n", diag, D2);
v[i] = (i % dimP1 == 0) ? diag : 0.0;
}
} else {
double
D2nsq = D2 / normSq0,
D1n = D / norm[0],
D1n3 = D / (norm[0] * normSq0),
delta = a * (D2nsq - D1n3),
Vector = b * ((D2nsq - D1n3) * normL2 + dim * D1n);
int k,l;
diag = Vector + a * D1n;
for (i=k=0; k<dim; k++) {
for (l=0; l<dim; l++, i++) {
v[i] = (i % dimP1 == 0) ? diag : 0.0;
v[i] += delta * x[k] * x[l];
}
}
}
}
/* Vector - Delta Delta^T */
void vectorAniso(double *x, cov_model *cov, double *v) {
/*
operator: -0.5 * (a + 1) laplace + a * hessian, a in [-1,1]; a=1
*/
cov_model *next = cov->sub[0];
cov_fct *N = CovList + next->nr;
double *D, laplace,
a = cov->p[0][0],
b = - 0.5 * (1.0 + a);
int i, j, k, endfor,
dim = ((int *)cov->p[1])[0],
dimP1 = dim + 1,
dimsq = dim * dim,
xdim = cov->tsdim,
xdimsq = xdim * xdim,
xdimP1 = xdim + 1,
dimxdim = dim * xdim;
D = (double*) malloc(sizeof(double) * xdimsq);
//printf("%d %d\n", dim, xdim); assert(false);
// should pass back the hessian
N->hess(x, next, D);
laplace = 0.0;
// PrintModelInfo(next);
// printf("return matrix %d\n", xdimsq);
// for (i=0; i<xdimsq; i+=xdim) {
// for (j=i; j <i+xdim; j++) {
// printf("%f;%d ", D[j], j);
// }
// printf("\n");
// }
for (i=0; i<dimxdim; i+=xdimP1) {
laplace += D[i];
// printf("lplace %d %f\n", i, D[i]);
}
laplace *= b;
k = 0;
for (i=0; i<dimxdim; i+=xdim) {
endfor = i + dim;
for (j=i; j<endfor; j++) {
// printf(">> %d %d %f %f %f\n", k, j, D[j], a, laplace );
v[k++] = D[j] * a;
}
}
for (i=0; i<dimsq; i+=dimP1) {
v[i] += laplace;
}
free(D);
}
int checkvector(cov_model *cov) {
cov_model *next = cov->sub[0],
*genuinenext = next->sub[0];
// cov_fct *N = CovList + next->sub[0]->nr;
int err,
dim = cov->tsdim;
isotropy_type isotropy;
cov->nr = VECTOR; // wegen nr++ unten;
cov->normalmix = cov->finiterange = false;
cov->manipulating_x=true;
if ((err = check2X(next, dim, 1, STATIONARY, ISOTROPIC,
UNIVARIATE)) != NOERROR) {
// printf("err=%d\n", err);
if ((err = check2X(next, dim, dim, STATIONARY, ZEROSPACEISO,
UNIVARIATE)) != NOERROR) {
// printf("err2=%d\n", err);
if ((err = check2X(next, dim, dim, STATIONARY, ANISOTROPIC,
UNIVARIATE)) != NOERROR) {
// printf("err3=%d\n", err);
return err;
}
}
}
assert(cov->xdim == cov->tsdim);
if (genuinenext->derivatives < 2 && !genuinenext->hess) {
// PrintModelInfo(cov);
// printf("%s (%d): derivatives %d \n", CovList[genuinenext->nr].name,
// genuinenext->nr, genuinenext->derivatives);
// next is GATTER
ERR("2nd derivative of submodel not defined (for the given paramters)");
}
isotropy = next->isoIn;
if (isotropy != ISOTROPIC && isotropy != SPACEISOTROPIC) {
if (!genuinenext->hess) ERR("hess matrix not defined");
cov->nr++;
}
kdefault(cov, 0, 0.5); // a
kdefault(cov, 1,
(isotropy==SPACEISOTROPIC || isotropy==ZEROSPACEISO)
? dim - 1 : dim); // Dspace
if ( (isotropy==SPACEISOTROPIC || isotropy==ZEROSPACEISO)
&& ((int*) (cov->p[1]))[0] != dim - 1) {
// printf(".... %d %d\n", ((int*) (cov->p[1]))[0], dim);
ERR("for spatiotemporal models 'vector' must be applied to spatial part");
}
// PrintModelInfo(cov);
setbackward(cov, next);
cov->vdim = ((int*) (cov->p[1]))[0];
next->tsdim = DEL_COV;
return NOERROR;
}
void rangevector(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = -1.0;
range->max[0] = 1.0;
range->pmin[0] = -1.0;
range->pmax[0] = 1.0;
range->openmin[0] = false;
range->openmax[0] = false;
range->min[1] = 1.0;
range->max[1] = cov->tsdim;
range->pmin[1] = 1.0;
range->pmax[1] = cov->tsdim;
range->openmin[1] = false;
range->openmax[1] = false;
}
int zzz=0;
/* Rot - Delta Delta^T */
void curl(double *x, cov_model *cov, double *v) {
/*
r = \| x\|
d/d x_i f(r) = f'(r) x_i / r
d^2/d x_i^2 f(r) = f''(r) x_i^2/r^2 + f'(r) (r^2 - x_i^2) / r^3
d^2/d x_i x_j f(r) = f''(r) x_i x_j/r^2 + f'(r) (- x_i x_j) / r^3
r = 0:
d/d x_i f(r) |_0 = 0
d^2/d x_i^2 f(r) |_0 = f''(0)
d^2/d x_i x_j f(r) = 0
Rot (n = dimension of spatial field)
r=0
L = n * f''(0)
operator: -0.5 * (a + 1) rot + a * hessian, a in [-1,1]; a=1/2
*/
cov_model *next = cov->sub[0];
cov_fct *N = CovList + next->nr;
double norm[2], normSq0, normL2, normT2, D, D2, D3, diag,
a = -1.0, // curl free
b = - 0.5 * (1.0 + a);
int i,
td = cov->tsdim,
dim = td, // ((int *)cov->p[1])[0],
// currently no space-time allowed -> BA thesis
dimP1 = dim + 1,
dimP2 = dim + 2,
dimP3 = dim + 3,
dimP2sqM1 = dimP2 * dimP2 -1; // for three dimensions much
// more complicated
normL2 = normT2 = 0.0;
for (i=0; i<dim; i++) normL2 += x[i] * x[i];
for (; i<td; i++) normT2 += x[i] * x[i];
if (next->isoIn == ISOTROPIC) {
normSq0 = normL2 + normT2;
} else {
normSq0 = normL2;
norm[1] = sqrt(normT2);
}
norm[0] = sqrt(normSq0);
N->D(norm, next, &D);
N->D2(norm, next, &D2);
N->D3(norm, next, &D3);
if (normSq0 == 0.0) {
for (i=0; i<=dimP2sqM1; i++) v[i] = 0.0;
N->cov(norm, next, v); // v[0]
diag = (b * dim + a) * D2;
for (i=dimP3; i<dimP2sqM1; i+=dimP3) v[i] = diag; // diag 2,3
N->D2(norm, next, v + dimP1);
v[dimP1] *= 2.0;
v[dimP2 * dimP1] = v[dimP1];
N->D4(norm, next, v + dimP2sqM1);
v[dimP2sqM1] *= (8.0 / 3.0);
} else {
double
z[2],
D2nsq = D2 / normSq0,
D1n = D / norm[0],
D3n = D3 / norm[0],
D1n3 = D / (norm[0] * normSq0),
delta = a * (D2nsq - D1n3),
div = b * ((D2nsq - D1n3) * normL2 + dim * D1n);
int k,l;
N->cov(norm, next, v); // v[0]
z[0] = x[0];
z[1] = x[1];
for (i=1; i<=dim; i++) { // kante links und oben
v[i] = -(v[i * dimP2] = z[i-1] * D1n);
}
diag = div + a * D1n;
for (i= dimP3, k=0; k<dim; k++, i+=2) {
for (l=0; l<dim; l++, i++) {
v[i] = (i % dimP3 == 0) ? diag : 0.0;
v[i] += delta * x[k] * x[l];
}
}
v[dimP1 * dimP2] = v[dimP1] = - v[dimP3] - v[2 * dimP3]; // ecken
for (i=1; i<=dim; i++) { // kanten unten und hinten
v[dimP2 * dimP1 + i] =
- (v[dimP2 * (i+1) - 1] = z[i-1] * (D3n + D2nsq - D1n3));
}
N->D4(norm, next, v + dimP2sqM1); // rechts unten
v[dimP2sqM1] += 2.0 * D3n - D2nsq + D1n3;
}
// for (i=0; i<=15; i++) {
// if (i!=15) v[i] = 0.0;
// }
}
/* Div - Delta Delta^T */
void div(double *x, cov_model *cov, double *v) { // div -free !!
/*
r = \| x\|
d/d x_i f(r) = f'(r) x_i / r
d^2/d x_i^2 f(r) = f''(r) x_i^2/r^2 + f'(r) (r^2 - x_i^2) / r^3
d^2/d x_i x_j f(r) = f''(r) x_i x_j/r^2 + f'(r) (- x_i x_j) / r^3
r = 0:
d/d x_i f(r) |_0 = 0
d^2/d x_i^2 f(r) |_0 = f''(0)
d^2/d x_i x_j f(r) = 0
Div (n = dimension of spatial field)
r=0
L = n * f''(0)
operator: -0.5 * (a + 1) div + a * hessian, a in [-1,1]; a=1
*/
cov_model *next = cov->sub[0];
cov_fct *N = CovList + next->nr;
double norm[2], normSq0, normL2, normT2, D, D2, D3, diag,
a = 1.0, // divergenz free
b = - 0.5 * (1.0 + a);
int i,
td = cov->tsdim,
dim = td, // ((int *)cov->p[1])[0],
// currently no space-time allowed -> BA thesis
dimP1 = dim + 1,
dimP2 = dim + 2,
dimP3 = dim + 3,
dimP2sqM1 = dimP2 * dimP2 -1; // for three dimensions much
// more complicated
normL2 = normT2 = 0.0;
for (i=0; i<dim; i++) normL2 += x[i] * x[i];
for (; i<td; i++) normT2 += x[i] * x[i];
if (next->isoIn == ISOTROPIC) {
normSq0 = normL2 + normT2;
} else {
normSq0 = normL2;
norm[1] = sqrt(normT2);
}
norm[0] = sqrt(normSq0);
N->D(norm, next, &D);
N->D2(norm, next, &D2);
N->D3(norm, next, &D3);
if (normSq0 == 0.0) {
for (i=0; i<=dimP2sqM1; i++) v[i] = 0.0;
N->cov(norm, next, v); // v[0]
diag = (b * dim + a) * D2;
for (i=dimP3; i<dimP2sqM1; i+=dimP3) v[i] = diag; // diag 2,3
N->D2(norm, next, v + dimP1);
v[dimP1] *= 2.0;
v[dimP2 * dimP1] = v[dimP1];
N->D4(norm, next, v + dimP2sqM1);
v[dimP2sqM1] *= (8.0 / 3.0);
} else {
double
z[2],
D2nsq = D2 / normSq0,
D1n = D / norm[0],
D3n = D3 / norm[0],
D1n3 = D / (norm[0] * normSq0),
delta = a * (D2nsq - D1n3),
div = b * ((D2nsq - D1n3) * normL2 + dim * D1n);
int k,l;
N->cov(norm, next, v); // v[0]
z[0] = -x[1];
z[1] = x[0];
for (i=1; i<=dim; i++) { // kante links und oben
v[i] = -(v[i * dimP2] = z[i-1] * D1n);
}
diag = div + a * D1n;
for (i= dimP3, k=0; k<dim; k++, i+=2) {
for (l=0; l<dim; l++, i++) {
v[i] = (i % dimP3 == 0) ? diag : 0.0;
v[i] += delta * x[k] * x[l];
}
}
v[dimP1 * dimP2] = v[dimP1] = - v[dimP3] - v[2 * dimP3]; // ecken
for (i=1; i<=dim; i++) { // kanten unten und hinten
v[dimP2 * dimP1 + i] =
-(v[dimP2 * (i+1) - 1] = z[i-1] * (D3n + D2nsq - D1n3));
}
N->D4(norm, next, v + dimP2sqM1); // rechts unten
v[dimP2sqM1] += 2.0 * D3n - D2nsq + D1n3;
}
// for (i=0; i<=15; i++) {
// if (i!=15) v[i] = 0.0;
// }
}
int checkdivcurl(cov_model *cov) {
cov_model *next = cov->sub[0],
*nexttrue=next->sub[0]; // next is gatter and derivatives have not
// been programmed yet
// cov_fct *N = CovList + next->sub[0]->nr;
int err,
dim = cov->tsdim,
spacedim;
isotropy_type isotropy;
cov->normalmix = cov->finiterange = false;
cov->manipulating_x=true;
if ((err = check2X(next, dim, 1, STATIONARY, ISOTROPIC,
UNIVARIATE)) != NOERROR) {
if ((err = check2X(next, dim, 1, STATIONARY, SPACEISOTROPIC,
UNIVARIATE)) != NOERROR)
return err;
}
if (nexttrue->derivatives < 4) ERR("4th derivative of submodel not defined");
if (cov->tsdim != 2) ERR("currently coded only for dim=2");
isotropy = next->isoIn;
spacedim = dim - (isotropy == SPACEISOTROPIC);
if (isotropy != ISOTROPIC && isotropy != SPACEISOTROPIC)
ERR("submodel must be spaceisotropic");
if (spacedim != 2)
ERR("model currently only defined for the plane");
setbackward(cov, next);
cov->vdim = spacedim + 2; // only correct for spacedim == 2!!!
assert(spacedim == 2);
next->tsdim = DEL_COV;
return NOERROR;
}
void rangedivcurl(cov_model *cov, range_arraytype* ra){ }
void lp(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
double z,
p = cov->p[0][0];
int d,
dim = cov->tsdim;
for (z=0.0, d=0; d<dim; d++) {
z += pow(fabs(x[d]), p);
}
z = pow(z, 1 / p);
CovList[next->nr].cov(&z, next, v);
}
int checklp(cov_model *cov) {
bool skipchecks = GLOBAL.general.skipchecks;
cov_fct *C = CovList + cov->nr;
cov_model *next = cov->sub[0];
double p;
int err;
kdefault(cov, 0, 1.0);
p =cov->p[0][0];
if ((err = checkkappas(cov)) != NOERROR) return err;
cov->manipulating_x=true;
strcpy(ERROR_LOC, C->name);
if ((err = check2X(next, cov->tsdim + 1, 1, STATIONARY, ISOTROPIC,
UNIVARIATE))
!= NOERROR) return err;
next->tsdim = DEL_COV;
if (p==1) {
cov->maxdim = next->maxdim - 1; // true for p==1
if (cov->maxdim == 0) cov->maxdim = 1;
} else if (p==2) {
cov->maxdim = next->maxdim;
} else if (!skipchecks) ERR("p must be 1 or 2.");
cov->normalmix = next->normalmix;
cov->semiseparatelast = false;
updatepref(cov, next);
return NOERROR;
}
void rangelp(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = 0.0;
range->max[0] = 2.0;
range->pmin[0] = 0.01;
range->pmax[0] = 2.0;
range->openmin[0] = true;
range->openmax[0] = false;
}
void ma1(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
double z,
alpha = cov->p[0][0],
theta = cov->p[1][0];
CovList[next->nr].cov(x, next, &z);
*v = pow(theta / (1 - (1-theta) * z), alpha);
}
int checkma1(cov_model *cov) {
// bool skipchecks = GLOBAL.general.skipchecks;
cov_fct *C = CovList + cov->nr;
cov_model *next = cov->sub[0];
// double p;
int err;
kdefault(cov, 0, 1.0);
kdefault(cov, 1, 0.5);
if ((err = checkkappas(cov)) != NOERROR) return err;
strcpy(ERROR_LOC, C->name);
if ((err = check2X(next, cov->tsdim, cov->xdim, cov->statIn,
cov->isoIn, UNIVARIATE))
!= NOERROR) return err;
cov->semiseparatelast = false;
// updatepref(cov, next);
return NOERROR;
}
void rangema1(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = 0.0;
range->max[0] = RF_INF;
range->pmin[0] = 0.01;
range->pmax[0] = 10;
range->openmin[0] = true;
range->openmax[0] = true;
range->min[1] = 0.0;
range->max[1] = 1.0;
range->pmin[1] = 0.0;
range->pmax[1] = 1.0;
range->openmin[1] = true;
range->openmax[1] = true;
}
void ma2(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
double z,z0;
CovList[next->nr].cov(ZERO, next, &z0);
CovList[next->nr].cov(x, next, &z);
z = z0 - z;
*v = (z == 0) ? 1.0 : (1.0 - exp(-z)) / z;
}
int checkma2(cov_model *cov) {
// bool skipchecks = GLOBAL.general.skipchecks;
cov_fct *C = CovList + cov->nr;
cov_model *next = cov->sub[0];
// double p;
int err;
if ((err = checkkappas(cov)) != NOERROR) return err;
cov->manipulating_x = true;
strcpy(ERROR_LOC, C->name);
if ((err = check2X(next, cov->tsdim, cov->xdim, VARIOGRAM,
cov->isoIn, UNIVARIATE))
!= NOERROR) return err;
cov->semiseparatelast = false;
// updatepref(cov, next);
return NOERROR;
}
void rangema2(cov_model *cov, range_arraytype* ra){
}
void kappaM(int i, cov_model *cov, int *nr, int *nc){
*nc = 0;
*nr = i < 1 ? 0 : -1;
}
void M(cov_model *cov, double *Z, double *V) {
// M^t C M
cov_model *next = cov->sub[0];
int
*ncol = cov->ncol,
*nrow = cov->nrow;
double *MtZ,
alpha = 1.0,
beta = 0.0,
*M = cov->p[0];
// SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
// C := alpha*op( A )*op( B ) + beta*C,
// M : rows of op(A)
// N : cols of op(B)
// K : cols of op(A)= rows of op(B)
// LD-X : first dimension of X
if (next->vdim == 1) {
MtZ = (double*) malloc(sizeof(double) * *ncol * *ncol);
F77_CALL(dgemm)("T", "N", ncol, ncol, nrow, Z, M, nrow, M, nrow,
&beta, V, ncol);
} else {
MtZ = (double*) malloc(sizeof(double) * *nrow * *ncol);
F77_CALL(dgemm)("T", "N", ncol, nrow, nrow, &alpha, M, nrow, Z, nrow,
&beta, MtZ, ncol);
F77_CALL(dgemm)("N", "N", ncol, ncol, nrow, &alpha, MtZ, ncol, M, nrow,
&beta, V, ncol);
}
free(MtZ);
}
void Mstat(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
int ncol = cov->ncol[0];
double *z = (double*) malloc(sizeof(double) * ncol * ncol);
CovList[next->nr].cov(x, next, z);
M(cov, z, v);
free(z);
}
void Mnonstat(double *x, double *y, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
int ncol = cov->ncol[0];
double *z = (double*) malloc(sizeof(double) * ncol * ncol);
CovList[next->nr].nonstat_cov(x, y, next, z);
M(cov, z, v);
free(z);
}
int checkM(cov_model *cov) {
cov_model *next = cov->sub[0];
int err;
if ((err = check2X(next, cov->tsdim, cov->xdim, cov->statIn,
cov->isoIn, cov->nrow[0])) != NOERROR) return err;
setbackward(cov, next);
cov->vdim = cov->ncol[0];
err = checkkappas(cov);
return err;
}
void rangeM(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = RF_NEGINF;
range->max[0] = RF_INF;
range->pmin[0] = -1e10;
range->pmax[0] = 1e10;
range->openmin[0] = true;
range->openmax[0] = true;
range->finiterange = true;
}
void IdStat(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[0];
CovList[next->nr].cov(x, next, v);
}
void IdNonStat(double *x, double *y, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
CovList[next->nr].nonstat_cov(x, y, next, v);
}
int checkId(cov_model *cov) {
cov_model *next = cov->sub[0];
int err;
cov->vdim = (cov->nrow[0] > 0 && cov->nrow[0] > 0) ?
((int*) cov->p[0])[0] : SUBMODELM;
if ((err = check2X(next, cov->tsdim, cov->xdim, cov->statIn,
cov->isoIn, cov->vdim
)) != NOERROR) return err;
if (cov->vdim == SUBMODELM) cov->vdim = next->vdim;
setbackward(cov, next);
return NOERROR;
}
void DId(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
CovList[next->nr].D(x, next, v);
}
void DDId(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
CovList[next->nr].D2(x, next, v);
}
void TBM2Id(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
CovList[next->nr].tbm2(x, next, v);
}
int initspectralId(cov_model *cov) {
cov_model *next = cov->sub[0];
return CovList[next->nr].initspectral(next);
}
void spectralId(cov_model *cov, spectral_storage *s, double *e) {
cov_model *next = cov->sub[0];
return CovList[next->nr].spectral(next, s, e);
}
void coinitId(cov_model *cov, localinfotype *li) {
cov_model *next = cov->sub[0];
CovList[next->nr].coinit(next, li);
}
void ieinitId(cov_model *cov, localinfotype *li) {
cov_model *next = cov->sub[0];
CovList[next->nr].ieinit(next, li);
}
void rangeId(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = 1.0;
range->max[0] = RF_INF;
range->pmin[0] = 1.0;
range->pmax[0] = 10;
range->openmin[0] = false;
range->openmax[0] = true;
range->finiterange = true;
range->maxdim = INFDIM;
}
void Exp(double *x, cov_model *cov, double *v){
double v0, v1;
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
C->cov(ZERO, next, &v0);
C->cov(x, next, &v1);
*v = exp(v1-v0); // == - (v0 - v1)
}
void DExp(double *x, cov_model *cov, double *v){
double D;
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
C->D(x, next, &D);
Exp(x, cov, v);
*v *= -D;
}
void DDExp(double *x, cov_model *cov, double *v){
double D, D2;
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
C->D(x, next, &D);
C->D2(x, next, &D2);
Exp(x, cov, v);
*v *= D * D + D2; // Achtung + D2 nicht - D2, da D2 Ableitung einer Cov.
}
int checkExp(cov_model *cov) {
cov_model *next = cov->sub[0];
int err;
if (cov->statIn != STATIONARY && cov->statIn != VARIOGRAM)
return ERRORSTATVARIO;
if ((err = check2X(next, cov->tsdim, cov->xdim, cov->statIn,
cov->isoIn, UNIVARIATE)) != NOERROR) return err;
cov->manipulating_x=true;
setbackward(cov, next);
cov->normalmix = false;
return NOERROR;
}
void rangeExp(cov_model *cov, range_arraytype* ra) { }
void Pow(double *x, cov_model *cov, double *v){
double v0, v1,
alpha = cov->p[0][0];
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
C->cov(ZERO, next, &v0);
C->cov(x, next, &v1);
*v = v0 - pow(v0 - v1, alpha);
// printf("Pow %f %f %f v=%f\n", v0, v1, alpha, *v);
}
void DPow(double *x, cov_model *cov, double *v){
double v0, v1, gamma,
alpha = cov->p[0][0];
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
C->D(x, next, v);
if (alpha == 1.0) return;
C->cov(ZERO, next, &v0);
C->cov(x, next, &v1);
gamma = v0 - v1;
*v *= - alpha * pow(gamma, alpha -1.0);
// Achtung "-" , da covarianzfunktion angenommen
}
void DDPow(double *x, cov_model *cov, double *v){
double D, v0, v1, gamma,
alpha = cov->p[0][0];
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
C->D2(x, next, v);
if (alpha == 1.0) return;
C->D(x, next, &D);
C->cov(ZERO, next, &v0);
C->cov(x, next, &v1);
gamma = v0 - v1;
*v *= - alpha * pow(gamma, alpha - 2.0) * ((alpha - 1.0) * D + gamma * (*v));
// Achtung "-" , da covarianzfunktion angenommen
}
void invPow(double *x, cov_model *cov, double *v) {
// invPow only defined for submodel having variance 1 !!!!
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
C->cov(x, next, v);
double y = 1.0 - *v;
if (y < 0.0 || y > 1.0) {
if (y > -1e-14 && y < 0.0) y=0.0;
else if (y < 1.0 + 1e-14) y=1.0;
else {
//PRINTF("covariance value %e (1 - %e = %e) at %e\n", *v, *v, 1.0-*v, *x);
ERR("invPow valid only for non-negative covariance models with variance 1");
}
}
*v = 1.0 - pow(y, 1.0 / cov->p[0][0]);
// printf("%f %f %f\n", y, *v, cov->p[0][0]);
// printf("Inv %f %f %f v=%f\n", *x, cov->p[0][0],y, *v);
}
int checkPow(cov_model *cov) {
cov_model *next = cov->sub[0];
// cov_fct *C = CovList + next->nr;
int err;
//printf("OK %d %d\n", cov->statIn, cov->isoIn);
cov->manipulating_x=true;
if (cov->statIn != STATIONARY && cov->statIn != VARIOGRAM)
return ERRORSTATVARIO;
if ((err = check2X(next, cov->tsdim, cov->xdim, cov->statIn,
cov->isoIn, UNIVARIATE)) != NOERROR){
// printf("OK %d\n",err);
return err;
}
// printf("xOK %d\n",err);
setbackward(cov, next);
cov->normalmix = false;
if (cov->statIn == STATIONARY && next->normalmix) {
cov->normalmix = true;
}
return NOERROR;
}
void rangePow(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->finiterange = true;
range->maxdim = INFDIM;
range->min[0] = 0.0;
range->max[0] = 1.0;
range->pmin[0] = 0.01;
range->pmax[0] = 1.0;
range->openmin[0] = true;
range->openmax[0] = false;
}
void tbm3(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
double v1;
cov_fct *C = CovList + next->nr;
C->cov(x, next, v);
C->D(x, next, &v1);
(*v) += *x * v1 / cov->p[0][0];
}
int checktbm3(cov_model *cov) {
cov_model *next = cov->sub[0];
assert(next != NULL);
int err, newdim;
kdefault(cov, 0, 1);
newdim = (int) cov->p[0][0];
if (cov->tsdim > newdim) return ERRORWRONGDIM; /// ??????
if (cov->isoIn != ISOTROPIC && cov->isoIn != SPACEISOTROPIC)
return ERRORANISO;
if ((err = check2X(next, 2 + newdim, cov->xdim, STATIONARY, cov->isoIn,
UNIVARIATE))
!= NOERROR) return err;
if (cov->isoIn != SPACEISOTROPIC || next->sub[0]->isoIn == SPACEISOTROPIC)
next->tsdim = DEL_COV;
cov->derivatives = next->derivatives - 1;
if (cov->derivatives < 0) {
PrintModelInfo(cov); // ok
ERR("derivative for the submodel does not exist or is unknown");
}
cov->manipulating_x=true;
setbackward(cov, next);
cov->pref[CircEmbed] = next->pref[CircEmbed] + 1;
if (cov->pref[CircEmbed] > PREF_BEST) cov->pref[CircEmbed] = PREF_BEST;
return NOERROR;
}
void rangetbm3(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->finiterange = true;
range->maxdim = INFDIM;
range->min[0] = 1.0;
range->max[0] = RF_INF;
range->pmin[0] = 1.0;
range->pmax[0] = 1000000.0;
range->openmin[0] = false;
range->openmax[0] = true;
}
void tbm2(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
if (C->tbm2 == NULL)
ERR("tbm2 operator is unavailable for the specified submodel");
C->tbm2(x, next, v);
}
int checktbm2(cov_model *cov){
cov_model *next = cov->sub[0];
assert(next != NULL);
bool layer = cov->xdim == 2;
int err;
cov->nr = TBM2NR;
cov->manipulating_x=true;
if (cov->tsdim > 2 + layer) return ERRORWRONGDIM;
assert(cov->xdim == 1 + layer);
if ((err = check2X(next, 2 + layer, cov->xdim, STATIONARY, cov->isoIn,
UNIVARIATE)) != NOERROR) return err;
setbackward(cov, next);
if (next->pref[TBM2] <= 0) {
return ERRORCURRENTLY;
cov->nr ++; // numerical evaluation ...
}
return NOERROR;
}
void rangetbm2(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->finiterange = false;
range->maxdim = 1;
}
typedef struct TBM2_integr{
cov_model *cov;
double *x;
bool layer;
} TBM2_integr;
void TBM2NumIntegrFct(double *u, int n, void *ex) {
int i;
double z[2];
TBM2_integr *info;
info = (TBM2_integr*) ex;
cov_model *cov=info->cov;
double *x = info->x;
for (i=0; i<n; i++) {
z[0] = x[0] * sqrt(1 - u[i] * u[i]);
tbm3(z, cov, u + i);
}
}
void tbm2num(double *x, cov_model *cov, double *v){
TBM2_integr info;
#define MAXSUBDIVISIONS 100
static double a = 0,
b = 1,
eps = 1e-10; // 1e-15
static 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);
}
/* undefined model -- for technical reasons only */
void rangeundefined(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->maxdim = 0;
assert(false);
}
int checkundefined(cov_model *cov) {
assert(false);
// info->CEbad = true;
}
/* qam */
void kappaqam(int i, cov_model *cov, int *nr, int *nc){
*nc = 1;
*nr = (i==0) ? cov->nsub-1 : -1;
}
void qam(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
int i,
nsub = cov->nsub;
double sum, s, w,
*theta = cov->p[0];
sum = 0.0;
for (i=1; i<nsub; i++) {
cov_model * sub = cov->sub[i];
CovList[sub->nr].cov(x, sub, &s);
C->invSq(&s, next, &w);
sum += theta[i - 1] * w;
}
sum = sqrt(sum);
C->cov(&sum, next, v);
}
int checkqam(cov_model *cov) {
cov_model *next = cov->sub[0],
*sub;
int i, err,
nsub = cov->nsub,
nsubM1 = nsub - 1;
double v, sum;
if ((err = checkkappas(cov)) != NOERROR) return err;
sum = 0.0;
for (i=0; i<nsubM1; i++) {
sum += cov->p[0][i];
}
if (fabs(sum - 1.0) > 1e-14) ERR("theta must sum up to 1");
cov->manipulating_x=true;
cov->normalmix = true;
cov->finiterange = false;
if ((err = check2X(next, 1, 1, STATIONARY, ISOTROPIC,
UNIVARIATE)) != NOERROR)
return err;
if (!next->normalmix) ERR("phi is not completely monotone");
for (i=1; i<nsub; i++) {
sub = cov->sub[i];
if ((err = check2X(sub, cov->tsdim, cov->tsdim, STATIONARY, ANISOTROPIC,
UNIVARIATE)) != NOERROR)
return err;
CovList[sub->nr].cov(ZERO, sub, &v);
if (v != 1.0) return ERRORUNITVAR;
setbackward(cov, sub);
}
return NOERROR;
}
void rangeqam(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
range->min[0] = 0.0;
range->max[0] = 1.0;
range->pmin[0] = 0.0;
range->pmax[0] = 1.0;
range->openmin[0] = false;
range->openmax[0] = false;
}
/* mqam */
void kappamqam(int i, cov_model *cov, int *nr, int *nc) {
int nsub = cov->nsub - 1;
*nc = (i==0) ? 1 : nsub;
*nr = (i <= 1) ? nsub : -1;
// printf("kappa %d nsub=%d %d %d\n", i, nsub, *nc, *nr);
}
void mqam(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
int i, j, k, l,
nsub = cov->nsub,
nsubM1 = nsub - 1;
double s0,
theta = cov->p[0][0],
EMtheta = 1.0 - theta,
*rho = cov->p[1],
s[MAXSUB];
for (i=1; i<nsub; i++) {
cov_model * sub = cov->sub[i];
CovList[sub->nr].cov(x, sub, &s0);
C->invSq(&s0, next, s + i - 1);
}
for (j=0; j<nsubM1; j++) {
l = k = j * nsub; // nicht nsubM1
for (i=j; i<nsubM1; i++, k++, l+=nsubM1) {
s0 = theta * s[i] + EMtheta * s[j];
s0 = sqrt(s0);
C->cov(&s0, next, v + k);
v[l] = v[k] * rho[l];
if (k != l) v[k] *= rho[k];
assert (rho[l] == rho[k]) ;
// printf("%d %d %d %d %f %f %e %e\n", i, j, k, l,
// v[k], rho[l], v[l]-v[k], rho[k]-rho[l]);
}
}
}
int checkmqam(cov_model *cov) {
// cov_model *next = cov->sub[0];
int err,
nsub = cov->nsub,
nsubM1 = nsub - 1;
if ((err = checkqam(cov)) != NOERROR) return err;
cov->vdim = nsubM1;
cov->manipulating_x=true;
return NOERROR;
}
void rangemqam(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges;
double pmax = 2.0 / (double) (cov->nsub-1);
range->min[0] = 0.0;
range->max[0] = 1.0;
range->pmin[0] = 0.0;
range->pmax[0] = pmax;
range->openmin[0] = false;
range->openmax[0] = false;
range->min[1] = RF_NEGINF;
range->max[1] = RF_INF;
range->pmin[1] = -10.0;
range->pmax[1] = 10.0;
range->openmin[1] = true;
range->openmax[1] = true;
}
/////////////// NATSC
void natsc(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
double scale, y;
if (next->nr >= GATTER && next->nr <= LASTGATTER) {
/* diese Abfrage muss bleiben, da roher aufruf ohne
voriges DeleteGatter aus stein_st_check heraus
erfolgt.
Wenn Gatter gar nicht eingefuehrt wuerde,
kann natsc nicht ueber checkX2 das primitive aufrufen
*/
next = next->sub[0];
}
assert(CovList[next->nr].naturalscale != NULL);
scale = CovList[next->nr].naturalscale(next);
y = *x / scale;
// letzteres darf nur passieren wenn dim = 1!!
CovList[next->nr].cov(&y, next, v);
}
void Dnatsc(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
int i,
vdimSq = cov->vdim * cov->vdim;
double scale, y;
if (next->nr >= GATTER && next->nr <= LASTGATTER)
next = next->sub[0];
assert(CovList[next->nr].naturalscale != NULL);
scale = CovList[next->nr].naturalscale(next);
y = *x / scale;
CovList[next->nr].D(&y, next, v);
for (i=0; i<vdimSq; i++) v[i] /= scale;
}
void DDnatsc(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
int i,
vdimSq = cov->vdim * cov->vdim;
double scale, y, ScSq;
if (next->nr >= GATTER && next->nr <= LASTGATTER)
next = next->sub[0];
assert(CovList[next->nr].naturalscale != NULL);
scale = CovList[next->nr].naturalscale(next);
y = *x / scale;
ScSq = scale * scale;
CovList[next->nr].D2(&y, next, v);
for (i=0; i<vdimSq; i++) v[i] /= ScSq;
}
void invnatsc(double *x, cov_model *cov, double *v) {
cov_model *next = cov->sub[0];
int i,
vdimSq = cov->vdim * cov->vdim;
double
scale = CovList[next->nr].naturalscale(next);
CovList[next->nr].D(x, next, v);
for (i=0; i<vdimSq; i++) v[i] *= scale;
}
void coinitnatsc(cov_model *cov, localinfotype *li) {
cov_model *next = cov->sub[0]->sub[0];
if ( CovList[next->nr].coinit == NULL)
error("# cannot find coinit -- please inform author");
CovList[next->nr].coinit(next, li);
}
void ieinitnatsc(cov_model *cov, localinfotype *li) {
cov_model *next = cov->sub[0]->sub[0];
if ( CovList[next->nr].ieinit == NULL)
error("# cannot find ieinit -- please inform author");
CovList[next->nr].ieinit(next, li);
}
void tbm2natsc(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
cov_fct *C = CovList + next->nr;
double scale = CovList[next->nr].naturalscale(next),
y = x[0] / scale; // purely spatial ,
C->tbm2(&y, next, v);
}
int checknatsc(cov_model *cov) {
cov_model *next = cov->sub[0];
cov_fct *C = CovList + cov->nr;
int err;
strcpy(ERROR_LOC, C->name);
assert(cov->nr == NATSC);
cov->manipulating_x = true;
if ((err = check2X(next, cov->tsdim, cov->xdim, cov->statIn, cov->isoIn,
SUBMODELM))
!= NOERROR) {
return err;
}
if (next->sub[0]->statIn == cov->statIn &&
next->sub[0]->isoIn == cov->isoIn) {
next->tsdim = DEL_COV - 12;
}
// PrintModelInfo(cov);
// printf(" %d %d %d %d \n", next->sub[0]->statIn, cov->statIn,
// next->sub[0]->isoIn, cov->isoIn); // assert(false);
// now, next is shifted behind gatter !!:
if (next->nr >= GATTER && next->nr <= LASTGATTER) {
next = next->sub[0];
} else ERR("natsc: # expected -- contact author");
if (CovList[next->nr].naturalscale == NULL) {
sprintf(ERRORSTRING, "natural scaling is not defined for %s",
CovList[next->nr].name);
return ERRORFAILED;
}
cov->vdim = next->vdim;
setbackward(cov, next);
return NOERROR;
}
int initspectralnatsc(cov_model *cov) {
cov_model *next = cov->sub[0];
if (CovList[next->nr].initspectral == NULL) {
// PRINTF("natsc: %ld %s\n",
// (POINTER) CovList[next->nr].initspectral, CovList[next->nr].name);
PrintModelInfo(cov);
assert(CovList[next->nr].initspectral != NULL);
}
return CovList[next->nr].initspectral(next);
}
void spectralnatsc(cov_model *cov, spectral_storage *s, double *e) {
cov_model *next = cov->sub[0];
int d,
dim = cov->tsdim;
double
scale = CovList[next->nr].naturalscale(next);
CovList[next->nr].spectral(next, s, e);
for (d=0; d<dim; d++) e[d] /= scale;
}
void rangenatsc(cov_model *cov, range_arraytype* ra){ }
int initnatsc(method_type *meth){
cov_model *cov=meth->cov,
*next = cov->sub[0];
// globalparam *gp = meth->gp;
int err, i,
origdim = meth->loc->timespacedim,
total = origdim * meth->xdimout;
double
scale = CovList[next->nr].naturalscale(next);
// setAniso(meth); // schon im S-Teil ist die neue c-aniso gesetzt --
// S-Teil darf somit nicht mit c-aniso aufgerufen werden
// siehe doS, das sofort naechsten Teil aufruft.
meth->sub[0] = (method_type*) malloc(sizeof(method_type));
METHOD_NULL(meth->sub[0]);
meth->nsub = 1;
method_type *sub = meth->sub[0];
cpyMethod(meth, sub, false);
// cum_dollar(meth, meth->loc->timespacedim, cov, sub); :
if (meth->caniso != NULL) {
sub->caniso = (double *) malloc(total * sizeof(double));
for (i=0; i<total; i++) sub->caniso[i] = meth->caniso[i] / scale;
} else {
sub->cscale = meth->cscale / scale;
}
sub->cov = cov->sub[0];
err = initstandard(sub);
meth->compatible = sub->compatible;
meth->domethod = doS;
meth->S = NULL;
meth->nr = DOLLAR;
return err;
}
void donatsc(method_type *meth, res_type *res){
do_meth dometh = meth->sub[0]->domethod;
dometh(meth->sub[0], res);
}
////////////////////////////////////////////////////////////////////
void X(double *x, cov_model *cov, double *v){
cov_model *sub,
*next = cov->sub[0];
listoftype *list = (listoftype*) (cov->p[0]),
*sublist;
double var;
int i,
end=cov->vdim*cov->vdim;
if (LIST_ELEMENT < 0 || LIST_ELEMENT >= cov->nrow[0])
ERR("illegal call of 'X' (it may be used only within fitvario)");
if (cov->nrow[0] == 0) { // no X given
CovList[next->nr].cov(x, cov, v);
return;
} else {
if (!(bool*) cov->q) ERR("model 'const' expected"); // i.e. CONST
if (CovMatrixRow >= list->nrow[LIST_ELEMENT] ||
CovMatrixCol >= list->ncol[LIST_ELEMENT]) {
ERR("matrices are not sufficently large");
}
sub = next;
var = 1.0;
if (sub->nr >= GATTER && sub->nr <= LASTGATTER) sub=sub->sub[0];
if (sub->nr >= DOLLAR && sub->nr <= LASTDOLLAR) {
var = sub->p[DVAR][0];
sub=sub->sub[0];
if (sub->nr >= GATTER && sub->nr <= LASTGATTER) sub=sub->sub[0];
}
sublist = (listoftype*) (sub->p[0]);
if (list->ncol[LIST_ELEMENT] != sublist->nrow[LIST_ELEMENT])
ERR("mixed model: X and covariance matrix M do not match");
*v = XkCXtl(list->p[LIST_ELEMENT], sublist->p[LIST_ELEMENT],
list->nrow[LIST_ELEMENT], list->ncol[LIST_ELEMENT],
CovMatrixRow, CovMatrixCol);
for (i=0; i<end; v[i++] *= var);
}
}
void Xnonstat(double *x, double *y, cov_model *cov, double *v) {
if (cov->nsub != 0 && cov->nrow[0] == 0) { // no X given, but a submodel
cov_model *next = cov->sub[0];
CovList[next->nr].nonstat_cov(x, y, next, v);
return;
}
X(x, cov, v);
}
// ToDO : write simulation init, do !!!
int checkX(cov_model *cov) {
cov_model *sub,
*next = cov->sub[0];
int err,
nsub = cov->nsub,
*nrow = cov->nrow; // taken[MAX DIM],
// *ncol = cov->ncol,
// listoftype *list = (listoftype*) (cov->p[0]);
bool *Const;
CovMatrixRow = CovMatrixCol = INT_MAX;
cov->vdim = 1;
assert(nsub==1);
if (cov->q == NULL) cov->q = (double*) malloc(sizeof(bool));
Const = ((bool*) cov->q);
// PrintModelInfo(next);
sub = next;
if (sub->nr >= GATTER && sub->nr <= LASTGATTER) sub=sub->sub[0];
if (sub->nr >= DOLLAR && sub->nr <= LASTDOLLAR) {
sub=sub->sub[0];
if (sub->nr >= GATTER && sub->nr <= LASTGATTER) sub=sub->sub[0];
}
*Const = sub->nr == CONSTANT;
// printf("call: %s %s\n", CovList[next->nr].name, CovList[next->sub[0]->nr].name);
if ((err = check2X(next, cov->tsdim, cov->xdim, cov->statIn, cov->isoIn,
SUBMODELM)) != NOERROR) {
// printf("error\n");
return err;
}
if (*Const) {
next->tsdim = DEL_COV - 6;
if (next->nr >= DOLLAR && next->nr<=LASTDOLLAR)
next->sub[0]->tsdim = DEL_COV - 6;
// if (ncol[0]==0)
// ERR("if the cov. model is `const' in trend, then X must be given");
} else {
if (nrow[0]>0)
ERR("if the cov. model is not `const' in trend, then X may not be given");
}
cov->vdim=sub->vdim;
setbackward(cov, next);
if (cov->vdim > 1)
error("multivariate version not programmed yet");
// incorrect. but save !!
cov->semiseparatelast = false; // taken[xdim - 1] <= 1;
cov->separatelast = false; // taken[xdim - 1] <= 1; ??
return NOERROR;
}
void rangeX(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges + 0;
range->min[0] = RF_NEGINF;
range->max[0] = RF_INF;
range->pmin[0] = -1e10;
range->pmax[0] = 1e10;
range->openmin[0] = true;
range->openmax[0] = true;
range->finiterange = false;
}
void mixed(double *x, cov_model *cov, double *v){
ERR("do not use 'mixed' as a component for the error -- use model 'X' instead");
}
void mixed_nonstat(double *x, double *y, cov_model *cov, double *v){
mixed(x, cov, v);
}
void MLEmixed(double *x, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
int i, end = cov->vdim * cov->vdim;
if (cov->nsub == 0) {
for (i=0; i<end; i++) v[i] = 0.0;
} else {
// printf("mm %s\n", CovList[cov->nr].name);
CovList[next->nr].cov(x, next, v);
}
}
void MLEmixed_nonstat(double *x, double *y, cov_model *cov, double *v){
cov_model *next = cov->sub[0];
int i,
end = cov->vdim * cov->vdim;
if (cov->nsub == 0) {
for (i=0; i<end; i++) v[i] = 0.0;
} else {
CovList[next->nr].nonstat_cov(x, y, next, v);
}
}
int checkmixed(cov_model *cov) {
// cov_model *sub;
int i, err,
nsub = cov->nsub,
*ncol = cov->ncol,
*nrow = cov->nrow; // taken[MAX DIM],
char msg[50];
listoftype *list = (listoftype*) (cov->p[0]);
CovMatrixRow = CovMatrixCol = INT_MAX;
cov->vdim = 1;
if (ncol[1] > 0) { // b is given
if (ncol[0] == 0) ERR("if b is given X must be given");
if (nsub != 0) ERR("either b or a covariance model must be given in trend");
for (i=0; i<nrow[0]; i++) {
if (list->ncol[i] != nrow[1]) {
sprintf(msg,
"%dth set: (%d x %d) matrix X and (%d x %d) vector b do not match",
i, list->nrow[0], list->ncol[i], nrow[1], ncol[1]);
ERR(msg);
}
}
} else if (nsub == 0) { // only X is given -- then only a deterministic
// constant is given
if (ncol[1] == 0)
ERR("if no covariance model is given in mixed model then b must be given");
if (ncol[0] != 1) // deterministic effect
ERR("X must have one column");
} else { // submodel is given
cov_model *sub, *next = cov->sub[0];
bool *Const;
if (cov->q == NULL) cov->q = (double*) malloc(sizeof(bool));
Const = ((bool*) cov->q);
// PrintModelInfo(next);
sub = next;
if (sub->nr >= GATTER && sub->nr <= LASTGATTER) sub=sub->sub[0];
if (sub->nr >= DOLLAR && sub->nr <= LASTDOLLAR) {
sub=sub->sub[0];
if (sub->nr >= GATTER && sub->nr <= LASTGATTER) sub=sub->sub[0];
}
*Const = sub->nr == CONSTANT;
// printf("call: %s %s\n", CovList[next->nr].name, CovList[next->sub[0]->nr].name);
if ((err = check2X(next, cov->tsdim, cov->xdim, cov->statIn, cov->isoIn,
SUBMODELM)) != NOERROR) {
// printf("error\n");
return err;
}
if (*Const) {
next->tsdim = DEL_COV - 6;
if (next->nr >= DOLLAR && next->nr<=LASTDOLLAR)
next->sub[0]->tsdim = DEL_COV - 6;
// if (ncol[0]==0)
// ERR("if the cov. model is `const' in trend, then X must be given");
} else {
if (nrow[0]>0)
ERR("if the cov. model is not `const' in trend, then X may not be given");
}
cov->vdim=sub->vdim;
setbackward(cov, next);
}
if (cov->vdim > 1)
error("multivariate version not programmed yet");
// incorrect. but save !!
cov->semiseparatelast = false; // taken[xdim - 1] <= 1;
cov->separatelast = false; // taken[xdim - 1] <= 1; ??
return NOERROR;
}
void rangemixed(cov_model *cov, range_arraytype* ra){
range_type *range = ra->ranges + 0;
int i;
for (i=0; i<=1; i++) {
range->min[i] = RF_NEGINF;
range->max[i] = RF_INF;
range->pmin[i] = -1e10;
range->pmax[i] = 1e10;
range->openmin[i] = true;
range->openmax[i] = true;
}
range->finiterange = false;
}
typedef struct mixed_storage {
// bool genuine_dim[MAX362DIM];
key_type key;
double *b;
} mixed_storage;
void mixed_destruct(void ** S) {
if (*S!=NULL) {
mixed_storage *s = *((mixed_storage**) S);
if (s->b != NULL) free(s->b);
free(*S);
*S = NULL;
}
}
void mixed_NULL(mixed_storage* s) {
KEY_NULL(&(s->key));
s->b = NULL;
}
int initmixed(method_type *meth) {
return ERRORNOTPROGRAMMED;
/*
location_type *loc = meth->loc;
mixed_storage *s;
globalparam *gp = meth->gp;
cov_model *cov= meth->cov,
*next = cov->sub[0];
char errorloc_save[nErrorLoc];
int err;
// cholesky zerlegung + bereitstellung von b
SET_DESTRUCT(mixed_destruct);
meth->domethod = dotrend;
meth->compatible = true;
if (next==NULL) return NOERROR;
// submodel exists:
if ((meth->S = malloc(sizeof(mixed_storage)))==0){
err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
s = (mixed_storage*) meth->S;
mixed_NULL(s);
if ((s->b = malloc(sizeof(double) * loc->totalpoints)) == 0) {
err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
strcpy(errorloc_save, ERROR_LOC);
sprintf(ERROR_LOC, "%s%s: ", errorloc_save, METHODNAMES[method]);
err = internal_InitSimulateRF(loc->x, loc->T, loc->spatialdim,
loc->lx, loc->grid, loc->Time,
meth->simu.distribution, &(s->key),
&global, 1, next);
strcpy(ERROR_LOC, errorloc_save);
if (err == NOERROR) return;
ErrorHandling:
return err;
*/
}
void domixed(method_type *meth, res_type *res){
error("domixed not programmed yet");
/*
int i;
cov_model *cov = meth->cov,
*next = cov->sub[0];
double *b;
*/
}
