RFsimu.cc
/*
Authors
Yindeng, Jiang, jiangyindeng@gmail.com
Martin Schlather, martin.schlather@math.uni-goettingen.de
library for unconditional simulation of stationary and isotropic random fields
Copyright (C) 2001 -- 2003 Martin Schlather
Copyright (C) 2004 -- 2004 Yindeng Jiang & Martin Schlather
Copyright (C) 2005 -- 2006 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.
*/
/*
PRINTING LEVELS
===============
error messages: 1
forcation : 2
minor tracing information : 3--5
large debugging information: >10
*/
/*
calculate the transformation of the points only once and store the result in
a register (indexed by ncov) of key, not of s->. Advantages: points have to
calculated only once for several tries; if zonal anisotropy then specific
methods may be called;
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
//#include <sys/timeb.h>
#include <assert.h>
#include <string.h>
#include "RFsimu.h"
#include "RFCovFcts.h"
#include <unistd.h>
int FirstCheck_Cov(key_type *key, int m, bool MultiplyAndHyper)
{
int v, error;
SimulationType Method;
covinfo_type *kc;
methodvalue_type *meth;
meth = &(key->meth[m]);
Method = meth->unimeth;
meth->actcov=0;
error = NOERROR;
for (v=0; v<key->ncov; v++) {
ERRORMODELNUMBER = v;
cov_fct *cov;
kc = &(key->cov[v]);
cov = &(CovList[kc->nr]);
meth->covlist[meth->actcov] = v;
// printf("%d %d %d %d\n", v, kc->method, Method, (int) kc->left);
if (kc->method==Method && kc->left) {
assert((kc->nr >= 0) && (kc->nr < currentNrCov));
assert(kc->param[VARIANCE] >= 0.0);
if (cov->implemented[Method] != IMPLEMENTED) return ERRORNOTDEFINED;
if (cov->type==ISOHYPERMODEL || cov->type==ANISOHYPERMODEL) {
if (!MultiplyAndHyper) error = ERRORHYPERNOTALLOWED;
else {
error =
cov->checkNinit(kc, &(COVLISTALL[v]), key->ncov - v - 1, Method,
key->anisotropy);
}
} else {
error = cov->check(kc->param, kc->reduceddim, kc->dim, Method);
}
if (error != NOERROR) return error;
// printf("v=%d %d %d \n", v, key->ncov, key->cov[v].op);
if (v < key->ncov -1 && key->cov[v].op) {
if (key->cov[v+1].method != Method) {
if (MultiplyAndHyper) {
if (GENERAL_PRINTLEVEL>0)
PRINTF("severe error -contact author. %d %d %d %d (%s) %d (%s)n",
v, key->cov[v-1].op, key->ncov, key->cov[v-1].method,
METHODNAMES[key->cov[v-1].method],
Method, METHODNAMES[Method]);
return ERRORMETHODMIX;
} else {
PRINTF("severe error - contact author. %d %d %d %d (%s) %d (%s)n",
v, key->cov[v-1].op, key->ncov, key->cov[v-1],
METHODNAMES[key->cov[v-1].method],
Method, METHODNAMES[Method]);
assert(false);
}
}
if (!MultiplyAndHyper) return ERRORNOMULTIPLICATION;
} // meth->actcov > 0
int w, nsub;
nsub = (cov->type==ISOHYPERMODEL || cov->type==ANISOHYPERMODEL) ?
(int) kc->param[HYPERNR] : 0;
for (w=0; w<nsub; w++) {
kc->left = false;
kc++;
(meth->actcov)++;
v++;
meth->covlist[meth->actcov] = v;
}
kc->left = false;
(meth->actcov)++;
} // left
}
ERRORMODELNUMBER = -1;
return (meth->actcov==0
? (error==NOERROR ? NOERROR_ENDOFLIST : error)
: (error==NOERROR ? NOERROR : NOERROR_REPEAT));
}
void COVINFO_NULL(covinfo_type *cov){
int m, d;
double *param;
covinfo_type *kc;
for (m=0; m<MAXCOV; m++) {
kc = &(cov[m]);
kc->x = NULL;
kc->reduceddim = kc->op = 0;
kc->nr = currentNrCov; // undefined model
kc->method = Forbidden;
kc->genuine_last_dimension = kc->simugrid = false;
kc->left = true;
param = kc->param;
for (d=0; d<MAXDIM; d++) {
kc->genuine_dim[d] = false;
kc->length[d] = kc->idx[d] = -1;
}
for (d=MAXDIM * 3 - 1; d>=0; d--) kc->xsimugr[d] = RF_NAN;
for (d=0; d<MAXDIMSQ; d++) param[d] = kc->aniso[d] = RF_NAN;
for (; d<TOTAL_PARAM; param[d++]=RF_NAN);
}
}
void KEY_NULL(key_type *key)
{
int m, d;
key->stop = key->active = false;
key->ncov = key->n_unimeth = 0;
key->spatialdim = key->timespacedim = key->totalparam = -1;
key->totalpoints = key->spatialtotalpoints = 0;
for (m=0; m<MAXCOV; m++) {
key->meth[m].S = NULL;
key->meth[m].destruct = NULL;
}
COVINFO_NULL(key->cov);
for (d=0; d<MAXDIM; d++) {
key->x[d]=NULL;
key->length[d] = -1;
}
key->T[0] = key->T[1] = key->T[2] = 0.0;
// key->naturalscaling =-1;
key->TrendModus = -1;
key->TrendFunction = NULL;
key->LinearTrend = NULL;
key->mean = 0.0;
}
void DeleteKeyNotTrend(key_type *key)
{
int d, m;
if (key->x[0]!=NULL) free(key->x[0]);
for (d=0; d<MAXDIM; d++) {key->x[d]=NULL;}
for (m=0; m<MAXCOV; m++) {
if (key->cov[m].x != NULL) {
free(key->cov[m].x);
key->cov[m].x = NULL;
}
if (key->meth[m].destruct!=NULL) {
key->meth[m].destruct(&(key->meth[m].S));
key->meth[m].destruct=NULL;
}
}
key->active = false;
key->ncov = key->n_unimeth = 0;
key->spatialdim = key->timespacedim = key->totalparam = -1;
key->totalpoints = key->spatialtotalpoints = 0;
}
void DeleteKeyTrend(key_type *key)
{
key->TrendModus = -1;
if (key->TrendFunction!=NULL) free(key->TrendFunction);
key->TrendFunction=NULL;
if (key->LinearTrend!=NULL) free(key->LinearTrend);
key->LinearTrend=NULL;
key->lTrendFct = 0;
key->lLinTrend = 0;
}
void DeleteKey(int *keyNr)
{
if (GENERAL_PRINTLEVEL>=4) {
PRINTF("deleting stored parameters of register %d ...\n", *keyNr);
}
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {return;}
key_type *key = &(KEY[*keyNr]);
DeleteKeyNotTrend(key);
DeleteKeyTrend(key);
key->active=false;
}
void DeleteAllKeys() {
int i;
for (i=0;i<MAXKEYS;i++) {DeleteKey(&i);}
}
void printkey(key_type *key) {
int i,j;
covinfo_type *keycov;
methodvalue_type *meth;
PRINTF("\ngrid=%d, active=%d, anisotropy=%d, sp_dim=%d,\ntotalpoints=%d, distr=%s Time=%d [%f %f %f]\ntimespacedim=%d compatible=%d ncov=%d\n",
key->grid, key->active, key->anisotropy,
key->spatialdim, key->totalpoints,
DISTRNAMES[key->distribution],
key->Time, key->T[0], key->T[1], key->T[2],
key->timespacedim, key->compatible,
key->ncov);
PRINTF("\n");
for (i=0; i<key->ncov; i++) {
keycov = &(key->cov[i]);
PRINTF("%d: meth=%s, cov=%d (%s) left=%d, param=", i,
METHODNAMES[keycov->method],
keycov->nr, CovList[keycov->nr].name,
keycov->left);
for (j=0; j<key->totalparam; j++) PRINTF("%f,", keycov->param[j]);
if (i < key->ncov - 1) PRINTF("\n%s\n",OP_SIGN[keycov->op]);
else PRINTF("\n\n");
}
for (i=0; i<key->ncov; i++) {
meth = &(key->meth[i]);
PRINTF("%d: ^S=%d ^destruct=%d\n unimeth=%d (%s) [%d]\n", i,
meth->S, meth->destruct,
meth->unimeth, METHODNAMES[meth->unimeth],
meth->unimeth_alreadyInUse
);
}
//
for (i=0; i<MAXDIM; i++)
PRINTF("%d : len=%d \n", i, key->length[i]);
PRINTF("^x=%d ^y=%d ^z=%d \n", key->x[0], key->x[1], key->x[2]);
}
void printKEY(int *keyNr) {
key_type *key;
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
PRINTF("\nkeynr out of range!\n");
return;
}
key = &(KEY[*keyNr]);
PRINTF("\nNr=%d", *keyNr);
printkey(key);
}
void printAllKeys() {
int i;
for (i=0;i<MAXKEYS;i++) {printKEY(&i);PRINTF("\n");}
}
int Transform2NoGrid(key_type *key, int v) {
int d;
covinfo_type *kc;
kc = &(key->cov[v]);
if (kc->simugrid) {
for (d=0; d<key->timespacedim; d++) {
kc->length[d] = kc->genuine_dim[d] ? key->length[d] : 1;
// printf("%d %d %d %d \n", d,
// kc->length[d], kc->genuine_dim[d], key->length[d]);
}
Getxsimugr(key->x, kc->param, key->timespacedim, key->anisotropy,
kc->xsimugr);
}
return Transform2NoGrid(key, kc->aniso, kc->reduceddim,
kc->simugrid, &(kc->x));
}
int Transform2NoGrid(key_type *key, aniso_type aniso, int reduceddim,
bool simugrid, double **xx)
{
/*
this function transforms the coordinates according to the anisotropy
matrix given in param (usually param is the param the first anisotropy
matrix of the covariance product -- all the other anisotropy matrices
have to multiple of the first, and this parameter is put into the
covariance function as scale parameter
developtime only considered if !key->simugrid && key->Time
output : keycov->xx
isotropic : grid : triple devided by scale ==multip. with aniso[0]!!
arbitrary : x-coordinates devided by scale
anisotropic : if grid and aniso=diag then
triples are multiplied with orig param[ANISO...]
else
x-coordinates multiplied explicitely with aniso: x*a;
dim-reduced
*/
int i,j,k,n,d,w,dimM1;
long endtime, total;
double t, *x;
// printf("%d %d %d %d %f %f %f %d\n",
// key->timespacedim,key->anisotropy,key->Time,key->grid,
// key->T[0],key->T[1],key->T[2],key->length[0]
// );
dimM1 = key->timespacedim - 1;
// assert(*xx == NULL);
if (*xx != NULL) {
if (GENERAL_PRINTLEVEL>-4)
PRINTF("Note: x values had already been transformed.");
return NOERROR;
}
// total number of coordinates to store
total = simugrid
? (key->timespacedim * 3)
: (reduceddim * key->totalpoints);
// printf("2nogr %d %d %d %d \n", total, reduceddim, simugrid, key->totalpoints);
// printf("total %d %d %d %d\n", total, key->totalpoints, reduceddim, key->timespacedim);
if ((x = *xx = (double*) malloc(sizeof(double) * total)) == NULL)
return ERRORMEMORYALLOCATION;
if (key->anisotropy) {
/* determine points explicitly */
if (key->Time && !key->grid) { // time component, no grid
endtime = key->length[key->spatialdim];
for (k=j=0, t=key->T[XSTART]; j<endtime; j++, t += key->T[XSTEP]) {
for (i=0; i<key->length[0]; i++) {
for (n=d=0; d<reduceddim; d++, k++) {
x[k] = 0.0;
for(w=0; w<key->spatialdim; w++) {
// printf("k=%d n=%d, w=%d, i=%d\n", k, n, w, i);
// printf("xk=%f aniso=%f %f\n", x[k], aniso[n], key->x[w][i]);
x[k] += aniso[n++] * key->x[w][i];
}
x[k] += aniso[n++] * t; //auf keinen Fall nach vorne holen
}
}
}
} else {
if (key->grid) {/* grid; with or without time component */
if (simugrid) { // necessarily stemming from diag matrix, but
// compressed according to zeros in the diagonal
for (i=0; i<3; i++) {
k = i;
for (n=d=0; d<key->timespacedim; d++, k+=3) {
// printf("2nogrid %d %f\n", k, x[k]);
x[k] = 0.0;
for(w=0; w<key->timespacedim; w++)
x[k] += aniso[n++] * key->x[w][i];
}
}
} else { // grid but not simugrid
double y[MAXDIM]; /* current point within grid, but without
anisotropy transformation */
int yi[MAXDIM]; /* counter for the current position in the grid */
for (w=0; w<key->timespacedim; w++) {y[w]=key->x[w][XSTART]; yi[w]=0;}
for (k=0; k<total; ){
// printf("%f %f %f %d %d %d\n", y[0], y[1], y[2], yi[0], yi[1], yi[2]);
for (n=d=0; d<reduceddim; d++, k++) {
x[k] = 0.0;
for(w=0; w<key->timespacedim; w++) {
x[k] += aniso[n++] * y[w];
}
}
i = 0;
(yi[i])++;
y[i] += key->x[i][XSTEP];
while(yi[i]>=key->length[i]) {
yi[i]=0;
y[i] = key->x[i][XSTART];
if (i<dimM1) {
i++;
(yi[i])++;
y[i] += key->x[i][XSTEP];
} else {
assert(k==total);
}
}
}
}
} else { /* anisoptropic, no grid, no time component */
for (k=i=0; i<key->totalpoints; i++)
for (n=d=0; d<reduceddim; d++, k++) {
x[k] = 0.0;
for(w=0; w<key->timespacedim; w++) x[k] += aniso[n++] * key->x[w][i];
}
}
}
} else { /* not key->anisotropy, no time component */
// printf("transform2 %d %d %d %d\n", key->anisotropy, key->grid,
// reduceddim, key->timespacedim);
assert(reduceddim == key->timespacedim);
if (key->grid) {
assert(simugrid);
for (k=d=0; d<key->timespacedim; d++) {
for (i=0; i<3; i++) {
x[k++] = key->x[d][i] * aniso[0];
}
}
} else {
for (k=i=0; i<key->totalpoints; i++) {
for (j=0; j<key->timespacedim; j++) x[k++] = key->x[j][i] * aniso[0];
}
}
}
return NOERROR;
}
int InternalGetGridSize(double *x[MAXDIM], int *dim, int *lx)
{
int d;
for (d=0; d<*dim; d++) {
if ((x[d][XEND]<=x[d][XSTART]) || (x[d][XSTEP]<=0))
return ERRORCOORDINATES;
lx[d] = 1 + (int)((x[d][XEND]-x[d][XSTART])/x[d][XSTEP]);
}
return NOERROR;
}
// obsolete?!
void GetGridSize(double *x,double *y,double *z,int *dim,int *lx,int *ly,int *lz)
{
// should now be consistent with `seq' of R
// if not, please report!!!!!!!!!!!!
double *xx[MAXDIM];
int lxx[MAXDIM];
xx[0]=x; xx[1]=y; xx[2]=z;
if (InternalGetGridSize(xx,dim,lxx)) {
*lx = *ly = *lz = 0;
} else {
*lx = lxx[0]; *ly = lxx[1]; *lz=lxx[2];
}
}
double SndDerivCovFct(double *x, int dim, covinfo_arraytype keycov,
covlist_type covlist, int ncov) {
// only for isotropic models and isotropy
// considers variance and the scale parameter
double fct[MAXCOV], abl[MAXCOV], snd[MAXCOV], zw, z, result = 0.0;
int v = 0, vold, w1, w2, i;
covinfo_type *kc;
cov_fct *cov=NULL; /* no meaning -- just avoids warning from -Wall */
assert(dim==1);
while (v<ncov) {
vold = v;
v++; while ((v<ncov) && keycov[covlist[v-1]].op) v++;
for (w1=vold; w1<v; w1++) {
kc = &(keycov[covlist[w1]]);
cov = &(CovList[kc->nr]);
z = fabs(x[0] * kc->aniso[0]);
fct[w1] = kc->param[VARIANCE] * cov->cov(&z, kc->param);
abl[w1] = kc->param[VARIANCE] * cov->derivative(&z, kc->param)
* kc->aniso[0];
snd[w1] = kc->param[VARIANCE] * cov->secondderivt(&z, kc->param)
* kc->aniso[0] * kc->aniso[0];
// printf("2nd: %d %f %f %f %f %s %d\n", w1,fct[w1], abl[w1], snd[w1],
// cov->secondderivt(&z, kc->param), cov->name, cov->secondderivt);
}
for (w1=vold; w1<v; w1++) {
for (w2=w1; w2<v; w2++) {
zw = (w1==w2 ? snd[w1] : (2.0 * abl[w1] * abl[w2]));
for (i=vold; i<v; i++)
if (i!=w1 && i!=w2) zw *= fct[i];
result += zw;
}
}
}
return result;
}
double DerivCovFct(double *x, int dim, covinfo_arraytype keycov,
covlist_type covlist, int ncov) {
// only for isotropic models and isotropy
// considers variance and the scale parameter
double fct[MAXCOV], abl[MAXCOV], zw, z, result = 0.0;
int v = 0, vold, i, w;
covinfo_type *kc;
cov_fct *cov=NULL; /* no meaning -- just avoids warning from -Wall */
assert(dim==1);
while (v<ncov) {
vold = v;
v++; while ((v<ncov) && keycov[covlist[v-1]].op) v++;
for (w=vold; w<v; w++) {
kc = &(keycov[covlist[w]]);
cov = &(CovList[kc->nr]);
z = fabs(x[0] * kc->aniso[0]);
// printf("phi1-C: %f %f %f %f\n", kc->param[VARIANCE], kc->aniso[0], x[0], z);
fct[w] = kc->param[VARIANCE] * cov->cov(&z, kc->param);
abl[w] = kc->param[VARIANCE] * cov->derivative(&z, kc->param)
* kc->aniso[0];
}
for (w=vold; w<v; w++) {
zw = abl[w];
for (i=vold; i<v; i++) if (i!=w) zw *= fct[i];
result += zw;
}
}
return result;
}
//int zaehler = 0;
double CovFct(double *x, int dim, covinfo_arraytype keycov,
covlist_type covlist, int ncov, bool anisotropy) {
int v, ani, k, j, endfor, ncovM1, nsub, subdim;
double var, zz, zw, result, d, z[MAXDIM + 2];
static double ZERO[MAXDIM] = {0,0,0,0};
covinfo_type *kc;
cov_fct *cov=NULL; /* no meaning -- just avoids warning from -Wall */
z[0] = z[1] = RF_NAN;
zw = result = 0.0;
ncovM1 = ncov - 1;
// printf("ani %d \n", anisotropy);
if (anisotropy) {
int cov_type;
for (v=0; v<ncov; v++) {
zw = var = 1.0;
for(; v<ncov; v++) {
kc = &(keycov[covlist[v]]);
cov = &(CovList[kc->nr]);
var *= kc->param[VARIANCE];
cov_type=cov->type;
for (k=0; k<kc->reduceddim; k++) z[k]=0.0;
for (ani=0, k=0; k<kc->reduceddim; k++) {
for (j=0; j<dim; j++) {
z[k] += kc->aniso[ani++] * x[j];
}
}
// printf("%d %f %f\n", dim, z[0], z[1]);
if (cov_type==ANISOTROPIC)
zw *= cov->cov(z, kc->param);
/* derzeit ist dim in cov->FCT bis auf assert-checks unbenutzt */
else {
/* the following definition allows for calculating the covariance */
/* function for spatialdim=0 and SPACEISOTROPIC model, namely as a */
/* purely temporal model; however, in CheckAndBuildCov such kind of */
/* calculations are prevend as being TRIVIAL(generally, independent*/
/* of being SPACEISOTROPIC or not) -- unclear whether relaxation in*/
/* CheckAndBuildCov possible */
if (cov_type != ANISOHYPERMODEL) {
//ISOTROPIC, ISOHYPERMODEL, spaceisotropic
zz = 0.0;
// printf("red=%d %f %f %f %f %f %f -- ",
// kc->reduceddim, z[0], z[1], z[2],
// kc->aniso[6], kc->aniso[7], kc->aniso[8]);
endfor = kc->reduceddim - 1;
for (j=0; j<endfor; j++) zz += z[j] * z[j];
if (cov_type==ISOTROPIC || cov_type==ISOHYPERMODEL)
zz += z[endfor] * z[endfor];
else z[1]=z[endfor]; // spaceisotropic
z[0] = sqrt(zz);
}
if (cov_type==ISOHYPERMODEL || cov_type==ANISOHYPERMODEL) {
nsub = (int) kc->param[HYPERNR];
subdim = (int) kc->param[EFFECTIVEDIM];
// subdim + 1 : C(0)
// subdim : C(x)
// printf("A %d %f\n", subdim, z[0]);
z[subdim] =CovFct(x, dim, keycov, &(covlist[v+1]), nsub, anisotropy);
// printf("B\n");
z[subdim+1] = CovFct(ZERO, dim, keycov, &(covlist[v+1]), nsub,
anisotropy);
// printf("%f %f %f %f %d \n",zw, z[0], z[1], z[2], subdim);
// printf("zw=%f ", zw);
zw *= cov->cov(z, kc->param);
// printf("%f %f %f var=%f %f %f %f subd=%d %s %f %f\n",
// z[0], z[1], z[2], kc->param[VARIANCE],
// kc->param[KAPPA1], kc->param[KAPPA2], kc->param[KAPPA3],
// subdim, cov->name, cov->cov(z, kc->param, subdim), zw);
//
// assert(zw < 100000000.00); // nicht loeschen
v += nsub;
kc = &(keycov[covlist[v]]);
} else {
zw *= cov->cov(z, kc->param);
}
}
if ((v<ncovM1) && (kc->op==0)) break;
}
result += var * zw; //
}
} else {
if (dim==1) d=fabs(x[0]);
else {
for (d=0.0, j=0; j<dim; j++) {d += x[j] * x[j];}
d = sqrt(d);
}
for (v=0; v<ncov; v++) {
zw = var = 1.0;
for(; v<ncov; v++) {
//printf("covfct %d %d\n", v, covlist[v]);
kc = &(keycov[covlist[v]]);
// printf("covfct %d %d %s %f\n", v, kc->nr, CovList[kc->nr].name,
// kc->aniso[0]);
cov = &(CovList[kc->nr]); /* cov needed in FORMULA for Vario */
if (cov->type==ISOHYPERMODEL) {
nsub = (int) kc->param[HYPERNR];
z[0] = d * kc->aniso[0];
// subdim + 1 : C(0)
// subdim : C(x)
z[1] = CovFct(x, dim, keycov, &(covlist[v+1]), nsub,
anisotropy);
z[2] = CovFct(ZERO, dim, keycov, &(covlist[v+1]), nsub,
anisotropy);
// if (zaehler < 100) printf("----- %f %f %f %f\n", z[0], z[1],
// d, kc->aniso[0]);
zw *= cov->cov(z, kc->param);
v += nsub;
kc = &(keycov[covlist[v]]);
} else {
zz = d * kc->aniso[0];
var *= kc->param[VARIANCE];
zw *= cov->cov(&zz, kc->param);
}
if ((v<ncovM1) && (kc->op==0)) break;
}
result += var * zw; //
// if ( zaehler++ < 20)
// printf("-->> %f %f %f %f %f %d %d %d\n", result, var, zw, zz,
// kc->param[KAPPA], v, covlist[v], kc->nr);
}
}
// if ( zaehler++ < 100)
// printf(" %f (%f, %f) %f d=%f %f %d %d %d -- zz=%f, var=%f zw=%f\n",
// *x, z[0], z[1], result, d, kc->aniso[0],
// ncov, keycov[0].nr, keycov[ncov==1 ? 0 : 1].nr
// ,zz, var, zw
// );
// if (dim==2)
// printf("X%d %f %f, %f %f: %f %f \n",dim, x[0], x[1], z[0], z[1], zw, result);
// else
// printf("-%d %f, %f : %f %f \n",dim, x[0], z[0], zw, result);
// assert(result > 0.4);
return result;
}
char STR[30];
char *strround(double x){
if (x == floor(x + 0.5) ) sprintf(STR, "%d", (int) x);
else sprintf(STR, "%f", x);
return STR;
}
void addleft(char *str, char *sign, double x) {
sprintf(str, "%s%s%s", str, strround(x), sign);
}
void addright(char *str, char *sign, double x) {
sprintf(str, "%s%s%s", str, sign, strround(x));
}
void addone(char *str, double x) {
sprintf(str, "%s, %s", str, strround(x));
}
void addpair(char *str, double x, double y) {
if (x == floor(x + 0.5) && y==floor(y + 0.5))
sprintf(str, "%s, (%d,%d)", str, (int) x, int(y));
else
sprintf(str, "%s, (%f,%f)", str, x, y);
}
int check_within_range(param_type param, cov_fct *cov, int timespacedim,
char *txt) {
int error, index, i4, i, kappas, maxdim;
double range[LASTKAPPA * 4];
bool truesegm;
rangefct getrange;
if (GENERAL_SKIPCHECKS) return NOERROR;
getrange = cov->range;
kappas = cov->kappas(timespacedim);
index = 1;
for(;;) {
getrange(timespacedim, &index, range);
if (index==RANGE_INVALIDDIM) {
maxdim = timespacedim;
while (index==RANGE_INVALIDDIM) { // check when dim starts being valid
maxdim--;
assert(maxdim > 0);
index = 1;
getrange(maxdim, &index, range); // range not of interest anymore
}
sprintf(ERRORSTRING_WRONG, "genuine dimension = %d%s",
timespacedim, txt == NULL ? "" : txt);
sprintf(ERRORSTRING_OK, "dim <= %d", maxdim);
return ERRORCOVFAILED;
}
error = -1;
truesegm = true;
for (i=0; i<kappas; i++) {
i4 = i * 4;
if (range[i4] == floor(range[i4]) + OPEN) {
if (param[KAPPA1 + i] <= range[i4] - OPEN) error=i;
} else {
if (param[KAPPA1 + i] < range[i4]) error=i;
}
i4++;
if (range[i4] + OPEN == ceil(range[i4])) {
if (param[KAPPA1 + i] >= range[i4] + OPEN) error=i;
} else {
if (param[KAPPA1 + i] > range[i4]) error=i;
}
truesegm &= range[i4] != range[i4-1] || param[KAPPA1 + i] == range[i4];
/*
// printf("%d %d %d %d %d %f %f %f %d\n", i,
param[KAPPA1 + i] < range[--i4],
param[KAPPA1 + i] == range[i4] &&range[i4] == floor(range[i4]) + OPEN ,
param[KAPPA1 + i] > range[++i4] ,
param[KAPPA1 + i] == range[i4] && range[i4] + OPEN == ceil(range[i4]),
range[i4-1], range[i4], param[KAPPA1 + i], truesegm
);
*/
}
if (truesegm) { // exakte indizierung darf nur 1x auftauchen
// kleine Einschr\"ankung, aber bisher kein Problem
if (error>=0) {
i = error;
i4 = i * 4;
sprintf(ERRORSTRING_WRONG, "kappa%d=%f", i+1, param[KAPPA1 + i]);
strcpy(ERRORSTRING_OK, "");
if (range[i4] > RF_NEGINF) {
if (range[i4] == floor(range[i4]) + OPEN)
addleft(ERRORSTRING_OK, "<", floor(range[i4]));
else
addleft(ERRORSTRING_OK, "<=", range[i4]);
}
i4++;
sprintf(ERRORSTRING_OK, "%skappa%d", ERRORSTRING_OK, i+1);
if (range[i4] < RF_INF) {
printf("%f %f %d %e\n", range[i4] + OPEN, ceil(range[i4]),
range[i4] + OPEN == ceil(range[i4]),
range[i4] + OPEN - ceil(range[i4]));
if (range[i4] + OPEN == ceil(range[i4]))
addright(ERRORSTRING_OK, "<", ceil(range[i4]));
else
addright(ERRORSTRING_OK, "<=", range[i4]);
}
sprintf(ERRORSTRING_OK, "%s when %s dim=%d%s",
ERRORSTRING_OK,
(cov->type == SPACEISOTROPIC) ? "genuine spatial" :
(cov->type == ISOTROPIC || cov->type==ISOHYPERMODEL)
? "genuine" : "",
timespacedim - (int) (cov->type == SPACEISOTROPIC),
txt == NULL ? "" : txt);
// printf("%e %e\n", range[i4-1], range[i4]);
// assert(false);
return ERRORCOVFAILED;
} else return NOERROR;
} else {
if (index==RANGE_LASTELEMENT) {
int n, idx[LASTKAPPA];
n = 0;
for (i4=i=0; i<kappas; i++, i4+=4)
if (range[i4] == range[i4+1]) {
idx[n] = i;
n++;
}
index = 1;
assert(n>0);
if (n==1) {
sprintf(ERRORSTRING_OK, "kappa%d =", idx[0]+1);
} else {
sprintf(ERRORSTRING_OK,"(kappa%d, kappa%d) =", idx[0]+1, idx[1]+1);
}
while (index > 0) {
cov->range(timespacedim, &index, range);
if (n==1) {
addone(ERRORSTRING_OK, range[idx[0] * 4]);
} else {
addpair(ERRORSTRING_OK, range[idx[0] * 4], range[idx[1] *4]);
}
}
if (n==1) {
sprintf(ERRORSTRING_WRONG, "kappa%d =", idx[0]+1);
addone(ERRORSTRING_WRONG, param[KAPPA1]);
} else {
sprintf(ERRORSTRING_WRONG,"(kappa%d, kappa%d) =",
idx[0]+1, idx[1]+1);
addpair(ERRORSTRING_WRONG,
param[KAPPA1+idx[0]], param[KAPPA1+idx[1]]);
}
return ERRORCOVFAILED;
}
}
}
}
// checks a *single* basic covariance function
// p and param must be of the same kind (double)!!!
int CheckAndBuildCov(int *covnr, int *op, int ncov,
double *ParamList, int nParam,
int naturalscaling, bool Time, bool grid,
int timespacedim, bool anisotropy,
covinfo_arraytype keycov, bool *equal)
{
// IMPORTANT!: p is the parameter vector of the R interface !!
// output param : equals p, except naturalscaling
int error, v, i, totkappas, pAniso, d;
double newscale;
covinfo_type *kc;
cov_fct *cov;
param_type param;
if (currentNrCov==-1) InitModelList(); assert(CovList!=NULL);
if (ncov<0 || ncov>MAXCOV) return ERRORNCOVOUTOFRANGE;
if (timespacedim < 1 || timespacedim > MAXDIM) return ERRORDIM;
pAniso = anisotropy ? timespacedim * timespacedim : 1;
if (Time && !anisotropy) return ERRORTIMENOTANISO;
totkappas = 0;
for (v=0; v<ncov; v++){
ERRORMODELNUMBER = v;
if (covnr[v]<0 || covnr[v]>=currentNrCov) return ERRORCOVNROUTOFRANGE;
*equal &= keycov[v].nr == covnr[v];
keycov[v].nr = covnr[v];
kc = &(keycov[v]);
kc->dim = timespacedim;
cov = &(CovList[kc->nr]);
if (GENERAL_PRINTLEVEL > 7)
PRINTF("Check&Build %d %d covnr=%d op=%d\n", v, ncov, covnr[v],
(v<ncov-1) ? op[v] : -1);
// Variance and Kappa
int kappas, kappasP1;
kappas = cov->kappas(timespacedim);
kappasP1 = kappas + 1;
*equal &= !memcmp(kc->param, ParamList, sizeof(double) * kappasP1);
memcpy(kc->param, ParamList, sizeof(double)* kappasP1);
ParamList += kappasP1;
totkappas += kappas;
if (kc->param[VARIANCE]<0.0) return ERRORNEGATIVEVAR;
switch(cov->type) {
case ISOTROPIC : case ISOHYPERMODEL :
kc->param[EFFECTIVEDIM] = 1.0;
break;
case SPACEISOTROPIC :
kc->param[EFFECTIVEDIM] = 2.0;
if (!Time) return ERRORWITHOUTTIME;
break;
case ANISOTROPIC : case ANISOHYPERMODEL :
kc->param[EFFECTIVEDIM] = (double) timespacedim;
break;
}
// op
if (v < ncov - 1) {
*equal &= kc->op == op[v];
kc->op = op[v];
// printf(" type %d %d %d %d %d %d %d\n", v, cov->type, ISOHYPERMODEL,
// kc->op, false xor false, false xor true, true xor true);
if ((cov->type == ISOHYPERMODEL || cov->type == ANISOHYPERMODEL)
xor (kc->op == 2))
return ERRORPARENTHESIS;
}
// not hypermodel: ?
// natural scaling
GetNaturalScaling(&(kc->nr), &(kc->param[KAPPA]),
&naturalscaling, &newscale, &error);
// printf("&naturalscaling = %d \n", naturalscaling);
if (error != NOERROR) return error;
if (anisotropy) newscale = 1.0 / newscale;
for (i=pAniso - 1; i>=0; i--)
param[ANISO + i] = newscale * ParamList[i];
if (!anisotropy) {
if (cov->type==SPACEISOTROPIC || cov->type==ANISOTROPIC ||
cov->type==ANISOHYPERMODEL)
return ERRORNOTANISO;
if (cov->cov==nugget && param[SCALE]==0.0) param[SCALE] = 1.0;
if (param[SCALE] <= 0.0) return ERRORNEGATIVESCALE;
}
ParamList += pAniso;
*equal &= !memcmp(&(kc->param[ANISO]), &(param[ANISO]),
sizeof(double)* pAniso);
memcpy(&(kc->param[ANISO]), &(param[ANISO]), sizeof(double) * pAniso);
// reduced anisotropy matrix and true grid, time and dimension
///////
// for hypermodels it is very unclear, how the dimension should be taken
GetTrueDim(anisotropy, timespacedim, kc->param, cov->type,
&(kc->genuine_last_dimension),
&(kc->reduceddim), kc->aniso);
if (kc->reduceddim==0 || // e.g. used in RFnugget and Spectral
(cov->type==SPACEISOTROPIC &&
(!kc->genuine_last_dimension || kc->reduceddim<2))
) return ERRORLOWRANK;
// check value of parameters
error = check_within_range(kc->param, cov, kc->reduceddim, NULL);
if (error != NOERROR) return error;
int maxdim, dummy;
cov->info(kc->param, &maxdim, &dummy);
// printf("%d %d\n", maxdim, kc->reduceddim);
if (maxdim < kc->reduceddim && !GENERAL_SKIPCHECKS)
return ERRORCOVNOTALLOWED;
kc->simugrid = // kc->param[ANISO] !!
grid && (!anisotropy || is_diag(&(kc->param[ANISO]), timespacedim));
if (kc->simugrid) {
if (anisotropy) {
int n;
for (i=n=d=0; d<timespacedim; d++, n += timespacedim + 1) {
if ((kc->genuine_dim[d] = param[ANISO + n] != 0.0)) {
// *** funktioniert nur deshalb weil nachgeprueft
// worden ist das diag matrix ***
kc->idx[d] = i++;
// WICHTIG! Dadurch wird Ordnung
// der Dimensionen eingehalten!!
// siehe auch Transform2nogrid
} else {
kc->idx[d] = kc->reduceddim + d - i;
}
}
} else {
for (d=0; d<timespacedim; d++) {
kc->genuine_dim[d] = true;
kc->idx[d] = d;
}
}
}
if (cov->type == ISOHYPERMODEL || cov->type == ANISOHYPERMODEL) {
////////////////////////
kc->simugrid = false; // not programmed yet -- save version !!!!!!
///////////////////////
} else {
error = cov->check(kc->param, kc->reduceddim, kc->dim, Nothing);
if (error != NOERROR) return error;
}
}
ERRORMODELNUMBER = -1;
for (v=ncov-1; v>=0; v--) {
int w;
ERRORMODELNUMBER = v;
// this loop cannot be intregrated above. The subsequent
// submodels must be initialised first, since
// 1) checkNinit relies on the knowledge of the submodels
// 2) parameters are taken over from the first submodel
// this loop must be in reverse order for the same reasons.
kc = &(keycov[v]);
cov = &(CovList[kc->nr]);
if (cov->type == ISOHYPERMODEL || cov->type == ANISOHYPERMODEL) {
// take over the anisotropy structure from the first submodel
// 28.9.05: not a good idea !!!
// memcpy(&(kc->param[ANISO]), &(keycov[v+1].param[ANISO]),
// sizeof(double) * pAniso);
/*
memcpy(&(kc->param[ANISO]), &(keycov[v+1].param[ANISO]),
sizeof(double) * pAniso);
GetTrueDim(anisotropy, timespacedim, kc->param,
CovList[keycov[v+1].nr].type,
&(kc->genuine_last_dimension), //
&(kc->reduceddim), kc->aniso);//
kc->aniso[0] = RF_NAN;
*/
for (w = v + (int) (kc->param[HYPERNR]); w>v; w--) {
if (CovList[keycov[w].nr].type != ISOTROPIC) {
// z.B. tbm viel complizierter
error = ERRORHYPERNOTISO;
}
}
error = cov->checkNinit(kc, &(COVLISTALL[v]), ncov-v-1, Nothing,
anisotropy);
if (error != NOERROR) {
return error;
}
}
}
ERRORMODELNUMBER = -1;
if (ncov * (1 + pAniso) + totkappas != nParam) return ERRORPARAMNUMBER;
return NOERROR;
}
static int user_dim=-1, user_ncov=-1;
static bool user_anisotropy=false;
static covinfo_arraytype user_keycov;
void CheckAndCompleteParameters(int *covnr, double *p, int *np,
int *dim, /* timespacedim ! */
int *ncov, int *anisotropy, int *op,
double *q, int* error) {
int v;
bool dummyequal = false;
// note: column of x is point !!
COVINFO_NULL(user_keycov);
user_dim = *dim;
user_ncov=*ncov;
user_anisotropy=(bool) *anisotropy;
*error = CheckAndBuildCov(covnr, op, user_ncov, p, *np, GENERAL_NATURALSCALING,
user_anisotropy /* Time */, false /* grid*/,
user_dim, user_anisotropy, user_keycov,
&dummyequal);
if (*error != NOERROR) {
if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing, *error);
return;
}
for (v=0; v<user_ncov; v++, q += TOTAL_PARAM) {
memcpy(q, user_keycov[v].param, sizeof(param_type));
}
}
int CheckCovariance(int *covnr, double *p, int np,
int logicaldim, /* timespacedim ! */
int xdim, int ncov, bool anisotropy, int *op,
int naturalscaling, covinfo_arraytype keycov)
{
int v, error;
bool dummyequal = false;
cov_fct *cov;
if ((logicaldim!=xdim) && (anisotropy || (xdim!=1))) return ERRORDIMMISMATCH;
error = CheckAndBuildCov(covnr, op, ncov, p, np, naturalscaling,
anisotropy /* Time */, false /* grid*/,
logicaldim, anisotropy, keycov, &dummyequal);
if (error != NOERROR) return error;
for (v=0; v<ncov; v++) {
ERRORMODELNUMBER = v;
covinfo_type *kc;
kc = &(keycov[v]);
cov = &(CovList[kc->nr]);
if (cov->cov==NULL) return ERRORNOTPROGRAMMED;
if (cov->type==ISOHYPERMODEL || cov->type==ANISOHYPERMODEL) {
v += (int) kc->param[HYPERNR];
} else {
if (cov->variogram) return ERRORNOTDEFINED;
}
}
ERRORMODELNUMBER = -1;
return NOERROR;
}
void CalculateCovariance(double *x, int lx, covinfo_arraytype keycov, int xdim,
int ncov, bool anisotropy, double *result) {
int i;
for (i=0; i<lx; i++, x += xdim)
result[i] = CovFct(x, xdim, keycov, COVLISTALL, ncov, anisotropy);
}
// note: number of parameters is not checked whether it is correct!!
// done by R function `PrepareModel'
void CovarianceNatSc(double *x, int *lx, int *covnr, double *p, int *np,
int *logicaldim, /* timespacedim ! */
int *xdim,
int *ncov, int *anisotropy, int *op,
double *result, int *naturalscaling)
{
int error, i;
// note: column of x is point !!
COVINFO_NULL(user_keycov);
user_dim = *logicaldim;
user_ncov=*ncov;
user_anisotropy=(bool) *anisotropy;
if ((error = CheckCovariance(covnr, p, *np, user_dim, *xdim, user_ncov,
user_anisotropy,
op, *naturalscaling, user_keycov)) != NOERROR) {
if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing, error);
for (i=0; i<*lx; i++) result[i] = RF_NAN;
return;
}
// die nachfolgende Verschraenkung von logicaldim und xdim
// funktioniert, siehe code von CheckAndBuildCov und GetTrueDim
CalculateCovariance(x, *lx, user_keycov, *xdim, user_ncov, user_anisotropy,
result);
}
void Covariance(double *x,int *lx, int *covnr, double *p, int *np,
int *logicaldim, /* timespacedim ! */
int *xdim, int *ncov, int *anisotropy, int *op,
double *result)
{
CovarianceNatSc(x, lx, covnr, p, np, logicaldim, xdim,
ncov, anisotropy, op, result, &GENERAL_NATURALSCALING);
}
void CalculateCovMatrix(double *dist, int lx, covinfo_arraytype keycov,
int xdim, int ncov, bool anisotropy, double *result) {
// note: column of x is point !!
int i, ii, endfor, lxP1;
long ve, ho;
double var;
var = CovFct(ZERO, xdim, keycov, COVLISTALL, ncov, anisotropy);
lxP1 = lx + 1;
for (ii=lx, i=0; ii>0; i+=lxP1, ii--) {
result[i] = var;
endfor = i + ii;
for (ve = i + 1, ho = i + lx; ve < endfor; ve++, ho += lx, dist += xdim) {
result[ve] = result[ho] =
CovFct(dist, xdim, keycov, COVLISTALL, ncov, anisotropy);
}
}
}
void CovarianceMatrixNatSc(double *dist,int *lx, int *covnr,double *p,
int *np, int *logicaldim, /* timespacedim ! */
int *xdim, int *ncov, int *anisotropy, int *op,
double *result, int *naturalscaling)
{
// note: column of x is point !!
int error;
long i, lxq;
COVINFO_NULL(user_keycov);
user_dim = *logicaldim;
user_ncov=*ncov;
user_anisotropy=(bool) *anisotropy;
error = CheckCovariance(covnr, p, *np, user_dim, *xdim, user_ncov,
user_anisotropy, op, *naturalscaling, user_keycov);
if (error != NOERROR) {
if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing, error);
lxq = *lx * *lx;
for (i=0; i<lxq; i++) result[i] = RF_NAN;
return;
}
CalculateCovMatrix(dist, *lx, user_keycov, *xdim, user_ncov,
user_anisotropy, result);
}
void CovarianceMatrix(double *dist, int *lx,int *covnr, double *p, int *np,
int *logicaldim, /* timespacedim ! */
int *xdim,
int *ncov, int *anisotropy, int *op,
double *result)
{
CovarianceMatrixNatSc(dist, lx, covnr, p, np, logicaldim, xdim, ncov,
anisotropy, op, result, &GENERAL_NATURALSCALING);
}
static int Unchecked_dim=-1, Unchecked_ncov=-1;
static bool Unchecked_anisotropy=false;
static covinfo_arraytype Unchecked_keycov;
void InitUncheckedCovFct(int *covnr, double *p, int *np,
int *logicaldim, /* timespacedim ! */
int *xdim, int *ncov, int *anisotropy, int *op,
int *naturalscaling, int *error)
{
if (Unchecked_dim==-1) COVINFO_NULL(Unchecked_keycov);
if ((*error = CheckCovariance(covnr, p, *np, *logicaldim, *xdim, *ncov,
(bool) *anisotropy, op, *naturalscaling,
Unchecked_keycov)) != NOERROR) {
if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing, *error);
return;
}
Unchecked_ncov = *ncov;
Unchecked_anisotropy = (bool) *anisotropy;
Unchecked_dim = *xdim;
return;
}
void UncheckedCovFct(double *x, int *lx, double *result)
{
CalculateCovariance(x, *lx, Unchecked_keycov, Unchecked_dim, Unchecked_ncov,
Unchecked_anisotropy, result);
}
void UncheckedCovMatrix(double *dist, int *lx, double *result)
// corresponds to CovarianceMatrixNatSc
{
CalculateCovMatrix(dist, *lx, Unchecked_keycov, Unchecked_dim, Unchecked_ncov,
Unchecked_anisotropy, result);
}
int CheckVariogram(int *covnr, double *p, int np,
int logicaldim, /* timespacedim ! */
int xdim, int ncov, bool anisotropy, int *op,
int naturalscaling, covinfo_arraytype keycov)
{
int v, error;
bool dummyequal = false;
cov_fct *cov;
if ((logicaldim!=xdim) && (anisotropy || (xdim!=1))) return ERRORDIMMISMATCH;
// printf("Cv %d %d %d\n", logicaldim, xdim, anisotropy);
error = CheckAndBuildCov(covnr, op, ncov, p, np, naturalscaling,
anisotropy, /* Time */ false, /* grid*/
logicaldim, anisotropy, keycov, &dummyequal);
if (error != NOERROR) return error;
for (v=0; v<ncov; v++) {
ERRORMODELNUMBER = v;
cov = &(CovList[covnr[v]]);
if (cov->cov==NULL) return ERRORNOTPROGRAMMED;
// below is too restrictive since there could be a hyper model
// with a multiplicative model as submodel
if (v>0 && op[v-1] && (cov->variogram || CovList[covnr[v-1]].variogram))
return ERRORNOMULTIPLICATION;
}
ERRORMODELNUMBER = -1;
return NOERROR;
}
void CalculateVariogram(double *x, int lx, covinfo_arraytype keycov, int xdim,
int ncov, bool anisotropy, double *result) {
int i;
double C0;
C0 = CovFct(ZERO, xdim, keycov, COVLISTALL, ncov, anisotropy);
for (i=0; i<lx; i++, x += xdim)
result[i] = C0 - CovFct(x, xdim, keycov, COVLISTALL, ncov, anisotropy);
}
void VariogramNatSc(double *x,int *lx,int *covnr,double *p,int *np,
int *logicaldim, /* timespacedim ! */
int *xdim,
int *ncov, int *anisotropy, int *op,
double *result, int *naturalscaling)
{
int error, i;
// printf("Cv %d %d %d\n", *logicaldim, *xdim, *anisotropy);
COVINFO_NULL(user_keycov);
user_dim = *logicaldim;
user_ncov=*ncov;
user_anisotropy=(bool) *anisotropy;
if ((error = CheckVariogram(covnr, p, *np, user_dim, *xdim, user_ncov,
user_anisotropy,
op, *naturalscaling, user_keycov)) != NOERROR) {
if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing, error);
for (i=0; i<*lx; i++) result[i] = RF_NAN;
return;
}
CalculateVariogram(x, *lx, user_keycov, *xdim, user_ncov, user_anisotropy,
result);
}
void Variogram(double *x,int *lx,int *covnr,double *p,int *np,
int *logicaldim, /* timespacedim ! */
int *xdim,
int *ncov, int *anisotropy, int *op,
double *result)
{
VariogramNatSc(x, lx, covnr, p, np, logicaldim, xdim, ncov, anisotropy, op,
result, &GENERAL_NATURALSCALING);
}
void VariogramMatrix(double *dist, int *lx, int *covnr, double *p, int *np,
double *result)
{
// not programmed yet
assert(false);
}
/*
int insert_method(int* last_incompatible, key_type *key,
SimulationType m, bool incompatible) {
int idx;
assert(m != Nugget || !incompatible);
assert(key->meth[key->n_unimeth].S==NULL);
assert(key->meth[key->n_unimeth].destruct==NULL);
if (incompatible) {
// first the non-compatible (i.e. which do not allow direct
// adding to the random field), then those which allow for
// direct adding
// --if attention wouldn't be payed to the ordering of the methods,
// the simulation algorithm had to create an intermediate field
// more frequently -- sometimes the fields are pretty large,
// and then this procedure is quite helpful
//
// da waehrend der for-schleife unten eingefuegt wird:
// immer moeglichst weit hinten einfuegen!
(*last_incompatible)++;
key->meth[key->n_unimeth].unimeth_alreadyInUse =
key->meth[*last_incompatible].unimeth_alreadyInUse;
key->meth[key->n_unimeth].unimeth=key->meth[*last_incompatible].unimeth;
key->meth[key->n_unimeth].S=key->meth[*last_incompatible].S; // necessary
// only if further methods are tried due to partial failure
key->meth[key->n_unimeth].destruct=key->meth[*last_incompatible].destruct;
key->meth[key->n_unimeth].incompatible =
key->meth[*last_incompatible].incompatible;
key->meth[*last_incompatible].S = NULL;
key->meth[*last_incompatible].destruct = NULL;
idx = *last_incompatible;
} else idx=key->n_unimeth;
key->meth[idx].unimeth = m;
key->meth[idx].incompatible = incompatible;
key->meth[idx].unimeth_alreadyInUse = false;
key->n_unimeth++;
return idx;
}
void delete_method(int *M, int *last_incompatible, key_type *key) {
// note that key->n_unimeth had been decreased before this function
// is called
int idx;
if (key->meth[*M].destruct != NULL) key->meth[*M].destruct(&(key->meth[*M].S));
assert(key->meth[*M].S==NULL);
key->meth[*M].destruct=NULL;
if (key->meth[*M].incompatible) {
assert(*last_incompatible>=0);
assert(*M<=*last_incompatible);
key->meth[*M].unimeth = key->meth[*last_incompatible].unimeth;
key->meth[*M].unimeth_alreadyInUse =
key->meth[*last_incompatible].unimeth_alreadyInUse;
key->meth[*M].S = key->meth[*last_incompatible].S;
key->meth[*last_incompatible].S=NULL;
key->meth[*M].incompatible = key->meth[*last_incompatible].incompatible;
key->meth[*M].destruct=key->meth[*last_incompatible].destruct;
key->meth[*last_incompatible].destruct=NULL;
idx = *last_incompatible;
(*last_incompatible)--;
} else idx = *M;
key->n_unimeth--;
if (idx==key->n_unimeth){ // deleting of the very last entry
assert(*M==key->n_unimeth || *last_incompatible + 1 == key->n_unimeth);
} else {
assert(idx>=0 && idx<key->n_unimeth);
key->meth[idx].unimeth = key->meth[key->n_unimeth].unimeth;
key->meth[idx].unimeth_alreadyInUse =
key->meth[key->n_unimeth].unimeth_alreadyInUse;
key->meth[idx].S = key->meth[key->n_unimeth].S;
key->meth[key->n_unimeth].S=NULL;
key->meth[idx].incompatible = key->meth[key->n_unimeth].incompatible;
key->meth[idx].destruct = key->meth[key->n_unimeth].destruct;
key->meth[key->n_unimeth].destruct=NULL;
}
key->meth[key->n_unimeth].unimeth_alreadyInUse = false;
(*M)--;
}
*/
int insert_method(int* last_incompatible, key_type *key,
SimulationType m, bool incompatible) {
int idx;
assert(m != Nugget || !incompatible);
if (key->meth[key->n_unimeth].destruct != NULL)
key->meth[key->n_unimeth].destruct(&(key->meth[key->n_unimeth].S));
assert(key->meth[key->n_unimeth].S==NULL);
key->meth[key->n_unimeth].destruct=NULL;
if (incompatible) {
// first the non-compatible (i.e. which do not allow direct
// adding to the random field), then those which allow for
// direct adding
// --if attention wouldn't be payed to the ordering of the methods,
// the simulation algorithm had to create an intermediate field
// more frequently -- sometimes the fields are pretty large,
// and then this procedure is quite helpful
//
// da waehrend der for-schleife unten eingefuegt wird:
// immer moeglichst weit hinten einfuegen!
(*last_incompatible)++;
key->meth[key->n_unimeth].unimeth_alreadyInUse =
key->meth[*last_incompatible].unimeth_alreadyInUse;
key->meth[key->n_unimeth].unimeth=key->meth[*last_incompatible].unimeth;
key->meth[key->n_unimeth].S=key->meth[*last_incompatible].S; // necessary
// only if further methods are tried due to partial failure
key->meth[key->n_unimeth].destruct=key->meth[*last_incompatible].destruct;
key->meth[key->n_unimeth].incompatible =
key->meth[*last_incompatible].incompatible;
key->meth[*last_incompatible].S = NULL;
key->meth[*last_incompatible].destruct = NULL;
idx = *last_incompatible;
} else idx=key->n_unimeth;
key->meth[idx].unimeth = m;
key->meth[idx].incompatible = incompatible;
key->meth[idx].unimeth_alreadyInUse = false;
key->n_unimeth++;
return idx;
}
void delete_method(int *M, int *last_incompatible, key_type *key) {
// note that key->n_unimeth had been decreased before this function
// is called
int idx, i;
methodvalue_type dummy;
dummy.S = key->meth[*M].S;
dummy.destruct= key->meth[*M].destruct;
dummy.unimeth = key->meth[*M].unimeth;
dummy.incompatible = key->meth[*M].incompatible;
if (key->meth[*M].incompatible) {
assert(*last_incompatible>=0);
assert(*M<=*last_incompatible);
key->meth[*M].unimeth = key->meth[*last_incompatible].unimeth;
key->meth[*M].unimeth_alreadyInUse =
key->meth[*last_incompatible].unimeth_alreadyInUse;
key->meth[*M].S = key->meth[*last_incompatible].S;
key->meth[*last_incompatible].S=NULL;
key->meth[*M].incompatible = key->meth[*last_incompatible].incompatible;
key->meth[*M].destruct=key->meth[*last_incompatible].destruct;
key->meth[*last_incompatible].destruct=NULL;
idx = *last_incompatible;
(*last_incompatible)--;
} else idx = *M;
key->n_unimeth--;
if (idx==key->n_unimeth){ // deleting of the very last entry
assert(*M==key->n_unimeth || *last_incompatible + 1 == key->n_unimeth);
} else {
assert(idx>=0 && idx<key->n_unimeth);
key->meth[idx].unimeth = key->meth[key->n_unimeth].unimeth;
key->meth[idx].unimeth_alreadyInUse =
key->meth[key->n_unimeth].unimeth_alreadyInUse;
key->meth[idx].S = key->meth[key->n_unimeth].S;
key->meth[key->n_unimeth].S=NULL;
key->meth[idx].incompatible = key->meth[key->n_unimeth].incompatible;
key->meth[idx].destruct = key->meth[key->n_unimeth].destruct;
key->meth[key->n_unimeth].destruct=NULL;
}
key->meth[key->n_unimeth].unimeth_alreadyInUse = false;
key->meth[key->n_unimeth].S = dummy.S;
key->meth[key->n_unimeth].destruct = dummy.destruct;
key->meth[key->n_unimeth].unimeth = dummy.unimeth;
key->meth[key->n_unimeth].incompatible = dummy.incompatible;
key->meth[key->n_unimeth].actcov = key->meth[*M].actcov;
for (i=key->meth[*M].actcov-1; i>=0; i--)
key->meth[key->n_unimeth].covlist[i] = key->meth[*M].covlist[i];
(*M)--;
}
/*
speed, exactness, stationarity
direct (lower -- always, upper -- never)
grid: circembed
nongrid: spectralTBM, TBM2, TBM3,
Ausnahmen: dimension beschraenkt
*/
void init_method_list(key_type *key)
{
unsigned int maxmem=500000000;
int best_dirct=DIRECTGAUSS_BESTVARIABLES,
max_dirct=DIRECTGAUSS_MAXVARIABLES,
i, v, nr, nml;
bool anyFirstDirect = false, anyFirstCE = false;
#define nLastDefault 2
SimulationType LastDefault[nLastDefault] = {AdditiveMpp, Hyperplane};
SimulationType DirectFst = key->totalpoints <= best_dirct ? Direct : Forbidden;
SimulationType DirectLst = key->totalpoints <= max_dirct ? Direct : Forbidden;
#define nStandard 6
#define TBM3pos 5
#define TBM2pos 4
SimulationType Standard[nStandard] =
{CircEmbed, CircEmbedIntrinsic, CircEmbedCutoff, SpectralTBM, TBM2, TBM3};
bool Allowed[SimulationTypeLength], allowed[SimulationTypeLength];
for (i=0; i<Nothing; i++) Allowed[i] = true;
if (DECISION_PARAM.stationary_only==DECISION_TRUE)
Allowed[CircEmbedIntrinsic]=false;
if (!key->grid) Allowed[CircEmbed] = Allowed[CircEmbedCutoff] =
Allowed[CircEmbedIntrinsic] = false;
if (DECISION_PARAM.exactness==DECISION_TRUE)
Allowed[TBM2] = Allowed[TBM3] = Allowed[SpectralTBM] =
Allowed[AdditiveMpp] = Allowed[Hyperplane] = false;
#define nraw 2 * SimulationTypeLength
SimulationType raw[nraw];
for (v=0; v<key->ncov; v++) {
cov_fct *cov;
cov = &(CovList[key->cov[v].nr]);
key->IML[v] = -1;
// special case: Nugget (part 1)
if (cov->cov==nugget) continue;
Allowed[TBM2] = Allowed[SpectralTBM] =Allowed[Hyperplane] =
key->cov[v].reduceddim <= 2;
if (key->cov[v].reduceddim == 3 && key->grid) {
Allowed[TBM2] = true;
Standard[TBM2pos] = TBM3;
Standard[TBM3pos] = TBM2;
}
for (i=0; i<nraw; i++) raw[i] = Forbidden;
int *implemented=cov->implemented, maxdim, CEbadlybehaved;
for (i=0; i<Nothing; i++)
allowed[i] = Allowed[i] && (implemented[i] == IMPLEMENTED ||
implemented[i] == NUM_APPROX);
if (DECISION_PARAM.stationary_only==DECISION_CASESPEC &&
!cov->variogram) allowed[CircEmbedIntrinsic] = false;
if (DECISION_PARAM.exactness != DECISION_FALSE &&
implemented[i]==NUM_APPROX)
allowed[TBM2] = false;
// multiplicative models
if (v<key->ncov-1 && key->cov[v].op) {
// nur die Reihenfolge fuer den ersten Faktor wird verwendet --
// der letzte Faktor wird gesetzt wie wenn nicht Faktor (spaeter ignoriert)
// alle anderen Faktoren werden als Faktoren gesetzt aber Methodenfolge
// wiederum ignoriert
// next_element_of_preference_list achtet auf die anderen Faktoren
for (i=0; i<Nothing; i++)
if (i!=TBM3 && i!=CircEmbed && i!=Direct) allowed[i] = false;
}
cov->info(key->cov[v].param, &maxdim, &CEbadlybehaved);
// printf("timespacedim %d %d %s\n", key->timespacedim, key->cov[v].reduceddim, cov->name);
assert(key->cov[v].reduceddim <= maxdim);
if (key->grid && CEbadlybehaved && DECISION_PARAM.exactness==DECISION_TRUE) {
// preference to Direct whenever possible?
DirectFst = DirectLst;
}
nr = 0;
raw[nr++] = DirectFst;
if ((((key->ncov==1 || DECISION_PARAM.exactness==DECISION_FALSE) &&
CEbadlybehaved) || CEbadlybehaved > 1) && key->grid) {
// key->ncov > 1 laesst hoffen, dass die andere Kovarianzfunktion
// sich freundlicher verhaelt (z.B. wenn Nugget), so dass die Matrix
// des Gesamtmodells sich gut bzgl. CE verhaelt.
if (key->cov[v].reduceddim==2 || CEbadlybehaved > 1) {
raw[nr++] = SpectralTBM;
if (DECISION_PARAM.exactness==DECISION_FALSE) raw[nr++] = TBM2;
}
if (key->cov[v].reduceddim==3 && CEbadlybehaved>1) raw[nr++] = TBM3;
}
if (key->grid && DECISION_PARAM.exactness==DECISION_FALSE &&
key->totalpoints * (1 << key->timespacedim) * 2*sizeof(double) > maxmem){
raw[nr++] = SpectralTBM;
raw[nr++] = TBM2;
raw[nr++] = TBM3;
}
for (i=0; i<nStandard; i++) raw[nr++] = Standard[i];
raw[nr++] = DirectLst;
for (i=0; i<nLastDefault; i++) raw[nr++] = LastDefault[i];
nml = 0;
for (i=0; i<nr; i++) {
if (raw[i]!=Forbidden && allowed[raw[i]]) {
key->ML[v][nml++] = raw[i];
allowed[raw[i]] = false;
}
}
key->NML[v] = nml;
switch (raw[0]) {
case Direct : anyFirstDirect = true; break;
case CircEmbed : anyFirstCE = true; break;
default: break;
}
}
// nugget, part 2
for (v=0; v<key->ncov; v++) if (CovList[key->cov[v].nr].cov==nugget){
nml = 0;
if (!key->anisotropy) {
// CE first considered, since nugget may have a positive effect on
// the CE simuation.
if (anyFirstCE) key->ML[v][nml++] = CircEmbed;
if (anyFirstDirect) key->ML[v][nml++] = Direct;
}
key->ML[v][nml++] = Nugget;
key->NML[v] = nml;
}
}
int next_element_of_method_list(key_type *key)
{
int v;
for (v=0; v<key->ncov; v++) {
if (key->cov[v].left) {
if (v>0 && key->cov[v-1].op) key->cov[v].method = key->cov[v-1].method;
else {
(key->IML[v])++;
// printf("next %d %d %d %d \n",
// v, key->IML[v], key->NML[v], key->ML[v][key->IML[v]]);
if (key->IML[v] >= key->NML[v]) return v + 1;
key->cov[v].method = key->ML[v][key->IML[v]];
}
}
}
return 0;
}
void InitSimulateRF(double *x, double *T,
int *spatialdim, /* spacial dim only ! */
int *lx, int *grid,
int *Time,
int *covnr, double *ParamList, int *nParam,
int *ncov, int *anisotropy, int *op,
int *method,
int *distr, /* still unused */
int *keyNr , int *error)
{
// keyNr : label where intermediate results are to be stored;
// can be chosen freely between 0 and MAXKEYS-1
key_type *key;
key = NULL;
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
*error=ERRORREGISTER;
goto ErrorHandling;
}
key = &(KEY[*keyNr]);
if (x==NULL) {*error=ERRORCOORDINATES; goto ErrorHandling;}
strcpy(ERROR_LOC, "");
*error =
internal_InitSimulateRF(x, T, *spatialdim,
*lx, (bool) *grid, (bool) *Time,
covnr, ParamList, *nParam,
*ncov, (bool) *anisotropy, op, method,
*distr, key,
GENERAL_NATURALSCALING, CovFct);
if (*error != NOERROR) goto ErrorHandling;
return;
ErrorHandling:
if (GENERAL_PRINTLEVEL>0) {
ErrorMessage(Nothing, *error);
PRINTF("\n");
}
}
int internal_InitSimulateRF(double *x, double *T,
int spatialdim, /* spatial dim only ! */
int lx, bool grid, bool Time,
int *covnr, double *ParamList, int nParam,
int ncov, bool anisotropy, int *op,
int *method,
int distr, /* still unused */
key_type *key,
int naturalscaling,
CovFctType covFct)
{ // grid : boolean
// NOTE: if grid and Time then the Time component is treated as if
// it was an additional space component
// On the other hand, SPACEISOTROPIC covariance function will
// always assume that the second calling component is a time component!
// printf("internal %d %d \n", covnr[0], covnr[1]);
bool user_defined, simple_definition;
int i, last_incompatible, d, v, M, err_occurred, error, timespacedim;
unsigned long totalBytes;
int method_used[(int) Nothing + 1];
// preference lists, distinguished by grid==true/false and dimension
// lists must end with Nothing!
SimulationType Merr=Nothing;
char errorloc_save[nErrorLoc];
strcpy(errorloc_save, ERROR_LOC);
timespacedim = spatialdim + (int) Time; // may not set to key->timespacedim,
// since key->timespacedim possibly by DeleteKeyNotTrend
// printf("timespacedim %d %d %d\n", timespacedim , spatialdim , (int) Time);
if (spatialdim<1 || timespacedim>MAXDIM){
error=ERRORDIM; goto ErrorHandling;}
user_defined = method[0] >= 0;
totalBytes = sizeof(double) * lx * spatialdim;
key->active &= !key->stop;
// also checked by R function `PrepareModel'
if ((error = CheckAndBuildCov(covnr, op, ncov, ParamList, nParam,
naturalscaling, Time, grid,
timespacedim, anisotropy,
key->cov, &(key->active))) != NOERROR)
goto ErrorHandling;
if (user_defined) {
for (v=0; v<ncov; v++) {
SimulationType meth;
covinfo_type *kc;
meth = (SimulationType) method[v];
// printf("method %d\n", method[v]);
kc = &(key->cov[v]);
if (CovList[kc->nr].cov == nugget && meth!=Nugget && ncov>1 &&
((meth != CircEmbed && meth != Direct) ||
kc->reduceddim < timespacedim)) {
if (GENERAL_PRINTLEVEL > 5)
PRINTF("method for nugget (nr. %d) set to nugget (orig=%s)",
v, METHODNAMES[meth]);
meth = Nugget;
}
key->active &= kc->method == meth;
key->cov[v].method = meth;
// printf("user %d %d %d \n", v, method[v], method);
}
}
if (key->active) {
assert(key->x[0]!=NULL);
key->active =
(covFct == key->covFct) &&
(key->grid == grid) &&
((key->grid && lx==3) || (!key->grid && key->totalpoints==lx)) &&
(key->spatialdim==spatialdim) &&
(key->distribution==distr) &&
(key->ncov==ncov) &&
(key->anisotropy==anisotropy) &&
(Time == key->Time &&
(!Time || !memcmp(key->T, T, sizeof(double) * 3))) &&
(!memcmp(key->x[0], x, totalBytes));
// if (GENERAL_PRINTLEVEL>5 && !key->active) {
// PRINTF("failure with previous initialisation at spatialdim=%d grid=%d distrib=%d ncov=%d aniso=%d time=%d key->T=%d op=%d x=%d\n",
// key->spatialdim==spatialdim, key->grid == grid,
// key->distribution==distr, key->ncov==ncov,
// key->anisotropy==anisotropy,
// Time == key->Time,
// !(Time) || ! memcmp(key->T, T, sizeof(double) * 3),
// !memcmp(key->x[0], x, totalBytes));
// }
}
if (key->active) {// detected that all relevant parameters agree,
// no initialization is necessary
// CHECK:
for (d = key->timespacedim; d<MAXDIM; d++) {assert(key->x[d]==NULL);}
if (GENERAL_PRINTLEVEL>=2) ErrorMessage(Nothing, USEOLDINIT);
return NOERROR; /* ***** END ***** */
} else {// ! key->active
DeleteKeyNotTrend(key);
// if not active, all the parameters have to be
// stored again
// save all the parameters if key has not been active
key->timespacedim = timespacedim;
key->spatialdim = spatialdim;
key->grid= grid;
key->anisotropy = anisotropy;
key->distribution=distr;
key->ncov=ncov;
// coordinates
if (key->x[0]!=NULL) free(key->x[0]);
if ((key->x[0]=(double*) malloc(totalBytes))==NULL){
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
memcpy(key->x[0], x, totalBytes);
for (d=1; d<key->spatialdim; d++) key->x[d]= &(key->x[0][d * lx]);
for (; d<MAXDIM; d++) key->x[d]=NULL;
//Time
if (key->Time = Time) {
if (!key->anisotropy) { error = ERRORTIMENOTANISO; goto ErrorHandling; }
memcpy(key->T, T, sizeof(double) * 3);
if (key->grid) {
key->x[key->spatialdim] = key->T;
}
}
}
key->covFct = covFct; // CovFct, CovFctTBM2, CovFctTBM2Num, CovFctTBM3
// not literally "total"param; there might be gaps in the sequence of the
// parameters (if not all seven KAPPA parameters are used)
key->totalparam = anisotropy
? ANISO + timespacedim * timespacedim
: SCALE + 1; // 0..SCALE
// Delete stored, intermediate result if there are any
// folgende Zeilen unschoen, da weiter vorne bereits DeleteKeyNotTrend
// aufgerufen wurde
// for (v=0; v<MAXCOV; v++) {
// if (key->meth[v].destruct!=NULL) {
// key->meth[v].destruct(&(key->meth[v].S));
// key->meth[v].destruct = NULL;
// }
// }
// minor settings, usually overtaken from previous run, if all the
// other parameters are identical
if (key->grid) {
if ((lx!=3) ||
(InternalGetGridSize(key->x,&key->timespacedim,key->length))) {
error=ERRORCOORDINATES; goto ErrorHandling;
}
for (key->spatialtotalpoints=1, d=0; d<key->spatialdim; d++) {
key->spatialtotalpoints *= key->length[d];
}
key->totalpoints = key->spatialtotalpoints;
if (key->Time) key->totalpoints *= key->length[key->spatialdim];
} else { // not grid
key->totalpoints = key->spatialtotalpoints = key->length[0]= lx;
if (key->Time) {
int Tdim; double *Tx[MAXDIM];
Tdim = 1;
Tx[0] = key->T;
if (InternalGetGridSize(Tx,&Tdim,&(key->length[key->spatialdim]))) {
error=ERRORCOORDINATES; goto ErrorHandling;
}
key->totalpoints *= key->length[key->spatialdim];
}
// just to say that considering these values does not make sense
for (d=1; d<key->spatialdim; d++) key->length[d]=0;
// correct value for higher dimensions (never used)
}
for (d=key->timespacedim; d<MAXDIM; d++) key->length[d] = (int) RF_NAN; // 1
// simple definition of the covariance function?
// (the only definition used up to version 1.0 of RandomFields)
simple_definition =
key->ncov==1 ||
(key->ncov==2 && !key->anisotropy && op[0]==0 &&
(key->cov[0].method==Nugget || key->cov[1].method==Nugget));
err_occurred = NOERROR;
error = ERROROUTOFMETHODLIST;
key->n_unimeth = 0;
last_incompatible = -1;
for (v=0; v<key->ncov; v++) key->cov[v].left=true;
// indicates which covariance function are left to be initialized
// method list must be treated separately since method is integer coded
// when passed to InitSimulate !
if (user_defined) {
//either all or none should be user defined
for (v=0; v<key->ncov; v++) {
ERRORMODELNUMBER = v;
covinfo_type *kc;
kc = &(key->cov[v]);
if (kc->method >= (int) Nothing) {
// printf("kc->method=%d %d %d \n", v, kc->method, (int) Nothing);
error=ERRORNOTDEFINED; goto ErrorHandling;}
if (CovList[kc->nr].type == ISOHYPERMODEL ||
CovList[kc->nr].type == ANISOHYPERMODEL) {
int i, nc = (int) kc->param[HYPERNR];
if (nc <= 0 || nc + v >= key->ncov) {
error=ERRORHYPERNR; goto ErrorHandling;}
for (i=v+1; i<nc+v; i++)
if (kc->method != key->cov[i].method) {
error=ERRORHYPERMETHOD; goto ErrorHandling;}
}
if (v>0 && op[v-1] && kc->method!=key->cov[v-1].method) {
error=ERRORMETHODMIX; goto ErrorHandling;}
if (DECISION_PARAM.stationary_only==DECISION_TRUE &&
kc->method==CircEmbedIntrinsic) {
error=ERRORMETHODEXCLUDED; goto ErrorHandling;}
key->NML[v] = 1;
key->IML[v] = -1;
key->ML[v][0] = kc->method;
// sonst: ML info wird nicht verwendet...
}
ERRORMODELNUMBER = -1;
} else {
for (v=0; v<key->ncov; v++) {
ERRORMODELNUMBER = v;
if (CovList[key->cov[v].nr].type == ISOHYPERMODEL ||
CovList[key->cov[v].nr].type == ANISOHYPERMODEL) {
error=ERRORHYPERMETHOD; goto ErrorHandling;
}
}
ERRORMODELNUMBER = -1;
init_method_list(key); // nicht weiter vorne,
// da es auf Daten von key zugreift, die weiter vorne noch nicht
// vorhanden sind
}
if (GENERAL_PRINTLEVEL>2) {
PRINTF("Lists of methods for the simulation\n===================================\n");
for (v=0; v<key->ncov; v++) {
PRINTF("%s: ", CovList[key->cov[v].nr].name);
for (i=0; i < key->NML[v]; i++) PRINTF("%s, ", METHODNAMES[key->ML[v][i]]);
PRINTF("\n");
}
}
while (error != NOERROR) {
int error_covnr, m;
error_covnr = NOERROR;
if (error != ERRORANYLEFT) {
error_covnr = next_element_of_method_list(key);
}
if (error_covnr) {
if (GENERAL_PRINTLEVEL>2)
PRINTF("covariance nr %d (%s) run out of method list.\n",
error_covnr, CovList[key->cov[error_covnr - 1].nr].name);
if (!user_defined) error = ERROROUTOFMETHODLIST;
break;
}
error = err_occurred = NOERROR;// necessary in case ncov=0, since then
// the following loop is not entered
// create list of methods used
for (i=0; i<=(int) Nothing; i++) method_used[i] = 0;
for (v=0; v<key->ncov; v++)
method_used[key->cov[v].method] += key->cov[v].left;
for (m=CircEmbed; m<=Nothing; m++) {
if (method_used[m]) {
insert_method(&last_incompatible, key, (SimulationType) m,
do_incompatiblemethod[m]!=NULL);
if (method_used[m]==1) {
v=0;
while(key->cov[v].method!=m || !key->cov[v].left) v++;
key->ML[v][key->IML[v]] = Nothing;
// es lohnt sich nicht, die Methoden durchzuleiern,
// die auf die Kovarianzfunktion bereits angewandt wurden und
// diese Kovarianzfunktion dabei die einzige war
}
}
}
for (M=0; M<key->n_unimeth; M++) {
if (!key->meth[M].unimeth_alreadyInUse) {
// unimeth_alreadyInUse: to avoid that method is called again
// when next while loop, which would leed to double use of key->meth[M]
bool left[MAXCOV];
methodvalue_type *meth;
meth = &key->meth[M];
meth->unimeth_alreadyInUse = true;
for (v=0; v<key->ncov; v++) left[v] = key->cov[v].left; // sicherung,
// da init_method ueberschreibt und *danach* entscheidet, ob
// die methode funktioniert -- ueber key->method ist es zu schwierig
// da eine methode mehrmals verwendet werden kann
if (GENERAL_PRINTLEVEL>=4) {
PRINTF("\n%sInit: method %d (%s); %d method instance(s) in use\n",
errorloc_save, M, METHODNAMES[meth->unimeth], key->n_unimeth);
if (GENERAL_PRINTLEVEL>=5)
for (v=0; v<key->ncov; v++)
PRINTF("%d (%s): left=%s, method=%s\n",
v, CovList[key->cov[v].nr].name,
key->cov[v].left ? "true" : "false",
METHODNAMES[key->cov[v].method]);
}
assert(meth->destruct==NULL && meth->S==NULL);
if (meth->unimeth > Nothing) {
error = ERRORUNKNOWNMETHOD;
goto ErrorHandling;
}
err_occurred = init_method[meth->unimeth](key, M);
if (GENERAL_PRINTLEVEL>=5) {
ErrorMessage(Nothing, err_occurred);
}
switch (err_occurred) {
case NOERROR_REPEAT:
insert_method(&last_incompatible, key, meth->unimeth,
do_incompatiblemethod[meth->unimeth]!=NULL); // ??!!
err_occurred = NOERROR;
break;
case NOERROR:
break;
case NOERROR_ENDOFLIST:
delete_method(&M, &last_incompatible, key);
err_occurred = NOERROR;
break;
default:
Merr = meth->unimeth; // do not shift after delete_method, since
// the latter changes the value of M !
delete_method(&M, &last_incompatible, key);
for (v=0; v<key->ncov; v++) {
covinfo_type *kc;
kc = &(key->cov[v]);
kc->left = left[v];
}
if (simple_definition && GENERAL_PRINTLEVEL>2) {
ErrorMessage(Merr,err_occurred);
}
error = err_occurred; // error may occur anywhere
// and final err_occured might be 0
}// switch
} // if not already in Use
} // for (M...)
if (error==NOERROR) {
for (v=0; v<key->ncov; v++) {
if (key->cov[v].left) {
error = ERRORANYLEFT;
Merr = Nothing;
break;
}
} // if (error == ERRORANYLEFT) continue;
}
} // while error!=NOERROR
if (error!=NOERROR && !simple_definition) {
int zaehler, vplus;
if (GENERAL_PRINTLEVEL>1) {
PRINTF("algorithm partially failed. Retrying. Last error has been:\n");
ErrorMessage(Merr, error);
}
zaehler = 0;
for (v=0; v<key->ncov; v++) if (key->cov[v].left) {
key->cov[v].method = Nothing;
zaehler++;
}
assert(zaehler > 0);
for (v=0; v<key->ncov; v+=vplus) {
bool left[MAXCOV];
covinfo_type *kc, *kcdummy;
vplus = 1;
kc = &(key->cov[v]);
if (kc->left) {
int w;
if (key->NML[v]==0) {
if (error != ERRORANYLEFT) key->n_unimeth++;
error = ERROROUTOFMETHODLIST;
goto ErrorHandling;
}
if (GENERAL_PRINTLEVEL>=5)
PRINTF("'%s' not initiated yet\n", CovList[kc->nr].name);
for (i=0; i<key->NML[v]; i++) {
if (key->ML[v][i] == Nothing) {
if (GENERAL_PRINTLEVEL>=6)
PRINTF("method '%s' (#%d of %d) jumped\n",
METHODNAMES[i], i, key->NML[v]);
continue;
}
if (GENERAL_PRINTLEVEL>=6)
PRINTF("trying '%s' (#%d of %d)\n",
METHODNAMES[key->ML[v][i]], i, key->NML[v]);
for (w=0; w<key->ncov; w++) left[w] = key->cov[w].left; // sicherung,
kcdummy = kc;
while (v + vplus < key->ncov && kcdummy->op) {
if (CovList[kcdummy->nr].type==ISOHYPERMODEL ||
CovList[kcdummy->nr].type==ANISOHYPERMODEL) {
// too complicated for the moment; should be improved
if (error != ERRORANYLEFT) key->n_unimeth++;
error = ERRORFAILED;
goto ErrorHandling;
}
// multiplication is left
kcdummy->method = key->ML[v][i];
kcdummy = &(key->cov[v + vplus]);
vplus++;
}
kcdummy->method = key->ML[v][i];
M = insert_method(&last_incompatible, key, kcdummy->method,
do_incompatiblemethod[kcdummy->method]!=NULL);
error = init_method[key->meth[M].unimeth](key, M);
if (GENERAL_PRINTLEVEL>=3) {
PRINTF("%s has return code %d;",
METHODNAMES[kcdummy->method], error);
ErrorMessage(Nothing, error);
}
if (error != NOERROR && error != NOERROR_REPEAT) {
Merr = kc->method;
delete_method(&M, &last_incompatible, key);
for (w=0; w<key->ncov; w++) {
kcdummy = &(key->cov[w]);
kc->left = left[w];
}
} else {
error = NOERROR;
break;
}
} // NML
if (error != NOERROR) {
if (key->ncov>2 && GENERAL_PRINTLEVEL>0)
PRINTF("The given model seems to be complex; try specifying explicitely the simulation methods, see help('GaussRF') and help('PrintMethodList') or even split up in subsequent calls of GaussRF\n");
if (GENERAL_PRINTLEVEL>3) {
PRINTF("All algorithm for %s failed. The last error has been:\n",
CovList[kc->nr].name);
ErrorMessage(Merr, error);
}
break;
}
} // left
} // for v < key->ncov
} // error & not simple definition
if (!(key->active = error==NOERROR)) {
if (error != ERRORANYLEFT) key->n_unimeth++;
goto ErrorHandling;
}
if (GENERAL_PRINTLEVEL>=4)
PRINTF("%sSuccessfully initialised.\n", errorloc_save);
key->compatible = last_incompatible <= 0; // if at most one incompatible
// method, then no intermediate results must be stored
strcpy(ERROR_LOC, errorloc_save);
return NOERROR;
ErrorHandling:
key->active = false;
if (GENERAL_PRINTLEVEL>3) {
ErrorMessage(Merr, error);
PRINTF("\n");
}
return error;
}
void AddTrend(int *keyNr, int *n, double *res, int *error) {
int i;
key_type *key;
*error = NOERROR;
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
*error=ERRORKEYNROUTOFRANGE; goto ErrorHandling;
}
key = &(KEY[*keyNr]);
if (!key->active) {*error=ERRORNOTINITIALIZED; goto ErrorHandling;}
// printf("addtrend %f \n", key->mean);
switch (key->TrendModus) {
case TREND_MEAN :
if (key->mean!=0.0) {
long int endfor= *n * key->totalpoints;
for(i=0; i<endfor; i++) res[i] += key->mean;
}
break;
case TREND_LINEAR :
case TREND_FCT :
case TREND_PARAM_FCT :
default: assert(false); // not programmed yet
}
return;
ErrorHandling:
if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing,*error);
return;
}
void DoPoissonRF(int *keyNr, int *pairs, int *n, double *res, int *error)
{
// not programmed yet
assert(false);
}
int internal_DoSimulateRF(key_type *key, int nn, double *orig_res) {
// does not assume that orig_res[...] = 0.0, but it is set
int error;
long i, m;
double *part_result, *res, realeach=0.0;
char back[]="\b\b\b\b\b\b\b\b\b\b\b", format[20], prozent[]="%";
int ni, digits, each=0;
res = orig_res;
part_result = NULL;
if (!key->active) {error=ERRORNOTINITIALIZED; goto ErrorHandling;}
if (nn>1 && GENERAL_PCH[0] != '\0') {
if (GENERAL_PCH[0] == '!') {
digits = (nn<900000000) ? 1 + (int) trunc(log((double) nn) / log(10.0)) : 9;
back[digits] = '\0';
each = (nn < 100) ? 1 : nn / 100;
sprintf(format, "%ss%s%dd", prozent, prozent, digits);
} else if (GENERAL_PCH[0] == '%') {
back[4] = '\0';
realeach = (double) nn / 100.0;
each = (nn < 100) ? 1 : (int) realeach;
sprintf(format, "%ss%s%dd%ss", prozent, prozent, 3, prozent);
} else each = 1;
} else each = nn + 1;
GetRNGstate();
if (!key->compatible) {
if ((part_result=(double *) malloc(sizeof(double) * key->totalpoints))
== NULL) {error=ERRORMEMORYALLOCATION; goto ErrorHandling;}
}
for (ni=1; ni<=nn; ni++, res += key->totalpoints) {
if (key->stop) {error=ERRORNOTINITIALIZED; goto ErrorHandling;}
if (ni % each == 0) {
if (GENERAL_PCH[0] == '!')
PRINTF(format, back, ni / each);
else if (GENERAL_PCH[0] == '%')
PRINTF(format, back, (int) (ni / realeach), prozent);
else PRINTF("%s", GENERAL_PCH);
}
for(m=0; m<key->n_unimeth; m++) {
methodvalue_type *meth;
meth = &(key->meth[m]);
if (GENERAL_PRINTLEVEL>=7)
PRINTF("DoSimu %d %s\n",m, METHODNAMES[meth->unimeth]);
if (meth->incompatible) {
if (m==0) {
if (GENERAL_PRINTLEVEL>=7) PRINTF("incomp\n");
do_incompatiblemethod[meth->unimeth](key, m, res);
} else {
if (GENERAL_PRINTLEVEL>=7) PRINTF("extra %d\n", key->compatible);
do_incompatiblemethod[meth->unimeth](key, m, part_result);
for (i=0; i<key->totalpoints; i++) res[i] += part_result[i];
}
} else {
do_compatiblemethod[meth->unimeth] (key, m, res);
}
} // for m
} // for n
PutRNGstate();
if (part_result!=NULL) free(part_result);
if (nn>1 && GENERAL_PCH[0] != '\0') {
if (GENERAL_PCH[0] == '!' || GENERAL_PCH[0] == '%') PRINTF("%s", back);
else PRINTF("\n");
}
return NOERROR;
ErrorHandling:
PutRNGstate();
if (part_result!=NULL) free(part_result);
key->active = false;
return error;
}
void DoSimulateRF(int *keyNr, int *n, int *pairs, double *res, int *error) {
// does not assume that orig_res[...] = 0.0, but it is set
int i, internal_n;
key_type *key=NULL;
*error=NOERROR;
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
*error=ERRORKEYNROUTOFRANGE; goto ErrorHandling;
}
key = &(KEY[*keyNr]);
internal_n = *n / (1 + (*pairs!=0));
if (!key->active) {
*error=ERRORNOTINITIALIZED; goto ErrorHandling;
}
if (key->n_unimeth == 0 || !key->meth[0].incompatible) {
long total =internal_n * key->totalpoints ;
for (i=0; i<total; i++)res[i] = 0.0;
}
if ((*error = internal_DoSimulateRF(key, internal_n, res)) != NOERROR)
goto ErrorHandling;
if (*pairs) {
double* res_pair;
long endfor=key->totalpoints * *n / 2;
res_pair = res + endfor;
for (i=0; i<endfor; i++) res_pair[i] = -res[i];
}
AddTrend(keyNr, n, res, error);
if (*error) goto ErrorHandling;
if (!(key->active = GENERAL_STORING))
DeleteKey(keyNr);
return;
ErrorHandling:
if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing, *error);
if (key!=NULL) {
key->active = false;
DeleteKey(keyNr);
}
return;
}
int InternalSimulate(key_type *key, double *orig_res) {
if (key->n_unimeth == 0 || !key->meth[0].incompatible) {
long i, total = key->totalpoints;
double *sres = orig_res;
for (i=0; i<total; i++) sres[i] = 0.0;
}
return internal_DoSimulateRF(key, 1, orig_res);
}
SEXP GetExtModelInfo(SEXP keynr) {
int knr;
knr = INTEGER(keynr)[0];
if (knr>=0) {
if (knr < MAXKEYS) {
key_type *key;
key = &(KEY[knr]);
return GetModelInfo(key->cov, key->ncov, key->totalparam,
key->totalpoints);
}
} else if (knr == -1) {
if (user_ncov>0)
return GetModelInfo(user_keycov, user_ncov, user_anisotropy ?
ANISO + user_dim * user_dim : SCALE, 0);
} else if (knr == -2) {
if (Unchecked_ncov > -1)
return GetModelInfo(Unchecked_keycov, Unchecked_ncov,
Unchecked_anisotropy
? ANISO + Unchecked_dim * Unchecked_dim : SCALE, 0);
}
return allocVector(VECSXP, 0);
}