https://github.com/cran/RandomFields
Tip revision: 12173496c6b934b4bc9d81609e19178b9675310f authored by Martin Schlather on 25 July 2014, 00:00:00 UTC
version 3.0.30
version 3.0.30
Tip revision: 1217349
trend.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
Handling of the different possibilities to pass the trend
Note:
* Never use the below functions directly, but only by the functions indicated
in RFsimu.h, since there is no error check (e.g. initialization of RANDOM)
* VARIANCE, SCALE are not used here
* definitions for the random coin method can be found in MPPFcts.cc
* definitions for genuinely anisotropic or nondomain models are in
SophisticatedModel.cc; hyper models also in Hypermodel.cc
Copyright (C) 2011 - 2014 Marco Oesting & 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.
*/
/*Note:
the parameter 'polycoeff' is hidden in R; it is a vector consisting of the
coefficients of the corresponding trend polynomials:
the first choose(polydeg[1]+d,d) components belong to the first polynomial,
the next choose(polydeg[2]+d,d) to the second one,...
the corresponding monomial functions are of the following order:
1, x, x^2, ..., x^k, y, x y, ... x^(k-1) y, y^2, x y^2..., y^k,
z, x z, x^2 z, ...., x^(k-1) z, y z, x y z, x^2 y z, ..., z^k
*/
#include <R.h>
#include <Rdefines.h>
#include <Rinternals.h>
#include <math.h>
#include "RF.h"
#include "primitive.h"
#include <R_ext/Lapack.h>
#include <R_ext/Linpack.h>
#include "Covariance.h"
//////////////////////////////////////////////////////////////////////
// mixed
//////////////////////////////////////////////////////////////////////
#define MIXED_CONSTANT 0
void mixed(double *x, cov_model *cov, double *v) {
// bis auf den letzten Schluss alles schon auf bel vdim geschrieben.
cov_model *next = cov->sub[0];
location_type *loc = Loc(cov);
mixed_storage *s = cov->Smixed;
listoftype *X = PLIST(MIXED_X);
double //distTF,
var = 1.0,
*covmatrix=NULL;
int i, err, // Xnrow,
// err=NOERROR,
element = P0INT(MIXED_ELMNT),
vdim = cov->vdim2[0],
vdimsq = vdim * vdim;
if (cov->vdim2[0] != cov->vdim2[1]) BUG;
NotProgrammedYet("");
if (cov->nsub == 0) {
for (i=0; i<vdimsq; i++) v[i] = 0.0;
return;
}
if (X == NULL) { // no X given, but a submodel
COV(x, next, v);
return;
}
if (element < 0 || element >= cov->nrow[MIXED_X]) {
// realex element <0 ??
GERR1("Illegal call of the internal function '%s'", NICK(cov));
}
// PMI(cov->calling);
// crash(cov);
NotProgrammedYet(""); // hier irgendwo stehen geblieben
if (cov->q[MIXED_CONSTANT]) {
// Xu , u ~ N(0, C)
// classical random effect
cov_model *sub = next;
while (isDollar(sub)) {
var *= PARAM0(sub, DVAR);
sub=sub->sub[0];
}
listoftype *list= PARAMLIST(next, CONSTANT_M);
covmatrix = list->p[element];
if (X->ncol[element] != list->nrow[element]) {
GERR("mixed model: X and covariance matrix M do not match");
}
} else {
// Xu , u ~ random field
// random effect, geostatistical covariance model
CovarianceMatrix(cov->key, s->mixedcov);
covmatrix = s->mixedcov;
} // sub->nr != CONSTANT
// ab hier vdim =1
*v = XkCXtl(X->p[element], covmatrix,
X->nrow[element], X->ncol[element],
loc->i_row, loc->i_col);
*v *= var;
return;
ErrorHandling:
XERR(err);
}
void mixed_nonstat(double *x, double *y, cov_model *cov, double *v){
int
element = P0INT(MIXED_ELMNT);
NotProgrammedYet("");
if (element < 0 || element >= cov->nrow[MIXED_X]) {
// relax element <0 ??
char msg[100];
sprintf(msg, "Illegal call of the internal function '%s'",
NICK(cov));
error(msg);
}
if (cov->nsub != 0 && PisNULL(MIXED_X)) { // no X given, but a submodel
cov_model *next = cov->sub[0];
NONSTATCOV(x, y, next, v);
} else mixed(x, cov, v);
}
void covmatrix_mixed(cov_model *cov, double *v) {
cov_model *sub = cov->sub[0];
int element = P0INT(MIXED_ELMNT);
if (cov->ncol[MIXED_X] == 0 || element < 0) {
CovList[sub->nr].covmatrix(sub, v);
return;
}
if (element >= cov->nrow[MIXED_X]) ERR("value of 'element' is too large");
listoftype *X = PLIST(MIXED_X);
double *C=NULL;
int nrow = X->nrow[element],
ncol=X->ncol[element];
C = (double*) MALLOC(sizeof(double) * ncol * ncol);
if (C==NULL) {
StandardCovMatrix(cov, v);
return;
}
CovList[sub->nr].covmatrix(sub, C);
XCXt(X->p[element], C, v, nrow, ncol);
Loc(cov)->totalpoints = nrow;
// PMI(cov);
/*
int i,j,k, tot=Loc(cov)->totalpoints;
printf("\nStart mixed C %d ---- check within comment\n", element); //
for (k=i=0; i<ncol * ncol; i+=tot) {
for (j=0; j<ncol; j++) printf("%f ", C[k++]); //
printf("\n");//
}
printf("\nStart mixed t(X) %d\n", element); //
for (k=i=0; i<tot*tot; i+=tot) {
for (j=0; j<tot; j++) printf("%f ", X->p[element][k++]);//
printf("\n");//
}
printf("\nStart mixed t(v) %d\n", element); //
for (k=i=0; i<tot*tot; i+=tot) {
for (j=0; j<tot; j++) printf("%f ", v[k++]);//
printf("--- end check within comment \n");//
}
*/
free(C);
}
int set_mixed_constant(cov_model *cov) {
cov_model
*next = cov->sub[0],
*sub = next;
bool simple = true;
location_type *loc=Loc(cov);
int i,
totalpoints = loc->totalpoints,
*ncol = cov->ncol,
*nrow = cov->nrow;
listoftype *X = PLIST(MIXED_X);
cov->q = (double*) MALLOC(sizeof(double) * 1);
cov->qlen = 1;
// PMI(cov);
while (sub != NULL && isDollar(sub) &&
((simple = PARAMisNULL(sub, DPROJ)
&& (PARAMisNULL(sub, DSCALE) || PARAM0(sub, DSCALE) == 1.0)
&& PARAMisNULL(sub, DANISO)))) sub=sub->sub[0];
if ((cov->q[MIXED_CONSTANT] = sub != NULL && sub->nr == CONSTANT)) {
//next->delflag = DEL_COV - 6;
if (isDollar(next) && next->nrow[DVAR]==0) {
//next->delflag = DEL_COV - 6;
if (!simple)
SERR1("'%s' not allowed together with an anisotropic structrue",
NICK(cov));
}
int constant_size;
for (i=0; i<nrow[MIXED_X]; i++) {
constant_size = PARAMLIST(sub, CONSTANT_M)->nrow[i];
if (ncol[MIXED_X] > 0 && X->ncol[i] != constant_size) {
SERR5("%dth matrix 'X' (%d x %d) and (%d x %d) constant matrix 'M' do not match", i, X->nrow[i], X->ncol[i], constant_size, constant_size);
}
}
} else {
for (i=0; i<nrow[MIXED_X]; i++) {
if (X->nrow[i] != X->ncol[i])
SERR3("%dth matrix is not symmetric (%d x %d)",
i+1, X->nrow[i], X->ncol[i]);
}
}
if (false) // test ist wichtig, kollidiert z.Zt. aber mit SetAndGetModelInfo
for (i=0; i<nrow[MIXED_X]; i++) {
if (X->nrow[i] != totalpoints) {
//PMI(cov);
SERR3("number of rows of entry %d of 'X' (%d) are different from the number of locations (%d)", i+1, X->nrow[i], totalpoints);
}
}
return NOERROR;
}
char iscovmatrix_mixed(cov_model *cov) {
int err;
// printf("iscov %d %ld\n", cov->qlen, cov->q);
if (cov->q == NULL && (err = set_mixed_constant(cov)) != NOERROR) XERR(err);
// printf("iscov %d %ld\n", cov->qlen, cov->q);
return 2 * (int) (cov->nsub > 0) * (int) (cov->q[MIXED_CONSTANT] ||
cov->ncol[MIXED_X] > 0);
}
void kappamixed(int i, cov_model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc){
if (i==MIXED_DIM || i==MIXED_ELMNT) *nc = *nr = 1;
else if (i==MIXED_BETA || i==MIXED_DIST) {
*nc = 1;
*nr = SIZE_NOT_DETERMINED;
} else if (i==MIXED_X || i==MIXED_COORD) *nc = *nr = SIZE_NOT_DETERMINED;
else *nc = *nr = -1;
}
int checkmixed(cov_model *cov) {
// cov_model *sub;
//location_type *loc=Loc(cov);
int i, err,
nkappa = CovList[cov->nr].kappas,
//totalpoints = loc->totalpoints,
nsub = cov->nsub,
*ncol = cov->ncol,
*nrow = cov->nrow; // taken[MAX DIM],
char msg[300];
listoftype *X = PLIST(MIXED_X);
ROLE_ASSERT(ROLE_GAUSS || cov->role == ROLE_COV);
for (i=0; i<nkappa; i++) // NEVER; cf kriging.R, for instance
if (cov->kappasub[i] != NULL) SERR("parameters may not be random")
cov->vdim2[0] = cov->vdim2[1] = 1; //falls kein submodel vorhanden (Marco)
cov->maxdim=INFDIM;
cov->matrix_indep_of_x = true;
kdefault(cov, MIXED_ELMNT, 0);
if (ncol[MIXED_BETA] > 0) { // b is given
if (nsub != 0)
SERR("b and a covariance model may not be given at the same time");
if (ncol[MIXED_X] == 0) SERR("if 'b' is given 'X' must be given");
for (i=0; i<nrow[MIXED_X]; i++) {
if (X->ncol[i] != nrow[MIXED_BETA]) {
sprintf(msg,
"%dth set: (%d x %d) matrix X and (%d x %d) vector b do not match",
i, X->nrow[0], X->ncol[i], nrow[MIXED_BETA], ncol[MIXED_BETA]);
SERR(msg);
}
}
} else if (nsub == 0) { // only X is given -- then only a deterministic
// constant is given
if (ncol[MIXED_BETA] == 0)
SERR("if no covariance model is given then 'b' must be given");
if (ncol[MIXED_X] != 1) // deterministic effect
SERR("X must have one column");
kdefault(cov, MIXED_BETA, 1);
} else { // submodel is given
cov_model
*next = cov->sub[0];
// double var = 1.0;
if (cov->tsdim != cov->xdimprev || cov->tsdim != cov->xdimown)
return ERRORDIM;
if ((err = CHECK(next, cov->tsdim, cov->xdimown, PosDefType,
cov->domown,
cov->isoown, SUBMODEL_DEP, ROLE_COV)) != NOERROR) {
// print("error\n");
return err;
}
if (cov->q == NULL && (err = set_mixed_constant(cov)) != NOERROR) {
return err;
}
// warning("some checks in model `mixed' are missing");
// ob X mit C zusammengeht.
setbackward(cov, next);
}
if (cov->vdim2[0] > 1)
SERR("multivariate version of mixed not programmed yet");
if (PisNULL(MIXED_DIST) xor PisNULL(MIXED_DIM))
SERR("if 'dim' and 'dist' must be given at the same time");
if (PisNULL(MIXED_DIST) xor PisNULL(MIXED_COORD))
SERR("'dist' and 'coord' may not be given together");
// incorrect. but save !!
cov->semiseparatelast = false; // taken[tsxdim - 1] <= 1;
cov->separatelast = false; // taken[tsxdim - 1] <= 1; ??
return NOERROR;
}
void rangemixed(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){
int i;
range->min[MIXED_ELMNT] = 0;
range->max[MIXED_ELMNT] = MAXELEMENTS;
range->pmin[MIXED_ELMNT] = 0;
range->pmax[MIXED_ELMNT] = MAXELEMENTS;
range->openmin[MIXED_ELMNT] = false;
range->openmax[MIXED_ELMNT] = false;
for (i=MIXED_X; i<=MIXED_COORD; 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;
}
i=MIXED_DIST;
range->min[i] = 0;
range->max[i] = RF_INF;
range->pmin[i] = 1e-10;
range->pmax[i] = 1e10;
range->openmin[i] = false;
range->openmax[i] = true;
}
int initmixed(cov_model *cov, storage VARIABLE_IS_NOT_USED *S) {
location_type *loc = Loc(cov);
mixed_storage *s;
char errorloc_save[nErrorLoc];
double
*b = P(MIXED_BETA),
*coord = P(MIXED_COORD),
*dist = P(MIXED_DIST);
int
Cn = -1,
Cdim = -1,
list_element=0,
err = NOERROR,
vdim = cov->vdim2[0],
dim = coord!=NULL ? cov->ncol[MIXED_COORD] : P0INT(MIXED_DIM),
totalpoints = loc->totalpoints,
total = vdim * totalpoints;
listoftype *X = PLIST(MIXED_X);
bool distTF;
assert(cov->vdim2[0] == cov->vdim2[1]);
return ERRORFAILED; // muss in zerlegt werden in init und struct
// und unten struct aufrufen richtig praeparieren.
ROLE_ASSERT_GAUSS;
// cholesky zerlegung + bereitstellung von b
strcpy(errorloc_save, ERROR_LOC);
sprintf(ERROR_LOC, "%s%s: ", errorloc_save, "init mixed model");
assert(cov->nr = MIXEDEFFECT);
// submodel exists:
if ((cov->Smixed = (mixed_storage*) MALLOC(sizeof(mixed_storage)))==NULL){
err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
s = cov->Smixed;
MIXED_NULL(s);
NotProgrammedYet("");
if (s->Xb != NULL) free(s->Xb);
if (cov->ncol[MIXED_BETA] > 0) { // b is given
// X is given, but no covariance model
if (cov->nrow[MIXED_X] > 1) {
warning("using first element of list X only");
}
if ((s->Xb = (double*) MALLOC(sizeof(double) * X->nrow[list_element])) ==
NULL) { err=ERRORMEMORYALLOCATION; goto ErrorHandling; }
Ax(X->p[list_element], b, X->nrow[list_element], X->ncol[list_element],
s->Xb);
} else { // submodel is given
NotProgrammedYet(""); //if (cov->q[MIXED_CONSTANT]) error("not 'constant'").
cov_model *sub, *next = cov->sub[0];
if ((s->Xb = (res_type*) MALLOC(sizeof(res_type) * total)) == 0) {
err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
if ((distTF = dist!=NULL)) {
// Xu , u ~ geostatmodel, u i.a. viel laenger als X (tierzucht)
// stehende Vektoren, aneinandergereiht
Cn = cov->ncol[MIXED_DIST];
Cn = (int) (0.5 + 0.5 * (1.0 + sqrt(1 + 8.0 * Cn)));
Cdim = cov->nrow[MIXED_DIST];
assert(Cdim==1);
if ((err = covcpy(&(cov->key), next, dist, NULL, dim, dim, Cn,
false /*Time */,
false, distTF)) != NOERROR) goto ErrorHandling;
} else if (coord != NULL) {
// Xu , u ~ geostatmodel, u i.a. viel laenger als X (tierzucht)
// stehende Vektoren, aneinandergereiht
Cn = cov->ncol[MIXED_COORD];
Cdim = cov->nrow[MIXED_COORD];
if (Cdim > Cn)
warning("The dimension of the coordinates is higher than the number of points");
if ((err = covcpy(&(cov->key), next)) != NOERROR) goto ErrorHandling;
} else {
// Xu , u ~ geostatmodel, X quadratisch
Cn = loc->totalpoints;
if ((err = covcpy(&(cov->key), next, coord, NULL, dim, dim, Cn, false,
false, distTF)) != NOERROR) goto ErrorHandling;
}
if (cov->key->nr != GAUSSPROC) addModel(&(cov->key), GAUSSPROC);
cov->key->calling = cov->key;
cov->stor = (storage *) MALLOC(sizeof(storage));
STORAGE_NULL(cov->stor);
if ((err = INIT(cov->key, 0, cov->stor)) != NOERROR) goto ErrorHandling;
int Xnrow, Xncol;
Xnrow = Xncol = Cn * sub->vdim2[0];
if (loc->i_row==0 && loc->i_col==0) {
if (s->mixedcov != NULL) free(s->mixedcov);
s->mixedcov = (double*) MALLOC(sizeof(double) * Xnrow * Xnrow);
if (s->mixedcov == NULL) {
err = ERRORMEMORYALLOCATION;
goto ErrorHandling;
}
}
} // end of submodel
FieldReturn(cov);
ErrorHandling:
cov->initialised = err == NOERROR;
return err;
}
static int keeprandomeffect = false;
void domixed(cov_model *cov, storage VARIABLE_IS_NOT_USED *S){
location_type *loc = Loc(cov);
mixed_storage *s = cov->Smixed;
double *res = cov->rf;
int i,
list_element=0,
vdim = cov->vdim2[0],
totalpoints = loc->totalpoints,
total = vdim * totalpoints;
assert(cov->vdim2[0] == cov->vdim2[1]);
listoftype
*X = PLIST(MIXED_X);
if (cov->ncol[MIXED_BETA] > 0) { // b is given
// X is given, but no covariance model
if (total == X->nrow[list_element]) {
for (i=0; i<total; i++) res[i] = s->Xb[i];
} else {
assert(X->nrow[list_element] == 1);
for (i=0; i<total; i++) res[i] = s->Xb[0];
}
// print("%f\n", res[0]); assert(false);
} else { // submodel is given
cov_model *key = cov->key;
if (!keeprandomeffect || !s->initialized) {
do_gaussprocess(key, cov->stor);
}
if (X != NULL) {
AxResType(X->p[list_element], cov->key->rf, X->nrow[list_element],
X->ncol[list_element], res);
} else {
res_type *rf = cov->key->rf;
for (i=0; i<total; i++) res[i] = rf[i];
}
}
}
//////////////////////////////////////////////////////////////////////
// trend
//////////////////////////////////////////////////////////////////////
void trend(double VARIABLE_IS_NOT_USED *x, cov_model *cov, double *v){
// todo : ueberlegen ob hier nicht der Trend-Wert berechnet werden sollte
///BUG; // sollte durch plusStat & plusNonStat abgeblock sein ?
// to do: stattdessen sollte hier der Trendwert zurueckgeliefert werden ?
// nein: model kann auch nur aus Trend bestehen!
int i,
vdim = cov->vdim2[0],
vSq = vdim * vdim;
for (i=0; i<vSq; i++) v[i]=0.0;
//print("trend %d %f\n", vSq, *v);
}
void trend_nonstat(double VARIABLE_IS_NOT_USED *x, double VARIABLE_IS_NOT_USED *y, cov_model *cov, double *v){
int i,
vdim = cov->vdim2[0],
vSq = vdim * vdim;
if (cov->role == ROLE_COV) for (i=0; i<vSq; i++) v[i]=0.0;
else ERR("trend is called unexpectately.");
}
void kappatrend(int i, cov_model *cov, int *nr, int *nc){
// i nummer des parameters
int k;
switch(i) {
case TREND_MEAN: //mu
*nr = SIZE_NOT_DETERMINED;
*nc = 1;
break;
case TREND_LINEAR: //plane
*nr = cov->tsdim;
*nc = SIZE_NOT_DETERMINED;
break;
case TREND_POLY: //polydeg
*nr = SIZE_NOT_DETERMINED;
*nc = 1;
break;
case TREND_PARAM_POLY: //polycoeff
if(PisNULL(TREND_POLY)) {
*nr = -1;
} else {
*nr = 0;
//APMI(cov);
for(k=0; k < cov->nrow[TREND_POLY]; k++) {
//printf("%d\n", *nr);
*nr += binomialcoeff(PINT(TREND_POLY)[k] +
cov->tsdim, cov->tsdim);
}
}
*nc = 1;
break;
case TREND_FCT: //arbitraryfct
*nr = 1;
*nc = 1;
break;
case TREND_PARAM_FCT: //fctcoeff
*nr = 1;
*nc = 1;
break;
default:
*nr = -1;
*nc = -1;
}
}
int checktrend(cov_model *cov){
// if (cov->ncol[TREND_LINEAR] > 0 || cov->ncol[TREND_FCT]>0
// || cov->ncol[TREND_PARAM_FCT] > 0)
// return(ERRORNOTPROGRAMMED);
double *mu = P(TREND_MEAN),
*plane = P(TREND_LINEAR);
int i,
nkappa = CovList[cov->nr].kappas,
*polydeg = PINT(TREND_POLY),
vdim = 0,
tsdim= cov->tsdim,
basislen = 0;
//SEXP Rx, fctbody, envir;
ROLE_ASSERT(ROLE_GAUSS || cov->role == ROLE_COV);
for (i=0; i<nkappa; i++) // NEVER; cf kriging.R, for instance
if (cov->kappasub[i] != NULL) SERR("parameters may not be random")
if (mu != NULL) {
vdim = cov->nrow[TREND_MEAN];
cov->matrix_indep_of_x = true;
}
if (plane != NULL) {
if(vdim>0 && (vdim != cov->ncol[TREND_LINEAR])) {
SERR("trend parameters have different multivariate dimensions");
} else vdim = cov->ncol[TREND_LINEAR];
cov->matrix_indep_of_x = false;
}
if (!PisNULL(TREND_POLY)) {
if (vdim>0) {
if (vdim != cov->nrow[TREND_POLY])
SERR("trend parameters have different multivariate dimensions");
SERR("polynomials and free functions in trend may not be mixed with other trend definitions. Please use a sum of trends.");
}
vdim = cov->nrow[TREND_POLY];
if (PisNULL(TREND_PARAM_POLY)) {
for (i=0; i<vdim; i++)
basislen += binomialcoeff(polydeg[i] + tsdim, tsdim);
PALLOC(TREND_PARAM_POLY, basislen, 1);
for(i=0; i<basislen; i++) P(TREND_PARAM_POLY)[i] = NA_REAL;
}
cov->matrix_indep_of_x = false;
}
if (!PisNULL(TREND_FCT)) { //brauche hier: construct.fct$vdim
cov->matrix_indep_of_x = false;
NotProgrammedYet("arbitrary function");
// kdefault(cov, TREND_PARAM_FCT, 1.0);
// PROTECT(envir = allocSExp(ENVSXP));
// SET_ENCLOS(envir, R_GlobalEnv);
//
// PROTECT(fctbody = allocSExp(CLOSXP));
// fctbody = BODY(*((SEXP *) arbitraryfct));
//
// PROTECT(Rx = allocVector(REALSXP,xlen));
// for(j=0; j<xlen; j++) REAL(Rx)[j] = 0;
// defineVar(install("x"), Rx, envir);
// defineVar(install("y"), ScalarReal(0), envir);
// defineVar(install("z"), ScalarReal(0), envir);
// defineVar(install("T"), ScalarReal(0), envir);
//
// int vdimproposal = length(eval(fctbody, envir));
// UNPROTECT(3);
//
// if (vdim > 0) {
// if (vdim != vdimproposal)
// GERR("trend parameters have different multivariate dimensions");
// GERR("polynomials and free functions in trend may not be mixed with other trend definitions. Please use a sum of trends.");
// }
// vdim = vdimproposal;
}
if (vdim <= 0) {
vdim = cov->calling->vdim2[0];
if (vdim <= 0)
SERR("multivariate dimension for trend cannot be determined.");
PALLOC(TREND_MEAN, vdim, 1);
for(i=0; i<vdim; i++) P(TREND_MEAN)[i] = 0.0;
}
cov->vdim2[0] = cov->vdim2[1] = vdim;
cov->isoown = cov->matrix_indep_of_x ? ISOTROPIC : CARTESIAN_COORD;
return NOERROR;
}
int checktrendproc(cov_model *cov){
return checktrend(cov);
}
void rangetrend(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){
//P(TREND_MEAN]: mu / mean
range->min[TREND_MEAN] = RF_NEGINF;
range->max[TREND_MEAN] = RF_INF;
range->pmin[TREND_MEAN] = -10^10;
range->pmax[TREND_MEAN] = 10^10;
range->openmin[TREND_MEAN] = true;
range->openmax[TREND_MEAN] = true;
//P(TREND_LINEAR]: plane
range->min[TREND_LINEAR] = RF_NEGINF;
range->max[TREND_LINEAR] = RF_INF;
range->pmin[TREND_LINEAR] = -10^10;
range->pmax[TREND_LINEAR] = 10^10;
range->openmin[TREND_LINEAR] = true;
range->openmax[TREND_LINEAR] = true;
//P(TREND_POLY]: polydeg / polynomial degree
range->min[TREND_POLY] = 0;
range->max[TREND_POLY] = RF_INF;
range->pmin[TREND_POLY] = 0;
range->pmax[TREND_POLY] = 10;
range->openmin[TREND_POLY] = false;
range->openmax[TREND_POLY] = false;
//P(TREND_PARAM_POLY]: polycoeff / coefficients of polynomial
range->min[TREND_PARAM_POLY] = RF_NEGINF;
range->max[TREND_PARAM_POLY] = RF_INF;
range->pmin[TREND_PARAM_POLY] = -10^10;
range->pmax[TREND_PARAM_POLY] = 10^10;
range->openmin[TREND_PARAM_POLY] = true;
range->openmax[TREND_PARAM_POLY] = true;
//P(TREND_FCT]: arbitraryfct / arbitrary function
range->min[TREND_FCT] = RF_NEGINF;
range->max[TREND_FCT] = RF_INF;
range->pmin[TREND_FCT] = -10^10;
range->pmax[TREND_FCT] = 10^10;
range->openmin[TREND_FCT] = true;
range->openmax[TREND_FCT] = true;
//P(TREND_PARAM_FCT]: fctcoeff / coefficient of arbitrary function
range->min[TREND_PARAM_FCT] = RF_NEGINF;
range->max[TREND_PARAM_FCT] = RF_INF;
range->pmin[TREND_PARAM_FCT] = -10^10;
range->pmax[TREND_PARAM_FCT] = 10^10;
range->openmin[TREND_PARAM_FCT] = true;
range->openmax[TREND_PARAM_FCT] = true;
}
void GetInternalMeanI(cov_model *cov, int vdim, double *mean){
// assuming that mean[i]=0.0 originally for all i
int i;
if (cov->nr == TREND) {
if (cov->ncol[TREND_MEAN]==1) {
if (cov->nrow[TREND_MEAN] != vdim) {
for (i=0; i<vdim; i++) mean[i] = RF_NA;
return; // only scalar allowed !
}
for (i=0; i<vdim; i++) mean[i] += P(TREND_MEAN)[i];
}
}
if (cov->nr == PLUS || cov->nr == TREND) {
for (i=0; i<cov->nsub; i++) GetInternalMeanI(cov->sub[i], vdim, mean);
}
}
void GetInternalMean(cov_model *cov, int vdim, double *mean){
int i;
for (i=0; i<vdim; i++) mean[i]=0.0;
GetInternalMeanI(cov, vdim, mean);
}
int init_trend(cov_model *cov, storage *S) {
long err = NOERROR;
trend_storage *s;
int i,
*polydeg = PINT(TREND_POLY),
basislen=0,
tsdim = cov->tsdim,
vdim = cov->vdim2[0];
//SEXP fctformals, argnames;
if (cov->vdim2[0] != cov->vdim2[1]) BUG;
//assert(false);
ROLE_ASSERT_GAUSS;
if(polydeg != NULL) {
for(i=0; i<vdim; i++)
basislen += binomialcoeff(polydeg[i]+tsdim,tsdim);
}
if (cov->Strend != NULL) free(cov->Strend);
if ((cov->Strend = (trend_storage *) MALLOC(sizeof(trend_storage)))==NULL) {
err=ERRORMEMORYALLOCATION;
goto ErrorHandling;
}
s = cov->Strend;
TREND_NULL(s);
if ((s->x = (double *) MALLOC(sizeof(double) * tsdim))==NULL ||
(s->xi = (int *) MALLOC(sizeof(int) * tsdim))==NULL ||
(s->evalplane = (double *) MALLOC(sizeof(double) * vdim))==NULL ||
(basislen > 0 &&
(s->powmatrix = (int *) MALLOC(sizeof(int) * basislen * tsdim))==NULL)){
err=ERRORMEMORYALLOCATION;
goto ErrorHandling;
}
if (basislen > 0) poly_basis(cov, S); //generates basis of monomials
//each row consists of one basis element
//the j-th column consists of the power of the j-th space-time-dimension
if (!PisNULL(TREND_FCT)) { //hier werden Argumente von arbitraryfct ueberprueft
NotProgrammedYet("");
// fctformals = getAttrib(FORMALS(*((SEXP *) arbitraryfct)), R_NamesSymbol);
// nargs = length(fctformals);
// PROTECT(argnames = allocVector(STRSXP,1));
//
// SET_STRING_ELT(argnames,0,mkChar("y"));
// matchind = 0;
// for(j=0; j<nargs; j++)
// matchind += (STRING_ELT(argnames,0) == STRING_ELT(fctformals,j));
// if (matchind>0) {
// if((meth->loc->spatialdim < 2) || (meth->loc->xvectorvalued))
// GERR("The variable y does not match to the locations.\n");
// }
//
// SET_STRING_ELT(argnames,0,mkChar("z"));
// matchind = 0;
// for(j=0; j<nargs; j++)
// matchind += (STRING_ELT(argnames,0) == STRING_ELT(fctformals,j));
// if (matchind>0) {
// if((meth->loc->spatialdim < 3) || (meth->loc->xvectorvalued))
// GERR("The variable z does not match to the locations.\n")
// }
//
// SET_STRING_ELT(argnames,0,mkChar("T"));
// matchind = 0;
// for(j=0; j<nargs; j++)
// matchind += (STRING_ELT(argnames,0) == STRING_ELT(fctformals,j));
// if (matchind>0) {
// if (meth->loc->Time == false)
// GERR("The variable T may be used for time only.\n")
// }
// UNPROTECT(1);
}
err = FieldReturn(cov);
return NOERROR;
ErrorHandling:
return err;
}
void do_trend(cov_model *cov, storage VARIABLE_IS_NOT_USED *s){
location_type *loc = Loc(cov);
char errorloc_save[nErrorLoc];
trend_storage *S = cov->Strend;
//PMI(cov->calling->calling);
assert(S != NULL);
double t,
*mu = P(TREND_MEAN),
*plane = P(TREND_LINEAR),
*polycoeff = P(TREND_PARAM_POLY),
// *fctcoeff = P(TREND_PARAM_FCT),
**xgr = loc->xgr,
*x = S->x,
*evalplane = S->evalplane;
int i, j, k, v, w,
basislen, startindex,
*polydeg = PINT(TREND_POLY),
totalpoints = loc->totalpoints,
vdim = cov->vdim2[0],
tsdim = cov->tsdim,
spatialdim = loc->spatialdim,
*len = loc->length,
total = totalpoints * vdim,
*xi = S->xi,
*powmatrix = S->powmatrix;
//SEXP fctbody, tempres, envir, Rx;
res_type *res = cov->rf;
assert(cov->vdim2[0] == cov->vdim2[1]);
strcpy(errorloc_save, ERROR_LOC);
sprintf(ERROR_LOC, "%s%s: ", errorloc_save, "add trend model");
// print("%s\n", ERROR_LOC);
if (mu != NULL) {
for (k=0; k<total; ) {
for (v=0; v<vdim; res[k++] = mu[v++]);
}
} else {
for (k=0; k<total; res[k++]=0.0);
}
if (plane != NULL) {
if (loc->grid) {
for (w=0; w<tsdim; w++) {
x[w]=xgr[w][XSTART];
xi[w]=0;
}
int endfor = totalpoints * vdim;
for (k=0; k<endfor; ) {
xA(x, plane, cov->nrow[TREND_LINEAR], cov->ncol[TREND_LINEAR],
evalplane);
for(v=0; v<vdim; v++) res[k++] += evalplane[v];
i = 0;
(xi[i])++;
x[i] += xgr[i][XSTEP];
while(xi[i]>=len[i]) {
xi[i] = 0;
x[i] = xgr[i][XSTART];
if (i<tsdim-1) {
i++;
(xi[i])++;
x[i] += xgr[i][XSTEP];
} else {
assert(k==endfor);
}
}
}
} else if (loc->Time) {
int m,
endfor= loc->length[loc->timespacedim-1];
for(k=m=0, t=loc->T[XSTART]; m<endfor; m++, t+=loc->T[XSTEP]) {
// for(t=loc->T[XSTART], i=0; t < loc->T[XEND]; t += loc->T[XSTEP]) {
for(m=0, j=0; m < loc->spatialtotalpoints; m++) {
for(w=0; w<spatialdim; w++) x[w] = (loc->x)[j++];
x[spatialdim] = t;
xA(x, plane, cov->nrow[TREND_LINEAR], cov->ncol[TREND_LINEAR],
evalplane);
for(v=0;v<vdim;v++) res[k++] += evalplane[v];
}
}
} else {
int m,
endfor = totalpoints * tsdim;
for (k=m=0; m<endfor; m+=tsdim) {
xA(loc->x + m, plane, cov->nrow[TREND_LINEAR],
cov->ncol[TREND_LINEAR], evalplane);
for(v=0;v<vdim;v++) res[k++] += evalplane[v];
}
}
}
if(polydeg != NULL) {
startindex=0;
int end_k = totalpoints * vdim;
for(v=0; v<vdim; v++) {
basislen = binomialcoeff(polydeg[v] + tsdim, tsdim);
if(isnan(polycoeff[0])) {
ERR("Error: cannot evaluate polynomial without coefficients.\n");
}
if (loc->grid) {
for (w=0; w<tsdim; w++) {
x[w]=xgr[w][XSTART];
xi[w]=0;
}
for(k=v; k<end_k; k+=vdim) {
//evaluation of trend polynomial
res[k] += evalpoly(x, powmatrix + startindex*tsdim,
polycoeff + startindex, basislen, tsdim);
i = 0;
(xi[i])++;
x[i] += xgr[i][XSTEP];
while(xi[i]>=len[i]) {
xi[i] = 0;
x[i] = xgr[i][XSTART];
if (i<tsdim-1) {
i++;
(xi[i])++;
x[i] += xgr[i][XSTEP];
} else {
assert(k==v + end_k - vdim);
}
}
}
} else if(loc->Time) {
for(t=loc->T[XSTART], k=v; k<end_k; t+=loc->T[XSTEP]) {
for(i=0, j=0; i<loc->spatialtotalpoints; i++, k+=vdim) {
for(w=0; w<spatialdim; w++) x[w] = (loc->x)[j++];
x[spatialdim] = t;
//evaluation of trend polynomial
res[k] += evalpoly(x, powmatrix + startindex*tsdim,
polycoeff + startindex, basislen, tsdim);
}
}
} else {
for(k=v; k<end_k; k+=vdim) {
//evaluation of trend polynomial
res[k] += evalpoly(x, powmatrix + startindex*tsdim,
polycoeff + startindex, basislen, tsdim);
}
}
startindex += basislen;
}
}
if (!PisNULL(TREND_FCT)) { //muss hier arbitraryfct auswerten
NotProgrammedYet("");
// if (isnan(fctcoeff[0])) {
// ERR("Error: cannot evaluate function without coefficient.\n");
// }
//
// PROTECT(envir = allocSExp(ENVSXP));
// SET_ENCLOS(envir, R_GlobalEnv);
//
// PROTECT(fctbody = allocSExp(CLOSXP));
// fctbody = BODY(*((SEXP *) arbitraryfct));
//
// PROTECT(tempres = allocVector(REALSXP, vdim));
// PROTECT(Rx = allocVector(REALSXP,xlen));
//
// if (loc->grid) {
// for (w=0; w<tsdim; w++) {
// x[w]=xgr[w][XSTART];
// xi[w]=0;
// }
// for(k=0; k<totalpoints; k++) {
// //evaluation of trend polynomial
// for(w=0; w<xlen; w++) REAL(Rx)[w] = x[w];
// defineVar(install("x"), Rx, envir);
// if (spatialdim>1)
// defineVar(install("y"), ScalarReal(x[1]), envir);
// if (spatialdim>2)
// defineVar(install("z"), ScalarReal(x[2]), envir);
// defineVar(install("T"), ScalarReal(x[tsdim-1]), envir);
// tempres = eval(fctbody, envir);
// for (v=0; v<vdim; v++)
// res[k*vdim+v] += fctcoeff[0]*REAL(tempres)[v];
// i = 0;
// (xi[i])++;
// x[i] += xgr[i][XSTEP];
// while(xi[i]>=len[i]) {
// xi[i] = 0;
// x[i] = xgr[i][XSTART];
// if (i<tsdim-1) {
// i++;
// (xi[i])++;
// x[i] += xgr[i][XSTEP];
// } else {
// assert(k==totalpoints-1);
// }
// }
// }
// } else if(loc->Time) {
// int k,
// endfor= loc->length[loc->timespacedim-1];
// for(k=0, t=loc->T[XSTART], i=0; k<endfor; k++, t += loc->T[XSTEP]) {
// // for(t=loc->T[XSTART], i=0; t < loc->T[XEND]; t+=loc->T[XSTEP]) {
// for(k=0, j=0; k < loc->spatialtotalpoints; k++, i++) {
// for(w=0; w<spatialdim; w++) x[w] = (loc->x)[j++];
// x[spatialdim] = t;
// //evaluation of trend polynomial
// for(w=0; w<xlen; w++) REAL(Rx)[w] = x[w];
// defineVar(install("x"), Rx, envir);
// if (spatialdim > 1)
// defineVar(install("y"), ScalarReal(x[1]), envir);
// if (spatialdim > 2)
// defineVar(install("z"), ScalarReal(x[2]), envir);
// if (loc->Time)
// defineVar(install("T"), ScalarReal(x[tsdim-1]), envir);
// tempres = eval(fctbody, envir);
// for (v=0; v<vdim; v++)
// res[i*vdim+v] += fctcoeff[0]*REAL(tempres)[v];
// }
// }
// } else {
// for(k=0; k<totalpoints; k++) {
// //evaluation of trend polynomial
// for(w=0; w<spatialdim; w++) x[w] = (loc->x)[k*spatialdim+w];
// for(w=0; w<xlen; w++) REAL(Rx)[w] = x[w];
// defineVar(install("x"), Rx, envir);
// if (spatialdim>1)
// defineVar(install("y"), ScalarReal(x[1]), envir);
// if (spatialdim>2)
// defineVar(install("z"), ScalarReal(x[2]), envir);
// defineVar(install("T"), ScalarReal(x[tsdim-1]), envir);
// tempres = eval(fctbody, envir);
// for (v=0; v<vdim; v++)
// res[k*vdim+v] += fctcoeff[0]*REAL(tempres)[v];
// }
// }
// UNPROTECT(4);
}
return;
}
int binomialcoeff(int n, int k) {
//programmed as in wikipedia
int i, res=1;
if((k < 0) || (k > n)) return 0;
if(k > n-k) k = n-k; //symmetry
for(i=0; i<k; i++) {
res *= n-i;
res /= i+1; //two steps because of integer division
}
return res;
}
void poly_basis(cov_model *cov, storage VARIABLE_IS_NOT_USED *s) {
trend_storage *S = cov->Strend;
int basislen=0, powsum, d, i, j, k, v,
dim = cov->tsdim,
vdim = cov->vdim2[0],
*powmatrix = S->powmatrix,
*dimi=NULL,
err=NOERROR;
int *polydeg = PINT(TREND_POLY);
assert(cov->vdim2[0] == cov->vdim2[1]);
dimi = (int *) MALLOC(dim * sizeof(int));
if (dimi == NULL) {
err = ERRORMEMORYALLOCATION;
goto ErrorHandling;
}
for(v=0, j=0; v<vdim; v++) {
basislen = binomialcoeff(polydeg[v] + dim, dim);
// print("v: %d\n", v);
for(d=0; d<dim; d++) dimi[d] = 0;
for(k=0; k<basislen; k++) {
for(d=0; d<dim; d++) powmatrix[j++] = dimi[d];
// for(d=0;d<dim;d++) print("%d ", powmatrix[j*dim+d]);
// print("\n");
i=0;
(dimi[i])++;
powsum = 0;
for(d=0; d<dim; d++) powsum += dimi[d];
while(powsum > polydeg[v]) {
dimi[i] = 0;
if (i < dim-1) {
i++;
(dimi[i])++;
}
powsum = 0;
for(d=0; d<dim; d++) powsum += dimi[d];
}
}
// print("\n");
}
//print("\n");
ErrorHandling:
if (dimi != NULL) free(dimi);
dimi = NULL;
if (err != NOERROR) XERR(err);
return;
}
double evalpoly(double *x, int *powmatrix, double *polycoeff, int basislen,
int dim) {
int i, d, j;
double res = 0, tempres;
for(j=i=0; i<basislen; i++) {
tempres=1;
for(d=0; d<dim; d++) tempres *= pow(x[d], powmatrix[j++]);
res += polycoeff[i] * tempres;
}
return(res);
}
void likelihood_trend(double VARIABLE_IS_NOT_USED *x, double VARIABLE_IS_NOT_USED invmatrix, cov_model *cov,
double VARIABLE_IS_NOT_USED *v){
if (cov->role == ROLE_GAUSS) {
NotProgrammedYet("");
//C = covarianzmatrix
// b = (X^T C^{-1} X)^{-1} X^top C^{-1} y
// Baysiean mit b ~ N(b_0, C_B) :
// b = (X^T C^{-1} X + C_b^{-1})^{-1} [X^top C^{-1} y + C_b^{-1} b_0
} else {
// deterministic is OK
NotProgrammedYet("");
}
}
