https://github.com/cran/RandomFields
Tip revision: 4877e49dad8ee6b04e79289f69ff7f2186f11506 authored by Martin Schlather on 20 January 2012, 00:00:00 UTC
version 2.0.54
version 2.0.54
Tip revision: 4877e49
nugget.cc
/*
Authors
Martin Schlather, martin.schlather@math.uni-goettingen.de
all around the nugget effect -- needs special treatment
Copyright (C) 2001 -- 2011 Martin Schlather,
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <R_ext/Lapack.h>
#include "RF.h"
#include "primitive.h"
void SetParamNugget(int *action, double *nuggettol, int *nuggetmeth)
{
nugget_param *gp = &(GLOBAL.nugget);
if (*action) {
gp->tol = *nuggettol;
if (gp->tol < 0) {
if (PL>=1)
PRINTF("negative tolerance for distance in nugget covariance not allowed; set to zero");
gp->tol = 0.0;
}
gp->meth = (bool) *nuggetmeth;
} else {
*nuggettol = gp->tol;
*nuggetmeth = (int) gp->meth;
}
}
void nugget_destruct(void ** S)
{
if (*S != NULL) {
nugget_storage *x;
x = *((nugget_storage**)S);
if (x->pos!=NULL) free(x->pos);
if (x->red_field!=NULL) free(x->red_field);
free(*S);
*S = NULL;
}
}
bool equal(int i, int j, double *X, int dim)
{
double *x, *y, v, dummy, dist;
cov_model cov;
int d;
x = X + i * dim;
y = X + j * dim;
for (dist=0.0, d=0; d<dim; d++) {
dummy = x[d]-y[d];
dist += dummy * dummy;
}
dist = sqrt(dist);
cov.vdim = 1;
nugget(&dist, &cov, &v);
return v==1.0;
}
// uses global RANDOM !!!
int init_nugget(method_type *meth){
cov_model *cov = meth->cov;
location_type *loc = meth->loc;
nugget_storage *s;
globalparam *gp = meth->gp;
nugget_param* lp = &(gp->nugget);
int err, i, vdim,
PL = gp->general.printlevel,
origdim = loc->timespacedim,
dim = cov->tsdim,
dimSq = origdim * origdim;
SET_DESTRUCT(nugget_destruct);
if ((meth->S=malloc(sizeof(mpp_storage)))==0){
err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
s = (nugget_storage*) meth->S;
s->pos = NULL;
s->red_field = NULL;
vdim = cov->calling->vdim;
if (vdim <= 0) {
err=ERRORMULTIMISMATCH; goto ErrorHandling;
}
cov->vdim = vdim;
// PrintMethodInfo(meth);
// PrintMethodInfo(meth);
if ((s->simple = origdim == dim)) {
double value[MAXNUGGDIM], ivalue[MAXNUGGDIM], dummy[5 * MAXNUGGDIM], *A;
int j,k, m,
ndummy = 5 * origdim, err;
char No = 'N';
// Upper = 'U';
if (dim > MAXNUGGDIM) {
ERR("dim larger then MAXNUGGDIM");
}
if (meth->cproj != NULL ){
ERR("projections are not programmed for nugget yet");
}
if (meth->caniso == NULL) {
s->simple = true;
} else {
A = (double*) malloc(dimSq * sizeof(double));
// memcpy(A, meth->caniso, sizeof(double) * origdim * origdim);
// because of numerical errors
// A = canio^T * caniso
for (k=i=0; i<dimSq; i+=origdim) {
for (j=0; j<dimSq; j+=origdim,k++) {
A[k] = 0.0;
for (m=0; m<origdim; m++) {
A[k] += meth->caniso[i+m] * meth->caniso[j+m];
}
}
}
// printf("%e %e %e %e\n", A[0], A[1], A[2], A[3]);
// A[0] = A[1] = A[2] = A[3] = 1.0;
F77_NAME(dgeev)(&No, &No, &origdim, A, &origdim,
value, ivalue, NULL, &origdim, NULL, &origdim,
dummy, &ndummy, &err);
if (err != 0) {
free(A);
ERR("dgeev failed in nugget.cc");
}
for (i=0; i<origdim; i++) {
// printf("simple %d %e %e %e %d\n",
// i, value[i], ivalue[i],EIGENVALUE_EPS,
// fabs(value[i]) + fabs(ivalue[i]) > EIGENVALUE_EPS);
if (!(s->simple = fabs(value[i]) + fabs(ivalue[i]) > EIGENVALUE_EPS))
break;
}
free(A);
}
}
s->simugrid = loc->grid && (meth->type == TypeDiag || meth->type == TypeIso);
// printf("%d\n", s->simple); assert(false);
if (!s->simple) {
int *pos, oldpos;
if (lp->tol==0.0 && PL>=1)
PRINTF("\nThe anisotropy matrix does not have full rank and RFparameters()$nugget.tol equals 0. From a theoretical point of view that's fine, but the simulations will probably be odd. Is this really what you want?\n");
if (s->simugrid) {
// TypeIso und nicht simple genau dann wenn identisch 0
int d, dimP1=dim + 1;
s->prod_dim[0] = 1;
if (meth->cproj != NULL ){
assert(false);
// only projection possible here, but not programmed yet
} else {
for (d=i=0; d<dim; d++, i+=dimP1) {
s->reduced_dim[d] =
fabs(meth->caniso[i]) < lp->tol ? 1 : loc->length[d];
s->prod_dim[d + 1] = s->prod_dim[d] * s->reduced_dim[d];
// printf("%d %d %d %d\n", i, s->prod_dim[0], s->prod_dim[1],
// s->prod_dim[2]);
}
}
if ((s->red_field=(res_type *) malloc(sizeof(res_type) * cov->vdim *
s->prod_dim[dim]))
== NULL){
err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
} else {
if ((pos = (int*) malloc(sizeof(int) * loc->totalpoints))==0) {
err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
Transform2NoGridExt(meth, false, true, &(meth->space), &(meth->sptime));
// assert(false);
ordering(meth->sptime, loc->totalpoints, dim, pos);
oldpos = pos[0];
for (i=1 /* ! */; i<loc->totalpoints; i++) {
if (equal(oldpos, pos[i], meth->sptime, cov->tsdim))
pos[i]= -1 - pos[i];
else oldpos=pos[i];
}
s->pos = pos;
}
}
// s->sqrtnugget = sqrt(meth->cvar);
// printf("%f %f\n", meth->cvar, s->sqrtnugget); assert(false);
return NOERROR;
ErrorHandling:
return err;
}
void do_nugget(method_type *meth, res_type *res) {
cov_model *cov = meth->cov;
location_type *loc = meth->loc;
/// double sqrtnugget;
nugget_storage* s;
long nx, endfor;
int v,
vdim = cov->vdim;
assert(meth->S != NULL);
s = (nugget_storage*) meth->S;
// sqrtnugget = s->sqrtnugget;
if (s->simple) {
for (nx=0, endfor=loc->totalpoints * vdim; nx<endfor; nx++) {
// if (nx==0) printf("nug %f ", res[nx]);
res[nx] += (res_type) GAUSS_RANDOM(1.0);
// if (nx==0)printf("%f \n", (double) res[nx]);
}
} else {
if (s->simugrid) {
int d, i, dim, dimM1, index[MAXNUGGDIM], *red_dim, *prod_dim, endfor;
long totpnts, idx;
res_type *field;
totpnts = loc->totalpoints;
dim = cov->tsdim;
dimM1 = dim -1;
field = s->red_field;
red_dim = s->reduced_dim;
prod_dim = s->prod_dim;
//PrintMethodInfo(meth);
// printf("%d %d %d %d\n", dim, s->prod_dim[0], s->prod_dim[1],
// s->prod_dim[2]);
endfor = vdim * s->prod_dim[dim-1];
for (i=0; i<endfor; i++) {
field[i] = (res_type) GAUSS_RANDOM(1.0);
}
for (d=0; d<dim; index[d++] = 0);
for (i=0; i<totpnts; i++) {
for(idx=d=0; d<dim; d++) idx += (index[d] % red_dim[d]) * prod_dim[d];
//printf("%d %d %f\n", i, idx, (double) field[idx]);
for (v=0; v<vdim; v++) {
res[i + v] += field[idx + v];
}
d = 0;
(index[d])++;
while (d < dimM1 && index[d] >= loc->length[d]) {
assert(d<dim); // assert(d>=0);
index[d] = 0;
d++;
(index[d])++;
}
}
} else {
int p;
double *dummy = (double*) malloc(sizeof(double) * vdim);
assert(s->pos[0]>=0);
for (v=0; v<vdim; v++) dummy[v] = RF_NAN; // just to avoid warnings
// from the compiler
for (nx=0; nx<loc->totalpoints; nx++) {
if ((p = s->pos[nx]) < 0) p = -1 - p; // if p<0 then take old variable
// and -p-1 is the true index
else {
for (v=0; v<vdim; v++)
dummy[v] = GAUSS_RANDOM(1.0);
}
for (v=0; v<vdim; v++)
res[p + v] += dummy[v];
}
free(dummy);
}
}
}