https://github.com/cran/RandomFields
Tip revision: 270c4ba0da8194f2f81416d85bae1ccbd3ed4933 authored by Martin Schlather on 08 August 1977, 00:00:00 UTC
version 1.0.4
version 1.0.4
Tip revision: 270c4ba
MPP.cc
/*
Authors
Martin Schlather, Martin.Schlather@uni-bayreuth.de
*********** this function is still under construction *********
Copyright (C) 2001 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 "MathlibCC.h"
#include <R_ext/Mathlib.h>
#include <assert.h>
#include "RFsimu.h"
#include "MPPFcts.h"
//#include "RFCovFcts.h"
#include <assert.h>
Real MPP_RADIUS =0.0;
Real MPP_APPROXZERO =0.001;
Real ADDMPP_REALISATIONS = 100.0; /* logically an integer, implementation allows for
positive real ! */
void SetMPP(int *action, Real *approxzero, Real *realisations, Real *radius)
{
switch(*action) {
case 0 :
if (*approxzero<=0.0) {
if (GENERAL_PRINTLEVEL>0) PRINTF("\nERROR! `approx.zero' not a positive number. Ignored.\n");
} else MPP_APPROXZERO = *approxzero;
if (*realisations<=0.0) {
if (GENERAL_PRINTLEVEL>0) PRINTF("\nERROR! `realisation' not a positive number. Ignored.\n");
} else ADDMPP_REALISATIONS = *realisations;
MPP_RADIUS = *radius;
if ((MPP_RADIUS>0.0) && (GENERAL_PRINTLEVEL>1))
PRINTF("\nWARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect.\n");
break;
case 1 :
*approxzero = MPP_APPROXZERO;
*realisations = ADDMPP_REALISATIONS;
*radius = MPP_RADIUS;
if (GetNotPrint) break;
case 2 :
PRINTF("\nmarked point process methods\n============================\napproxzero=%e,\nrealisations=%e\n",
MPP_APPROXZERO,ADDMPP_REALISATIONS);
break;
default : if (GENERAL_PRINTLEVEL>0) PRINTF("unknown action\n");
}
}
void printMPPinfo(int *keyNr)
{
mpp_storage *bs;
if (((KEY[*keyNr].method!=AdditiveMpp) && (KEY[*keyNr].method!=MaxMpp))
|| (!KEY[*keyNr].active)) {
PRINTF("no information available\n");
return;
}
assert(KEY[*keyNr].S!=NULL);
bs = (mpp_storage*) KEY[*keyNr].S;
PRINTF("integral=%e\n",bs->integral);
PRINTF("integralsq=%e\n",bs->integralsq);
PRINTF("integralpos=%e\n",bs->integralpos);
PRINTF("eff.Radius=%e\n",bs->effectiveRadius);
PRINTF("eff.area=%e\n",bs->effectivearea);
PRINTF("maxheight=%e\n",bs->maxheight);
PRINTF("add.Radius=%e\n",bs->addradius);
}
void GetMPPInfo(int *keyNr, Real *integral, Real *integralsq,
Real *integralpos, Real *maxheight)
{
mpp_storage *bs;
if (((KEY[*keyNr].method!=AdditiveMpp) && (KEY[*keyNr].method!=MaxMpp))
|| (!KEY[*keyNr].active)) {
*integral = *integralsq = *integralpos = *maxheight = RF_NAN;
return;
}
assert(KEY[*keyNr].S!=NULL);
bs = (mpp_storage*) KEY[*keyNr].S;
*integral = bs->integral;
*integralsq = bs->integralsq;
*integralpos = bs->integralpos;
*maxheight = bs->maxheight;
}
void mpp_destruct(void **S) {
// mpp_storage in RFsimu.h !!
if (*S!=NULL) {
free(*S);
*S = NULL;
}
}
int init_mpp(key_type * key)
{
int error, d, i;
Real max[MAXDIM];
mpp_storage *bs;
assert(key->active==false);
assert((key->covnr>=0) &&(key->covnr<currentNrCov));
//assert(key->cov->cov!=NULL); here is the first !!
// where cov is unknown, but should be simulated
assert(key->param[VARIANCE]>=0.0);
assert(key->param[SILL]>=key->param[VARIANCE]);
assert( fabs(key->param[SCALE]*key->param[INVSCALE]-1.0) < EPSILON);
assert((key->S==NULL) && (key->destruct==NULL));
assert((key->cov->add_mpp_scl!=NULL) && (key->cov->add_mpp_rnd!=NULL));
if ((key->S=malloc(sizeof(mpp_storage)))==0){
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
key->destruct = mpp_destruct;
// determine the minimum area where random field values are to be generated
bs = (mpp_storage*) key->S;
if (key->grid) {
for (d=0; d<key->dim; d++) {
bs->min[d] = key->x[d][XSTART];
bs->length[d]= key->x[d][XSTEP] * (Real) (key->length[d] - 1);
}
} else {
for (d=0;d<MAXDIM;d++) {bs->min[d]=RF_INF; max[d]=RF_NEGINF;}
for (i=0;i<key->totallength;i++) {
for (d=0;d<key->dim;d++) {
if (key->x[d][i]<bs->min[d]) bs->min[d]=key->x[d][i];
if (key->x[d][i]>max[d]) max[d]=key->x[d][i];
}
}
for (d=0;d<MAXDIM;d++) { bs->length[d] = max[d] - bs->min[d]; }
}
bs->dim = key->dim;
bs -> addradius = MPP_RADIUS;
if ((MPP_RADIUS>0.0) && (GENERAL_PRINTLEVEL>=2))
PRINTF("Note: area has been enlarged by fixed value (%f)",bs->addradius);
key->cov->add_mpp_scl(key);
return NOERROR;
ErrorHandling:
if (key->destruct!=NULL) key->destruct(&key->S);
key->destruct=NULL;
key->active = false;
return error;
}
void do_addmpp(key_type *key, Real *res END_WITH_GSLRNG)
{
int i, d, dim;
Real lambda, min[MAXDIM], max[MAXDIM], factor, average, poisson;
long segment[MAXDIM+1];
mpp_storage *bs;
mppmodel model;
MPPRandom MppFct;
assert(key->active);
#ifdef RF_GSL
assert(RANDOM!=NULL);
#endif
bs = (mpp_storage*) key->S;
MppFct = key->cov->add_mpp_rnd;
dim = key->dim;
switch (key->distribution) {
case DISTR_GAUSS :
lambda = ADDMPP_REALISATIONS * bs->effectivearea;
factor =
sqrt(key->param[VARIANCE] / (ADDMPP_REALISATIONS * bs->integralsq));
average = ADDMPP_REALISATIONS * bs->integral;
break;
case DISTR_POISSON :
lambda = key->param[VARIANCE] * bs->effectivearea;
factor = average = 0 /* replace by true value !! */; assert(false);
if(key->param[VARIANCE]!=key->param[MEAN]) {
// if (GENERAL_PRINTLEVEL>0)
ErrorMessage(AdditiveMpp,ERRORVARMEAN); // print it always
assert(false); // for debugging -- delete this later on
// finally, the value of MEAN is ignored, and that of VARIANCE
// is used for simulating the random field
}
break;
default : assert(false);
}
segment[0] = 1;
for (d=0; d<dim; d++)
segment[d+1] = segment[d] * key->length[d];
for (i=0; i<key->totallength; i++) res[i]=0.0;
for(poisson = -log(1.0 - UNIFORM_RANDOM);
poisson < lambda;
poisson -= log(1.0 - UNIFORM_RANDOM))
{
MppFct( (mpp_storage*) key->S, min, max, &model);
if (key->grid) {
int start[MAXDIM], end[MAXDIM], resindex, index[MAXDIM],
segmentdelta[MAXDIM], dimM1;
Real coord[MAXDIM], startcoord[MAXDIM];
// determine rectangle of indices, where the mpp function
// is different from zero
for (d=0; d<dim; d++) {
if (min[d]< key->x[d][XSTART]) {start[d]=0;}
else start[d] = (int) ( (min[d] - key->x[d][XSTART]) / key->x[d][XSTEP]);
// "end[d] = 1 + *" since inequalities are "* < end[d]"
// not "* <= end[d]" !
end[d] = 1 + (int) ((max[d] - key->x[d][XSTART]) / key->x[d][XSTEP]);
if (end[d] > key->length[d]) end[d] = key->length[d];
}
// prepare coord starting value and the segment vectors for res
resindex = 0;
for (d=0; d<dim; d++) {
index[d]=start[d];
segmentdelta[d] = segment[d] * (end[d]-start[d]);
resindex += segment[d] * start[d];
coord[d] = startcoord[d] =
key->x[d][XSTART] + (Real)start[d] * key->x[d][XSTEP];
}
// add mpp function to res
dimM1 = dim -1;
d = 0; // only needed for assert(d<dim)
while (index[dimM1]<end[dimM1]) {
assert(d<dim);
res[resindex] += model(coord);
d=0;
index[d]++;
coord[d]+=key->x[d][XSTEP];
resindex++;
while (index[d] >= end[d]) {
if (d>=dimM1) break; // below not possible
// as dim==1 will fail!
index[d]=start[d];
coord[d]=startcoord[d];
resindex -= segmentdelta[d];
d++; // if (d>=dim) break;
index[d]++;
coord[d]+=key->x[d][XSTEP];
resindex += segment[d];
}
}
} else {
Real y[MAXDIM];
// the following algorithm can greatly be improved !
// but for ease, just the stupid algorithm
for (i=0; i<key->totallength; i++) {
for (d=0; d<dim; d++) { y[d] = key->x[d][i]; }
res[i] += model(y);
}
}
}
if (key->distribution==DISTR_GAUSS) {
for (i=0; i<key->totallength;i++)
res[i] = factor * (res[i] - average);
}
}