/*
Authors
Martin Schlather, martin.schlather@cu.lu
all around the nugget effect -- needs special treatment
Copyright (C) 2001 -- 2004 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 <assert.h>
#include "RFsimu.h"
static Real NUGGET_TOL_SINGULAR=1e-12;
static Real NUGGET_TOL_POINTS=1e-12;
typedef struct nugget_storage {
double sqrtnugget;
bool simple;
int *pos;
} nugget_storage;
void nugget_destruct(void ** S)
{
if (*S != NULL) {
nugget_storage *x;
x = *((nugget_storage**)S);
if (x->pos!=NULL) free(x->pos);
free(*S);
*S = NULL;
}
}
/* nugget effect model */
Real nugget(Real *x, Real *p, int dim){
if (*x==0.0) return 1.0; return 0.0;
}
Real Scalenugget(Real *p, int scaling) { return 1.0; }//or better 0.0 => error?
int is_singular(Real *C, int dim, bool *singular)
{
double D[MAXDIM];
int Xerror, d;
if ((Xerror=eigenvalues(C, dim, D))!=NOERROR) goto ErrorHandling;
*singular = false;
for (d=0; d<dim; d++)
if (fabs(D[d])<NUGGET_TOL_SINGULAR) {*singular=true; break;}
return 0;
ErrorHandling:
return Xerror;
}
bool equal(int i, int j, Real *ORDERD, int ORDERDIM)
{
Real *x, *y;
int d;
x = ORDERD + i * ORDERDIM;
y = ORDERD + j * ORDERDIM;
for(d=0; d<ORDERDIM; d++) {
if(fabs(x[d]-y[d]) > fabs(x[d]>y[d] ? x[d] : y[d]) * NUGGET_TOL_POINTS) {
return false;
}
}
return true;
}
// uses global RANDOM !!!
int init_nugget(key_type *key, int m){
int Xerror;
nugget_storage *s;
param_type param;
Real nugget_effect, quotient[MAXCOV];
unsigned short int actcov;
int i,covnr[MAXCOV], nonzero_pos;
int multiply[MAXCOV], TrueDim;
bool singular, no_last_comp;
Real *xx;
xx=NULL;
nonzero_pos = -1;
SET_DESTRUCT(nugget_destruct);
FIRSTCHECK_COV_ANISO(Nugget, cov, param, false);
if ((key->S[m]=malloc(sizeof(mpp_storage)))==0){
Xerror=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
s = (nugget_storage*) key->S[m];
s->pos = NULL;
if (key->anisotropy) {
if ((Xerror = is_singular(&(param[0][ANISO]), key->timespacedim, &singular))
!=NOERROR)
goto ErrorHandling;
s->simple = !singular;
if (singular) {
int *pos, oldpos, start_param[MAXDIM], index_dim[MAXDIM];
GetTrueDim(key->anisotropy, key->timespacedim, param[0],
&TrueDim, &no_last_comp, start_param, index_dim);
if ((Xerror=Transform2NoGrid(key, param[0], TrueDim, start_param, &(xx)))
!=NOERROR)
goto ErrorHandling;
if ((pos = (int*) malloc(sizeof(int) * key->totalpoints))==0){
Xerror=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
ordering(xx, key->totalpoints, TrueDim, pos);
oldpos = pos[0];
for (i=1 /* ! */; i<key->totalpoints; i++) {
if (equal(oldpos, pos[i], xx, TrueDim)) pos[i]= -1 - pos[i];
else oldpos=pos[i];
}
s->pos = pos;
}
} else {
s->simple = true;
}
for (i=0; i<actcov; i++) assert(CovList[covnr[i]].cov==nugget);
nugget_effect = CovFct(ZERO,1, covnr, multiply, param, actcov, false);
s->sqrtnugget = sqrt(nugget_effect);
if (xx!=NULL) free(NULL);
if (key->anisotropy) return NOERROR_REPEAT;
else return 0;
ErrorHandling:
if (xx!=NULL) free(NULL);
return Xerror;
}
void do_nugget(key_type *key, bool add, int m, Real *res ) {
Real sqrtnugget;
nugget_storage* s;
long nx, endfor;
#define RESULT(OP)\
for (nx=0,endfor=key->totalpoints; nx<endfor; nx++)\
res[nx] OP GAUSS_RANDOM(sqrtnugget);\
s = (nugget_storage*) key->S[m];
sqrtnugget = s->sqrtnugget;
if (s->simple) {
if (add) { RESULT(+=) } else { RESULT(=) }
} else {
int p;
Real dummy;
assert(s->pos[0]>=0);
dummy = RF_NAN;
for (nx=0; nx<key->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 dummy = GAUSS_RANDOM(sqrtnugget);
if (add) res[p] += dummy;
else res[p] = dummy;
}
}
}
void rangenugget(int dim, int *index, Real* range){
*index = -1;
}